C++ MultiArray::begin方法代码示例

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++;
    }
}

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;
}

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());
}

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
}

//.........这里部分代码省略.........
        string inPath(argv[optind++]);
        size_t fileSep = inPath.find_last_of("/") + 1;
        size_t fileExt = inPath.find_last_of(".");
        if ((fileExt == string::npos) || (inPath.substr(fileExt) != ".fid")) {
            cerr << inPath << " is not a valid .fid directory" << endl;
        }
        string outPath = outPrefix + inPath.substr(fileSep, fileExt - fileSep) + ".nii";
        if (zip)
            outPath = outPath + ".gz";

        Agilent::FID fid(inPath);


        string apptype = fid.procpar().stringValue("apptype");
        string seqfil  = fid.procpar().stringValue("seqfil");

        if (apptype != "im3D") {
            cerr << "apptype " << apptype << " not supported, skipping." << endl;
            continue;
        }

        if (verbose) {
            cout << fid.print_info() << endl;
            cout << "apptype = " << apptype << endl;
            cout << "seqfil  = " << seqfil << endl;
        }

        /*
         * Assemble k-Space
         */
        MultiArray<complex<float>, 4> vols;
        if (seqfil.substr(0, 5) == "mge3d") {
            vols = reconMGE(fid);
        } else if (seqfil.substr(0, 7) == "mp3rage") {
            vols = reconMP2RAGE(fid);
        } else {
            cerr << "Recon for " << seqfil << " not implemented, skipping." << endl;
            continue;
        }

        /*
         * Build and apply filter
         */
        MultiArray<float, 3> filter;
        switch (filterType) {
        case Filters::None: break;
        case Filters::Hanning:
            if (verbose) cout << "Building Hanning filter" << endl;
            filter = Hanning3D(vols.dims().head(3), f_a);
            break;
        case Filters::Tukey:
            if (verbose) cout << "Building Tukey filter" << endl;
            filter = Tukey3D(vols.dims().head(3), f_a, f_q);
            break;
        }
        if (filterType != Filters::None) {
            if (verbose) cout << "Applying filter" << endl;
            for (int v = 0; v < vols.dims()[3]; v++) {
                MultiArray<complex<float>, 3> vol = vols.slice<3>({0,0,0,v},{-1,-1,-1,0});
                ApplyFilter3D(vol,filter);
            }
        }
        /*
         * FFT
         */
        if (!kspace) {
            for (int v = 0; v < vols.dims()[3]; v++) {
                if (verbose) cout << "FFTing vol " << v << endl;
                MultiArray<complex<float>, 3> vol = vols.slice<3>({0,0,0,v},{-1,-1,-1,0});
                phase_correct_3(vol, fid);
                fft_shift_3(vol);
                fft_X(vol);
                fft_Y(vol);
                fft_Z(vol);
                fft_shift_3(vol);
            }
        }

        if (verbose) cout << "Writing file: " << outPath << endl;
        list<Nifti::Extension> exts;
        if (procpar) {
            if (verbose) cout << "Embedding procpar" << endl;
            ifstream pp_file(inPath + "/procpar", ios::binary);
            pp_file.seekg(ios::end);
            size_t fileSize = pp_file.tellg();
            pp_file.seekg(ios::beg);
            vector<char> data; data.reserve(fileSize);
            data.assign(istreambuf_iterator<char>(pp_file), istreambuf_iterator<char>());
            exts.emplace_back(NIFTI_ECODE_COMMENT, data);
        }
        Affine3f xform  = scale * fid.procpar().calcTransform();
        ArrayXf voxdims = (Affine3f(xform.rotation()).inverse() * xform).matrix().diagonal();
        Nifti::Header outHdr(vols.dims(), voxdims, dtype);
        outHdr.setTransform(xform);
        Nifti::File output(outHdr, outPath, exts);
        output.writeVolumes(vols.begin(), vols.end(), 0, vols.dims()[3]);
        output.close();
    }
    return 0;
}