C++ MultiArray类代码示例

void testRoi() {
    using namespace vigra;
    typedef ChannelSelector<4, float>::V V;
    typedef vigra::TinyVector<int, 4> V4;
    
    int ch = 0;
   
    MultiArray<4, float> data(vigra::TinyVector<int, 4>(10,20,30,2));
    FillRandom<float, float*>::fillRandom(data.data(), data.data()+data.size());
    {
        HDF5File f("test.h5", HDF5File::Open);
        f.write("test", data);
    }
    
    SourceHDF5<4, float> source("test.h5", "test");
    source.setRoi(Roi<4>(V4(1,3,5,0), V4(7,9,30,2)));
    
    SinkHDF5<3, float> sink("channel.h5", "channel");
    
    ChannelSelector<4, float> cs(&source, V(10,10,10));
    
    cs.run(3, ch, &sink);

    HDF5File f("channel.h5", HDF5File::Open);
    MultiArray<3, float> r;
    f.readAndResize("channel", r);

    MultiArrayView<3, float> shouldResult = data.bind<3>(ch).subarray(V(1,3,5), V(7,9,30));
    
    shouldEqualSequence(r.begin(), r.end(), shouldResult.begin());
}

Fields FieldAlgorithms::fieldsByLaplasian(MultiArray<2, float> & image) 
{
    MultiArray<2, float> valley(image.shape());
    Pyramid pyramid(image);
    for (int i = 3; i > 0; i--)
    {
        MultiArray<2, float> resized(pyramid.get(i));
        MultiArray<2, float> tmpArr(resized.shape());
        laplacianOfGaussianMultiArray(resized, tmpArr, 1.0);
        valley += pyramid.toOriginalSize(tmpArr);
    }
    MultiArray<2, float> peak = valley * -1;
    //remove DC component
    float thrhld = valley[argMax(valley)] * 0.3;
    threshold(valley, valley, thrhld);
    thrhld = peak[argMax(peak)] * 0.3;
    threshold(peak, peak, thrhld);
    

    //Edge as Gradients
    MultiArray<2, float> edgeField(image.shape());
    gaussianGradientMagnitude(image, edgeField, 1.0);

    //Localize
    std::vector<Shape2> valleyLocals(localizePOI(image));

    //Result as Field Object
    Fields fields(valley, valleyLocals, peak, edgeField, image);
    //A heuristic initialization of the localization, as a priori known position
    Shape2 nextToIris = Shape2(image.width() / 2, image.height() / 2);
    fields.specializedIrisValley = localizeByFollowingLocalMaxima(image, nextToIris);
    return fields;
}

void test() {
    using namespace vigra;
    typedef ChannelSelector<4, float>::V V;
   
    for(int ch=0; ch<=1; ++ch) {
        std::cout << "* channel = " << ch << std::endl;
    
        MultiArray<4, float> data(vigra::TinyVector<int, 4>(10,20,30,2));
        FillRandom<float, float*>::fillRandom(data.data(), data.data()+data.size());
        {
            HDF5File f("test.h5", HDF5File::Open);
            f.write("test", data);
        }
        
        SourceHDF5<4, float> source("test.h5", "test");
        SinkHDF5<3, float> sink("channel.h5", "channel");
        sink.setBlockShape(V(10,10,10));
        
        ChannelSelector<4, float> cs(&source, V(10,10,10));
        
        cs.run(3, ch, &sink);
    
        HDF5File f("channel.h5", HDF5File::Open);
        MultiArray<3, float> r;
        f.readAndResize("channel", r);
    
        MultiArrayView<3, float> shouldResult = data.bind<3>(ch);
        
        shouldEqualSequence(r.begin(), r.end(), shouldResult.begin());
    } //loop through channels
}

std::vector<Shape2> FieldAlgorithms::localizePOIExample(const MultiArray<2, float> &image,  MultiArray<2, RGBValue<UInt8> > &rgb_array)
{
    std::vector<Shape2> pois(30);
    float thld = image[argMax(image)] * 0.2;
    for (int i = 0; i < pois.size(); i++)
    {
        int v1 = std::rand() % image.size() -1;         // v1 in the range 0 to image.size()
        Shape2 poi(image.scanOrderIndexToCoordinate(v1));
        poi = localizeByFollowingLocalMaxima(image, poi);
        if (image[poi] > thld)
        {
            pois[i] = poi;
            MultiArrayView<2, RGBValue<UInt8> > markAsStep(rgb_array.subarray(Shape2(poi[0] - 5, poi[1] -5), Shape2(poi[0] +5, poi[1] +5)));
            for (RGBValue<UInt8> & val : markAsStep)
            {
                val.setRed(200);
            }
        }
        else
        {
            i--;
        }
    }
    return pois;

}

void test() {
    using namespace vigra;
    typedef Thresholding<3, float>::V V;

    MultiArray<3, float> data(V(24,33,40));
    FillRandom<float, float*>::fillRandom(data.data(), data.data()+data.size());
    {
        HDF5File f("test.h5", HDF5File::Open);
        f.write("test", data);
    }

    SourceHDF5<3, float> source("test.h5", "test");
    SinkHDF5<3, vigra::UInt8> sink("thresh.h5", "thresh");
    sink.setBlockShape(V(10,10,10));

    Thresholding<3, float> bs(&source, V(6,4,7));

    bs.run(0.5, 0, 1, &sink);

    HDF5File f("thresh.h5", HDF5File::Open);
    MultiArray<3, UInt8> r;
    f.readAndResize("thresh", r);

    MultiArray<3, UInt8> t(data.shape());
    transformMultiArray(srcMultiArrayRange(data), destMultiArray(t), Threshold<float, UInt8>(-std::numeric_limits<float>::max(), 0.5, 1, 0));
    shouldEqual(t.shape(), r.shape());
    shouldEqualSequence(r.begin(), r.end(), t.begin());
}

MultiArray<complex<float>, 4> reconMGE(Agilent::FID &fid) {
    int nx = fid.procpar().realValue("np") / 2;
    int ny = fid.procpar().realValue("nv");
    int nz = fid.procpar().realValue("nv2");
    int narray = fid.procpar().realValue("arraydim");
    int ne = fid.procpar().realValue("ne");

    MultiArray<complex<float>, 4> vols({nx, ny, nz, narray*ne});
    int vol = 0;
    if (verbose) cout << "Reading MGE fid" << endl;
    for (int a = 0; a < narray; a++) {
        if (verbose) cout << "Reading block " << a << endl;
        shared_ptr<vector<complex<float>>> block = make_shared<vector<complex<float>>>();
        *block = fid.readBlock(a);
        int e_offset = 0;
        for (int e = 0; e < ne; e++) {
            if (verbose)  cout << "Reading echo " << e << endl;
            MultiArray<complex<float>, 3> this_vol({nx, ny, nz}, block, {1,ne*nx,ne*nx*ny}, e_offset);
            MultiArray<complex<float>, 3> slice = vols.slice<3>({0,0,0,vol},{-1,-1,-1,0});

            auto it1 = this_vol.begin();
            auto it2 = slice.begin();
            while (it1 != this_vol.end()) {
                *it2++ = *it1++;
            }
            vol++;
            e_offset += nx;
        }
    }
    return vols;
}

Shape2 FieldAlgorithms::localizeByFollowingLocalMaxima(const MultiArray<2, float> &image, Shape2 current)
{

    //Open viewBox of image with center at current
    int upperLeftX = current[0] - ((image.width() / 10) / 2);
    upperLeftX = upperLeftX > -1 ? upperLeftX : 0;
    int upperLeftY = current[1] - ((image.height() / 10) / 2);
    upperLeftY = upperLeftY > -1 ? upperLeftY : 0;
    Shape2 upperLeft(upperLeftX, upperLeftY);
     

    int lowerRightX = current[0] + ((image.width() / 10) / 2);
    lowerRightX = lowerRightX < image.width() ? lowerRightX : image.width() -1;
    int lowerRightY = current[1] + ((image.height() / 10) / 2);
    lowerRightY = lowerRightY < image.height() ? lowerRightY : image.height() -1;
    Shape2 lowerRight(lowerRightX, lowerRightY);
    MultiArrayView<2, float> box = image.subarray(upperLeft, lowerRight);
    //3: get local max of view as next 
    int maxIndex = argMax(box);
    if (maxIndex == -1)
    {
            std::cout << "Something went wrong: argMax returned -1";
        return current;
    }
    Shape2 max(box.scanOrderIndexToCoordinate(maxIndex));
    Shape2 next(upperLeftX + max[0], upperLeftY + max[1]);
    if (next == current)
    {
        return next;
    }
    return localizeByFollowingLocalMaxima(image, next);
}

void mark(
    MultiArray& m,const Sequence& pts,
    const point& p,int dx,int dy)
{
    for(int k=p.x<0||p.y<0||p.x>=m.shape()[0]||p.y>=m.shape()[1]?1:0,
            k_end=distZ2inf(p,pts.begin(),pts.end()); k<k_end; ++k) {
        int x=p.x-k*dx;
        int y=p.y-k*dy;
        if(x>=0&&y>=0&&x<m.shape()[0]&&y<m.shape()[1])m[x][y]=true;
    }
}

void fft_X(MultiArray<complex<float>, 3> &a) {
    FFT<float> fft;
    int nx = a.dims()[0];
    for (int z = 0; z < a.dims()[2]; z++) {
        for (int y = 0; y < a.dims()[1]; y++) {
            VectorXcf fft_in = a.slice<1>({0,y,z},{nx,0,0}).asArray();
            VectorXcf fft_out(nx);
            fft.fwd(fft_out, fft_in);
            a.slice<1>({0,y,z},{nx,0,0}).asArray() = fft_out;
        }
    }
}

void fft_Y(MultiArray<complex<float>, 3> &a) {
    FFT<float> fft;
    int ny = a.dims()[1];
    for (int z = 0; z < a.dims()[2]; z++) {
        for (int x = 0; x < a.dims()[0]; x++) {
            VectorXcf fft_in = a.slice<1>({x,0,z},{0,ny,0}).asArray();
            VectorXcf fft_out(ny);
            fft.fwd(fft_out, fft_in);
            a.slice<1>({x,0,z},{0,ny,0}).asArray() = fft_out;
        }
    }
}

void fft_Z(MultiArray<complex<float>, 3> &a) {
    FFT<float> fft;
    int nz = a.dims()[2];
    for (int x = 0; x < a.dims()[0]; x++) {
        for (int y = 0; y < a.dims()[1]; y++) {
            VectorXcf fft_in = a.slice<1>({x,y,0},{0,0,nz}).asArray();
            VectorXcf fft_out(nz);
            fft.fwd(fft_out, fft_in);
            a.slice<1>({x,y,0},{0,0,nz}).asArray() = fft_out;
        }
    }
}

void ApplyFilter3D(MultiArray<complex<float>, 3> ks, MultiArray<float, 3> filter) {
    if ((ks.dims() != filter.dims()).any()) {
        throw(runtime_error("K-space and filter dimensions do not match."));
    }

    auto k_it = ks.begin();
    auto f_it = filter.begin();
    while (k_it != ks.end()) {
        *k_it = (*k_it) * (*f_it);
        k_it++;
        f_it++;
    }
}

Fields FieldAlgorithms::fieldsByErosionDilation(MultiArray<2, float> & image) 
{
    Shape2 shape = image.shape();
    MultiArray<2, float> source(image);
    MultiArray<2, float> t1(shape);
    MultiArray<2, float> t2(shape);

    double radius = 9; //Sinnvoll?

    multiGrayscaleErosion(source, t1, radius);
    multiGrayscaleDilation(source, t2, radius);

    MultiArray<2, float> psi_e = (t1- t2) * -1;
    //Opening of source
    multiGrayscaleErosion(source, t1, radius);
    multiGrayscaleDilation(t1, t2, radius);

    MultiArray<2, float> psi_p = source - t2;
    MultiArray<2, float> psi_pSmooth(psi_p.shape());
    gaussianSmoothMultiArray(psi_p, psi_pSmooth, 8.0);
    //Opening of source
    MultiArray<2, float> inverse = source * -1;
    multiGrayscaleErosion(inverse, t1, 9);
    multiGrayscaleDilation(t1, t2, 9);
    MultiArray<2, float> psi_v = t2 - inverse; 
    MultiArray<2, float> psi_vSmooth(psi_v.shape());
    gaussianSmoothMultiArray(psi_v, psi_vSmooth, 3.0);
    std::vector<Shape2> valleyLocals(0);
    Fields fields(psi_vSmooth, valleyLocals, psi_pSmooth, psi_e, image);
    return fields;
};

// MAIN
int main(int argc, char** argv) {

	if(argc != 3) {
		std::cout << "usage: " << argv[0] << " infile.sif outfile" << std::endl << std::endl;
		std::cout << "Converts one sif file to tiff stack" << std::endl;
		return 2;
	}

	std::string infile = argv[1];
	std::string outfile = argv[2];
	
    try
    {

		MultiArray<3,float> in;
		typedef MultiArray<3, float>::difference_type Shape;


		int width, height, stacksize;
		HDF5ImportInfo info(infile.c_str(), "/data");
		
		width = info.shapeOfDimension(0);
		height = info.shapeOfDimension(1);
		stacksize = info.shapeOfDimension(2);

		// create a 3D array of appropriate size
		in.reshape(Shape(width,height,stacksize));
		readHDF5(info, in); //Eingabe Bild

		std::cout << "Images with Shape: " << Shape(width, height, stacksize) << std::endl;
		std::cout << "Processing a stack of " << stacksize << " images..." << std::endl;


		exportVolume(in, VolumeExportInfo(outfile.c_str(), ".tif"));
		


    }
    catch (vigra::StdException & e)
    {
        std::cout<<"There was an error:"<<std::endl;
        std::cout << e.what() << std::endl;
        return 1;
    }
	
	
}

MultiArray<2, float > FieldAlgorithms::morphologyByGradientPattern(MultiArray<2, float> & image, MultiArray<2,float> &mask) {
    
    Shape2 shape(image.shape());
    MultiArray<2, TinyVector<float, 2>> imageGradients(shape);
    MultiArray<2, TinyVector<float, 2>> maskGradients(mask.shape());
    gaussianGradientMultiArray(image, imageGradients, 2.0);
    gaussianGradientMultiArray(mask, maskGradients, 2.0);

    //Threshold gradients with small magnitude, 
    //because we only consider directions from now on.
    MultiArray<2, float> magnitudes(shape);
    gaussianGradientMagnitude(image, magnitudes, 3.0);
    float thrhld = magnitudes[argMax(magnitudes)] * 0.3;
    thresholdGrad(magnitudes, imageGradients, thrhld);
    //The actual machting:
    return matchGradients(imageGradients, maskGradients);
   };