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

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 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 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 phase_correct_3(MultiArray<complex<float>, 3> & a, Agilent::FID &fid) {
    float ppe = fid.procpar().realValue("ppe");
    float ppe2 = fid.procpar().realValue("ppe2");

    float lpe = fid.procpar().realValue("lpe");
    float lpe2 = fid.procpar().realValue("lpe2");

    float ph = -2*M_PI*ppe/lpe;
    float ph2 = -2*M_PI*ppe2/lpe2;

    for (int z = 0; z < a.dims()[2]; z++) {
        const complex<float> fz = polar(1.f, ph2*z);
        for (int y = 0; y < a.dims()[1]; y++) {
            const complex<float> fy = polar(1.f, ph*y);
            for (int x = 0; x < a.dims()[0]; x++) {
                complex<float> val = a[{x,y,z}];
                val = val * fy * fz;
                a[{x,y,z}] = val;
            }
        }
    }
}

void fft_shift_3(MultiArray<complex<float>, 3> & a) {
    int x2 = a.dims()[0] / 2;
    int y2 = a.dims()[1] / 2;
    int z2 = a.dims()[2] / 2;

    MultiArray<complex<float>, 3>::Index size_half{x2,y2,z2};
    MultiArray<complex<float>, 3>::Index start_1{0,0,0};
    MultiArray<complex<float>, 3>::Index start_2{x2,y2,z2};

    for (int i = 0; i < 4; i++) {
        auto local_1 = start_1;
        auto local_2 = start_2;
        if (i < 3) {
            local_1[i] = a.dims()[i] / 2;
            local_2[i] = 0;
        }

        auto octant_1 = a.slice<3>({local_1}, {size_half});
        auto octant_2 = a.slice<3>({local_2}, {size_half});

        std::swap_ranges(octant_1.begin(), octant_1.end(), octant_2.begin());
    }
}

MultiArray<uint64,1> calcula_histograma(MultiArray<uint8,2>& in){
	//Histograma
	MultiArray<uint64,1> r(256);

	//Inicializamos todos los valores a cero
	r = 0;

	//Tamaño del MultiArray
	imageplus::uint64 x=0,y=0;
	x = in.dims(0);
	y = in.dims(1);

	for (uint64 j=0; j < x; ++j)
	for (uint64 i=0; i < y; ++i)
	{
		 r[in[j][i]] += 1;
	}
	return r;
}

MultiArray<uint8,2> ecualiza_histograma (const MultiArray<float64,2>& in){
	//Tamaño del MultiArray
	uint64 x=0,y=0;
	x = in.dims(0);
	y = in.dims(1);

	//Calculamos el histograma y lo normalizamos
	MultiArray<uint64,1> histogram;
	MultiArray<uint8,2> in_conv = convert<uint8>(in);
	histogram = calcula_histograma(in_conv);

	//std::cout << histogram << std::endl;

	//Calculamos la acumulada
	MultiArray<int64,1> prob(256);
	prob[0] = histogram[0];
	for(int64 i=1; i<256; i++){
		prob[i] = histogram[i] + prob[i-1];
	}

	//Imagen resultante
	MultiArray<uint8,2> result(x,y);
	int64 min = minval(prob);
	//std::cout << "Max " << maxval(in) << " Max converted " << int64(maxval(in_conv)) << " Min converted " << int64(minval(in_conv)) << std::endl;
	for (uint64 j=0; j < x; ++j)
	for (uint64 i=0; i < y; ++i)
	{
		 float64 num = (float64)(prob[in_conv[j][i]] - min);
		 float64 den = 1 - min;
		 float64 a = (num/den) * 255.0;

		 if( a > 255.0 ){
			 a = 255.0;
		 }else if( a < 0.0 ){
			 a = 0.0;
		 }

		 result[j][i] = mnint<uint8>(a);
	}
	return result;
}

//.........这里部分代码省略.........
        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;
}