void maxAbsolute(const hoNDArray<T>& x, T& r, size_t& ind)
    {
        size_t N = x.get_number_of_elements();
        const T* pX = x.begin();

        ind = 0;
        if ( N == 0 ) return;

        long long n;

        typename realType<T>::Type v = abs(pX[0]);
        typename realType<T>::Type v2;

        ind = 0;
        for ( n=1; n<(long long)N; n++ )
        {
            v2 = std::abs(pX[n]);
            if ( v2 > v )
            {
                v = v2;
                ind = n;
            }
        }

        r = pX[ind];
    }

    void correct_time_stamp_with_fitting(hoNDArray<float>& time_stamp, size_t startE1, size_t endE1)
    {
        try
        {
            size_t E1 = time_stamp.get_size(0);
            size_t N = time_stamp.get_size(1);
            size_t rE1 = endE1 - startE1 + 1;

            size_t e1, n;

            size_t num_acq_read_outs = 0;
            for ( n=0; n<N; n++ )
            {
                for ( e1=0; e1<E1; e1++ )
                {
                    if ( time_stamp(e1, n) > 0 )
                    {
                        num_acq_read_outs++;
                    }
                }
            }

            GDEBUG_STREAM(" Number of acquired lines : " << num_acq_read_outs);

            float a, b; // y = a + b*x
            {
                std::vector<float> x(num_acq_read_outs), y(num_acq_read_outs);

                size_t ind = 0;
                for ( n=0; n<N; n++ )
                {
                    for ( e1=startE1; e1<=endE1; e1++ )
                    {
                        float acq_time = time_stamp(e1, n);
                        if ( acq_time > 0 )
                        {
                            x[ind] = (float)(e1-startE1 + n*rE1);
                            y[ind] = acq_time;
                            ind++;
                        }
                    }
                }

                Gadgetron::simple_line_fit(x, y, a, b);
            }

            for ( n=0; n<N; n++ )
            {
                for ( e1=startE1; e1<=endE1; e1++ )
                {
                    float x_v = (float)(e1-startE1 + n*rE1);
                    time_stamp(e1, n) = a + b*x_v;
                }
            }
        }
        catch(...)
        {
            GADGET_THROW("Exceptions happened in correct_time_stamp_with_fitting(...) ... ");
        }
    }

void grappa2d_calib_convolution_kernel(const hoNDArray<T>& acsSrc, const hoNDArray<T>& acsDst, size_t accelFactor, double thres, size_t kRO, size_t kNE1, hoNDArray<T>& convKer)
{
    size_t startRO = 0;
    size_t endRO = acsSrc.get_size(0) - 1;
    size_t startE1 = 0;
    size_t endE1 = acsSrc.get_size(1) - 1;

    grappa2d_calib_convolution_kernel(acsSrc, acsDst, accelFactor, thres, kRO, kNE1, startRO, endRO, startE1, endE1, convKer);
}

int CS_FOCUSS_2D::fRecon(hoNDArray<std::complex<float> >  &hacfInput, hoNDArray<std::complex<float> >  &hacfRecon){

	// input dimensions
	vtDim_ = *hacfInput.get_dimensions();

	// number of channels
	iNChannels_ = (int)vtDim_[2];

	// if Matlab is used, initialize singleton variables
	/*if (bMatlab_){
		for (int iI = 0; iI < 20; iI++){
			GlobalVar_FOCUSS::instance()->vbStatPrinc_.push_back(false);
			GlobalVar_FOCUSS::instance()->vfPrincipleComponents_.push_back(new hoNDArray<std::complex<float> > ());
		}
	}*/

	// const complex values
	const std::complex<float> cfZero(0.0);
	const std::complex<float> cfOne(1.0);

	// store incoming data in array
	hoNDArray<std::complex<float> >  hacfKSpace = hacfInput;

	// permute kSpace: x-y-c -> y-x-c
	std::vector<size_t> vtDimOrder; vtDimOrder.push_back(1); vtDimOrder.push_back(0); vtDimOrder.push_back(2);
	hacfKSpace = *permute(&hacfKSpace, &vtDimOrder,false);

	// update dim_ vector
	vtDim_.clear(); vtDim_ = *hacfKSpace.get_dimensions();

	//------------------------------------------------------------------------
	//-------------------------- sampling mask -------------------------------
	//------------------------------------------------------------------------
	hoNDArray<std::complex<float> >  hacfFullMask(hacfKSpace.get_dimensions()); hacfFullMask.fill(cfZero);
	hoNDArray<bool> habFullMask(hacfKSpace.get_dimensions()); habFullMask.fill(false);
	pcfPtr_ = hacfKSpace.get_data_ptr();
	pcfPtr2_ = hacfFullMask.get_data_ptr();
	pbPtr_ = habFullMask.get_data_ptr();
	for (int i = 0; i < hacfKSpace.get_number_of_elements(); i++)
		if (pcfPtr_[i] != cfZero){
			pcfPtr2_[i] = cfOne;
			pbPtr_[i] = true;
		}

	//-------------------------------------------------------------------------
	//---------------- iFFT x direction - x ky kz ^= v (n�) -------------------
	//-------------------------------------------------------------------------
	if (Transform_fftBA_->get_active()){
		if (!bMatlab_ && bDebug_)
			#if __GADGETRON_VERSION_HIGHER_3_6__ == 1
				GDEBUG("FFT in read direction..\n");
			#else
				GADGET_DEBUG1("FFT in read direction..\n");
			#endif
		else if(bMatlab_ && bDebug_){
			// mexPrintf("FFT in read direction..\n");mexEvalString("drawnow;");
		}
		Transform_fftBA_->FTransform(hacfKSpace);
	}

    void compute_phase_time_stamp(const hoNDArray<float>& time_stamp, const hoNDArray<float>& cpt_time_stamp, size_t startE1, size_t endE1, 
        hoNDArray<float>& phs_time_stamp, hoNDArray<float>& phs_cpt_time_stamp)
    {
        try
        {
            size_t E1 = time_stamp.get_size(0);
            size_t N = time_stamp.get_size(1);
            size_t rE1 = endE1 - startE1 + 1;

            size_t e1, n;

            for ( n=0; n<N; n++ )
            {
                // phase time stamp as the mean of all aquired lines
                size_t num = 0;
                float tt = 0.0f;
                for ( e1=startE1; e1<=endE1; e1++ )
                {
                    if(time_stamp(e1, n)>0)
                    {
                        tt += time_stamp(e1, n);
                        num++;
                    }
                }
                phs_time_stamp(n, 0) = tt/((num>0) ? num : 1);

                //// phase cpt time as the median of all acquired lines
                //std::vector<float> cpt_buf(rE1, 0);
                //for ( e1=startE1; e1<=endE1; e1++ )
                //{
                //    if(cpt_time_stamp(e1, n)>=0)
                //        cpt_buf[e1-startE1] = cpt_time_stamp(e1, n);
                //}

                //std::sort(cpt_buf.begin(), cpt_buf.end());
                //phs_cpt_time_stamp(n, 0) = cpt_buf[E1/2-startE1];

                // phase cpt time as the cpt time of center line
                phs_cpt_time_stamp(n, 0) = cpt_time_stamp(E1/2, n);
            }
        }
        catch(...)
        {
            GADGET_THROW("Exceptions happened in compute_phase_time_stamp(...) ... ");
        }
    }

void apply_unmix_coeff_kspace(const hoNDArray<T>& kspace, const hoNDArray<T>& unmixCoeff, hoNDArray<T>& complexIm)
{
    try
    {
        GADGET_CHECK_THROW(kspace.get_size(0) == unmixCoeff.get_size(0));
        GADGET_CHECK_THROW(kspace.get_size(1) == unmixCoeff.get_size(1));
        GADGET_CHECK_THROW(kspace.get_size(2) == unmixCoeff.get_size(2));

        hoNDArray<T> buffer2DT(kspace);
        GADGET_CATCH_THROW(Gadgetron::hoNDFFT<typename realType<T>::Type>::instance()->ifft2c(kspace, buffer2DT));

        std::vector<size_t> dim;
        kspace.get_dimensions(dim);
        dim[2] = 1;

        if (!complexIm.dimensions_equal(&dim))
        {
            complexIm.create(&dim);
        }

        Gadgetron::multiply(buffer2DT, unmixCoeff, buffer2DT);
        Gadgetron::sum_over_dimension(buffer2DT, complexIm, 2);
    }
    catch (...)
    {
        GADGET_THROW("Errors in apply_unmix_coeff_kspace(const hoNDArray<T>& kspace, const hoNDArray<T>& unmixCoeff, hoNDArray<T>& complexIm) ... ");
    }
}

void apply_unmix_coeff_aliased_image(const hoNDArray<T>& aliasedIm, const hoNDArray<T>& unmixCoeff, hoNDArray<T>& complexIm)
{
    try
    {
        GADGET_CHECK_THROW(aliasedIm.get_size(0) == unmixCoeff.get_size(0));
        GADGET_CHECK_THROW(aliasedIm.get_size(1) == unmixCoeff.get_size(1));
        GADGET_CHECK_THROW(aliasedIm.get_size(2) == unmixCoeff.get_size(2));

        std::vector<size_t> dim;
        aliasedIm.get_dimensions(dim);
        dim[2] = 1;

        if (!complexIm.dimensions_equal(&dim))
        {
            complexIm.create(&dim);
        }

        hoNDArray<T> buffer2DT(aliasedIm);

        Gadgetron::multiply(aliasedIm, unmixCoeff, buffer2DT);
        Gadgetron::sum_over_dimension(buffer2DT, complexIm, 2);
    }
    catch (...)
    {
        GADGET_THROW("Errors in apply_unmix_coeff_aliased_image(const hoNDArray<T>& aliasedIm, const hoNDArray<T>& unmixCoeff, hoNDArray<T>& complexIm) ... ");
    }
}

    hoNDArray<uint16_t>
    update_field_map(const hoNDArray<uint16_t> &field_map_index, const hoNDArray<uint16_t> &proposed_field_map_index,
                     const hoNDArray<float> &residuals_map, const hoNDArray<float> &lambda_map) {


        hoNDArray<float> residual_diff_map(field_map_index.get_dimensions());
        const auto X = field_map_index.get_size(0);
        const auto Y = field_map_index.get_size(1);
        const auto Z = field_map_index.get_size(2);

        for (size_t kz = 0; kz < Z; kz++) {
            for (size_t ky = 0; ky < Y; ky++) {
                for (size_t kx = 0; kx < X; kx++) {
                    residual_diff_map(kx, ky,kz) = residuals_map(field_map_index(kx, ky,kz), kx, ky,kz) -
                                                residuals_map(proposed_field_map_index(kx, ky,kz), kx, ky,kz);


                }
            }
        }


        std::vector<boost::default_color_type> color_map;
        if (Z == 1) {
            color_map = graph_cut<2>(field_map_index, proposed_field_map_index, lambda_map,
                                     residual_diff_map);
        } else {
            color_map = graph_cut<3>(field_map_index, proposed_field_map_index, lambda_map, residual_diff_map);
        }



        auto result = field_map_index;
        size_t updated_voxels = 0;
        for (size_t i = 0; i < field_map_index.get_number_of_elements(); i++) {
            if (color_map[i] != boost::default_color_type::black_color) {
                updated_voxels++;
                result[i] = proposed_field_map_index[i];
            }
        }

        return result;

    }

    void maxValue(const hoNDArray<T>& a, T& v)
    {
        typedef T ValueType;

        try
        {
            const ValueType* pA = a.begin();
            size_t n = a.get_number_of_elements();
            v = pA[0];

            size_t ii;
            for (ii=1; ii<n; ii++)
            {
                if (pA[ii]>v) v = pA[ii];
            }
        }
        catch(...)
        {
            GADGET_THROW("Errors in maxValue(const hoNDArray<T>& a, T& v) ... ");
        }
    }

    void sort(const hoNDArray<T>& x, hoNDArray<T>& r, bool isascending)
    {
        if ( &r != &x )
        {
            if ( r.get_number_of_elements()!=x.get_number_of_elements())
            {
                r = x;
            }
            else
            {
                memcpy(r.begin(), x.begin(), x.get_number_of_bytes());
            }
        }

        sort(x.get_number_of_elements(), x.begin(), r.begin(), isascending);
    }

void outputPlotIm(const hoNDArray<unsigned char>& im, bool trueColor, hoNDArray<float>& plotIm)
{
    size_t xsize = im.get_size(1);
    size_t ysize = im.get_size(2);

    plotIm.copyFrom(im);

    if (trueColor)
    {
        std::vector<size_t> dim_order(3);
        dim_order[0] = 1;
        dim_order[1] = 2;
        dim_order[2] = 0;

        hoNDArray<float> plotImPermuted;
        plotImPermuted.copyFrom(plotIm);

        plotIm.create(xsize, ysize, 3);

        Gadgetron::permute(plotImPermuted, plotIm, dim_order);
    }
    else
    {
        hoNDArray<float> plotIm2D;
        Gadgetron::sum_over_dimension(plotIm, plotIm2D, 0);

        plotIm2D.squeeze();

        std::vector<size_t> dim_order(2);
        dim_order[0] = 0;
        dim_order[1] = 1;

        plotIm.create(xsize, ysize);
        Gadgetron::permute(plotIm2D, plotIm, dim_order);
    }
}

    hoNDArray< std::complex<float> > fatwater_separation(hoNDArray< std::complex<float> >& data, FatWaterParameters p, FatWaterAlgorithm a)
    {

	//Get some data parameters
	//7D, fixed order [X, Y, Z, CHA, N, S, LOC]
        uint16_t X = data.get_size(0);
        uint16_t Y = data.get_size(1);
        uint16_t Z = data.get_size(2);
        uint16_t CHA = data.get_size(3);
        uint16_t N = data.get_size(4);
        uint16_t S = data.get_size(5);
        uint16_t LOC = data.get_size(6);

	GDEBUG("Size of my array: %d, %d, %d .\n", X,Y,Z);

	hoNDArray< std::complex<float> > out(X,Y,Z,CHA,N,2,LOC); // S dimension gets replaced by water/fat stuff

	float fieldStrength = p.fieldStrengthT_;
        std::vector<float> echoTimes = p.echoTimes_;
	bool precessionIsClockwise = p.precessionIsClockwise_;
        for (auto& te: echoTimes) {
          te = te*0.001; // Echo times in seconds rather than milliseconds
        }

	GDEBUG("In toolbox - Field Strength: %f T \n", fieldStrength);
        for (auto& te: echoTimes) {
	  GDEBUG("In toolbox - Echo time: %f seconds \n", te);
        }
	GDEBUG("In toolbox - PrecessionIsClockwise: %d \n", precessionIsClockwise);

	//Get or set some algorithm parameters
	//Gadgetron::ChemicalSpecies w = a.species_[0];
	//Gadgetron::ChemicalSpecies f = a.species_[1];

	//	GDEBUG("In toolbox - Fat peaks: %f  \n", f.ampFreq_[0].first);
	//	GDEBUG("In toolbox - Fat peaks 2: %f  \n", f.ampFreq_[0].second);

	// Set some initial parameters so we can get going
	// These will have to be specified in the XML file eventually
	std::pair<float,float> range_r2star = std::make_pair(0.0,0.0);
	uint16_t num_r2star = 1;
	std::pair<float,float> range_fm = std::make_pair(-80.0,80.0);
	uint16_t num_fm = 101;
	uint16_t num_iterations = 40;
	uint16_t subsample = 1;
	float lmap_power = 2.0;
	float lambda = 0.02;
	float lambda_extra = 0.02;

	//Check that we have reasonable data for fat-water separation


	//Calculate residual
	//
	float relAmp, freq_hz;
	uint16_t npeaks;
	uint16_t nspecies = a.species_.size();
	uint16_t nte = echoTimes.size();
	GDEBUG("In toolbox - NTE: %d \n", nte);

	hoMatrix< std::complex<float> > phiMatrix(nte,nspecies);
	for( int k1=0;k1<nte;k1++) {
	  for( int k2=0;k2<nspecies;k2++) {
	    phiMatrix(k1,k2) = 0.0;
	    npeaks = a.species_[k2].ampFreq_.size();
	    for( int k3=0;k3<npeaks;k3++) {
	      relAmp = a.species_[k2].ampFreq_[k3].first;
	      freq_hz = fieldStrength*GAMMABAR*a.species_[k2].ampFreq_[k3].second;
	      phiMatrix(k1,k2) += relAmp*std::complex<float>(cos(2*PI*echoTimes[k1]*freq_hz),sin(2*PI*echoTimes[k1]*freq_hz));
	    }
	    GDEBUG("Cur value phiMatrix = (%f,%f) \n", phiMatrix(k1,k2).real(), phiMatrix(k1,k2).imag());
	  }
	}
	//auto a_phiMatrix = as_arma_matrix(&phiMatrix);
	//auto mymat2 = mymat.t()*mymat;

	for(int ka=0;ka<phiMatrix.get_size(0);ka++) {
	  for(int kb=0;kb<phiMatrix.get_size(1);kb++) {
	    GDEBUG("Check phiMatrix(%d,%d) = %f + i %f \n", ka,kb,phiMatrix(ka,kb).real(),phiMatrix(ka,kb).imag());
	  }
	}




	hoMatrix< std::complex<float> > IdentMat(nte,nte);
	for( int k1=0;k1<nte;k1++) {
	  for( int k2=0;k2<nte;k2++) {
	    if( k1==k2 ) {
	      IdentMat(k1,k2) = std::complex<float>(1.0,0.0);
	    } else {
	      IdentMat(k1,k2) = std::complex<float>(0.0,0.0);
	    }
	  }
	}
	//	auto a_phiMatrix = as_arma_matrix(&IdentMat);

	float fm;
	std::vector<float> fms(num_fm);
	fms[0] = range_fm.first;
//.........这里部分代码省略.........

    void correct_heart_beat_time_stamp_with_fitting(hoNDArray<float>& cpt_time_stamp, hoNDArray<int>& ind_hb, size_t startE1, size_t endE1, 
                                                const std::vector<size_t>& start_e1_hb, const std::vector<size_t>& end_e1_hb, 
                                                const std::vector<size_t>& start_n_hb, const std::vector<size_t>& end_n_hb )
    {
        try
        {
            size_t E1 = cpt_time_stamp.get_size(0);
            size_t N = cpt_time_stamp.get_size(1);
            size_t rE1 = endE1-startE1+1;

            size_t e1, n, ind, ii;

            size_t num_acq_read_outs = 0;
            for ( n=0; n<N; n++ )
            {
                for ( e1=0; e1<E1; e1++ )
                {
                    if ( cpt_time_stamp(e1, n) >= 0 )
                    {
                        num_acq_read_outs++;
                    }
                }
            }

            size_t numOfHB = start_e1_hb.size();

            std::vector<size_t> ind_HB_start(numOfHB);
            std::vector<size_t> ind_HB_end(numOfHB);

            for ( ind=0; ind<numOfHB; ind++ )
            {
                ind_HB_start[ind] = start_e1_hb[ind] + start_n_hb[ind] * E1;
                ind_HB_end[ind] = end_e1_hb[ind] + end_n_hb[ind] * E1;
            }

            // --------------------------------------------------------
            // fit a line to every heart beat
            // --------------------------------------------------------
            float a, b;
            std::vector<float> A(numOfHB, 0.0f), B(numOfHB, 0.0f);
            for ( ind=0; ind<numOfHB; ind++ )
            {
                std::vector<float> x, y;

                size_t cpt;
                for ( cpt=ind_HB_start[ind]; cpt<=ind_HB_end[ind]; cpt++ )
                {
                    size_t n = cpt / E1;
                    size_t e1 = cpt - n*E1;

                    if(e1>=startE1 && e1<=endE1)
                    {
                        if ( cpt_time_stamp[cpt] > -1 )
                        {
                            size_t x_ind = (e1-startE1) + n*rE1;
                            x.push_back( (float)x_ind );
                            y.push_back(cpt_time_stamp[cpt]);
                        }
                    }
                }

                if ( !x.empty() )
                {
                    Gadgetron::simple_line_fit(x, y, a, b);
                    A[ind] = a;
                    B[ind] = b;
                }
            }

            // --------------------------------------------------------
            // compute cpt time stamp for every line
            // --------------------------------------------------------
            size_t num = cpt_time_stamp.get_number_of_elements();
            for ( ind=0; ind<num; ind++ )
            {
                n = ind / E1;
                e1 = ind - n*E1;

                if(e1>=startE1 && e1<=endE1)
                {
                    // find to which heart beat this line belongs
                    bool foundHB = false;
                    for ( ii=0; ii<numOfHB; ii++ )
                    {
                        size_t startHB = ind_HB_start[ii];
                        size_t endHB = ind_HB_end[ii];

                        if ( ii==0 && ind<=startHB )
                        {
                            foundHB = true;
                            break;
                        }

                        if ( ii==numOfHB-1 && ind>=endHB )
                        {
                            foundHB = true;
                            break;
                        }

                        if ( ind>=startHB && ind<=endHB )
//.........这里部分代码省略.........

 void dotu(const hoNDArray<T>& x, const hoNDArray<T>& y, T& r)
 {
     GADGET_DEBUG_CHECK_THROW(x.get_number_of_elements()==y.get_number_of_elements());
     dotu(x.get_number_of_elements(), x.begin(), y.begin(), r);
 }

 template<class T> size_t amax(const hoNDArray<T>& x)
 {
     return amax(x.get_number_of_elements(), x.begin());
 }

void grappa2d_convert_to_convolution_kernel(const hoNDArray<T>& ker, size_t kRO, const std::vector<int>& kE1, const std::vector<int>& oE1, hoNDArray<T>& convKer)
{
    try
    {
        long long srcCHA = (long long)(ker.get_size(2));
        long long dstCHA = (long long)(ker.get_size(3));
        long long kNE1 = (long long)(kE1.size());
        long long oNE1 = (long long)(oE1.size());

        long long kROhalf = kRO / 2;
        if (2 * kROhalf == kRO)
        {
            GWARN_STREAM("grappa2d_convert_to_convolution_kernel - 2*kROhalf == kRO " << kRO);
        }
        kRO = 2 * kROhalf + 1;

        //// fill the convolution kernels
        long long convKRO = 2 * kRO + 3;

        long long maxKE1 = std::abs(kE1[0]);
        if (std::abs(kE1[kNE1 - 1]) > maxKE1)
        {
            maxKE1 = std::abs(kE1[kNE1 - 1]);
        }
        long long convKE1 = 2 * maxKE1 + 1;

        //// allocate the convolution kernel
        convKer.create(convKRO, convKE1, srcCHA, dstCHA);
        Gadgetron::clear(&convKer);

        //// index
        long long oe1, kro, ke1, src, dst;

        //// fill the convolution kernel and sum up multiple kernels
        for (oe1 = 0; oe1<oNE1; oe1++)
        {
            for (ke1 = 0; ke1<kNE1; ke1++)
            {
                for (kro = -kROhalf; kro <= kROhalf; kro++)
                {
                    for (dst = 0; dst<dstCHA; dst++)
                    {
                        for (src = 0; src<srcCHA; src++)
                        {
                            convKer(-kro + kRO + 1, oE1[oe1] - kE1[ke1] + maxKE1, src, dst) = ker(kro + kROhalf, ke1, src, dst, oe1);
                        }
                    }

                }
            }
        }

        if (oE1[0] != 0)
        {
            for (dst = 0; dst<dstCHA; dst++)
            {
                convKer(kRO + 1, maxKE1, dst, dst) = 1.0;
            }
        }
    }
    catch (...)
    {
        GADGET_THROW("Errors in grappa2d_convert_to_convolution_kernel(...) ... ");
    }

    return;
}

    void detect_heart_beat_with_time_stamp(hoNDArray<float>& cpt_time_stamp, hoNDArray<int>& ind_hb, 
                                        std::vector<size_t>& start_e1_hb, std::vector<size_t>& end_e1_hb, 
                                        std::vector<size_t>& start_n_hb, std::vector<size_t>& end_n_hb )
    {
        try
        {
            size_t E1 = cpt_time_stamp.get_size(0);
            size_t N = cpt_time_stamp.get_size(1);

            size_t e1, n, ind, ii;

            size_t num_acq_read_outs = 0;
            for ( n=0; n<N; n++ )
            {
                for ( e1=0; e1<E1; e1++ )
                {
                    if ( cpt_time_stamp(e1, n) >= 0 )
                    {
                        num_acq_read_outs++;
                    }
                }
            }

            ind_hb.create(E1, N);
            Gadgetron::clear(ind_hb);

            // --------------------------------------------------------
            // cpt time stamps
            // --------------------------------------------------------

            std::vector<float> acquired_cpt(num_acq_read_outs);
            std::vector<size_t> ind_acquired_cpt(num_acq_read_outs);

            ind = 0;
            for ( n=0; n<N; n++ )
            {
                for ( e1=0; e1<E1; e1++ )
                {
                    if ( cpt_time_stamp(e1, n) > -1 )
                    {
                        acquired_cpt[ind] = cpt_time_stamp(e1, n);
                        ind_acquired_cpt[ind] = e1 + n*E1;
                        ind++;
                    }
                }
            }

            // --------------------------------------------------------
            // find the number of heart beats
            // --------------------------------------------------------
            size_t numOfHB = 0;

            // store the line indexes for every heart beat
            std::vector<size_t> ind_HB_start, ind_HB_end;
            ind_HB_start.push_back(0);

            for ( ind=1; ind<num_acq_read_outs; ind++ )
            {
                if ( acquired_cpt[ind] < acquired_cpt[ind-1] )
                {
                    // find a new heart beat
                    numOfHB++;

                    size_t end_ind_prev_HB = ind_acquired_cpt[ind-1];
                    size_t start_ind_curr_HB = ind_acquired_cpt[ind];

                    // if there is a gap between end and start ind, fill the gap
                    if ( end_ind_prev_HB+1 != start_ind_curr_HB )
                    {
                        long long gap = start_ind_curr_HB - end_ind_prev_HB - 1;
                        if ( gap % 2 == 0 )
                        {
                            end_ind_prev_HB += gap;
                        }
                        else
                        {
                            end_ind_prev_HB += gap;
                        }

                        if ( end_ind_prev_HB+1 != start_ind_curr_HB )
                        {
                            GWARN_STREAM("end_ind_prev_HB+1 ~= start_ind_curr_HB : " << end_ind_prev_HB << " " << start_ind_curr_HB);
                        }
                    }

                    ind_HB_end.push_back( end_ind_prev_HB );
                    ind_HB_start.push_back( start_ind_curr_HB );
                }
            }

            ind_HB_end.push_back( E1*N-1 );
            numOfHB = ind_HB_end.size();

            // --------------------------------------------------------
            // fill the start and end indexes
            // --------------------------------------------------------
            start_e1_hb.resize(numOfHB, 0);
            end_e1_hb.resize(numOfHB, 0);

            start_n_hb.resize(numOfHB, 0);
//.........这里部分代码省略.........

    void GenericReconCartesianGrappaGadget::prepare_down_stream_coil_compression_ref_data(
            const hoNDArray<std::complex<float> > &ref_src, hoNDArray<std::complex<float> > &ref_coil_map,
            hoNDArray<std::complex<float> > &ref_dst, size_t e) {

        if (!downstream_coil_compression.value()) {
            GDEBUG_CONDITION_STREAM(verbose.value(), "Downstream coil compression is not prescribed ... ");
            ref_dst = ref_src;
            return;
        }

        if (downstream_coil_compression_thres.value() < 0 && downstream_coil_compression_num_modesKept.value() == 0) {
            GDEBUG_CONDITION_STREAM(verbose.value(),
                                    "Downstream coil compression is prescribed to use all input channels ... ");
            ref_dst = ref_src;
            return;
        }

        // determine how many channels to use
        size_t RO = ref_src.get_size(0);
        size_t E1 = ref_src.get_size(1);
        size_t E2 = ref_src.get_size(2);
        size_t CHA = ref_src.get_size(3);
        size_t N = ref_src.get_size(4);
        size_t S = ref_src.get_size(5);
        size_t SLC = ref_src.get_size(6);

        size_t recon_RO = ref_coil_map.get_size(0);
        size_t recon_E1 = ref_coil_map.get_size(1);
        size_t recon_E2 = ref_coil_map.get_size(2);

        std::complex<float> *pRef = const_cast< std::complex<float> * >(ref_src.begin());

        size_t dstCHA = CHA;
        if (downstream_coil_compression_num_modesKept.value() > 0 &&
            downstream_coil_compression_num_modesKept.value() <= CHA) {
            dstCHA = downstream_coil_compression_num_modesKept.value();
        } else {
            std::vector<float> E(CHA, 0);
            long long cha;

#pragma omp parallel default(none) private(cha) shared(RO, E1, E2, CHA, pRef, E)
            {
                hoNDArray<std::complex<float> > dataCha;
#pragma omp for
                for (cha = 0; cha < (long long) CHA; cha++) {
                    dataCha.create(RO, E1, E2, pRef + cha * RO * E1 * E2);
                    float v = Gadgetron::nrm2(dataCha);
                    E[cha] = v * v;
                }
            }

            for (cha = 1; cha < (long long) CHA; cha++) {
                if (std::abs(E[cha]) < downstream_coil_compression_thres.value() * std::abs(E[0])) {
                    break;
                }
            }

            dstCHA = cha;
        }

        GDEBUG_CONDITION_STREAM(verbose.value(),
                                "Downstream coil compression is prescribed to use " << dstCHA << " out of " << CHA
                                                                                    << " channels ...");

        if (dstCHA < CHA) {
            ref_dst.create(RO, E1, E2, dstCHA, N, S, SLC);
            hoNDArray<std::complex<float> > ref_coil_map_dst;
            ref_coil_map_dst.create(recon_RO, recon_E1, recon_E2, dstCHA, N, S, SLC);

            size_t slc, s, n;
            for (slc = 0; slc < SLC; slc++) {
                for (s = 0; s < S; s++) {
                    for (n = 0; n < N; n++) {
                        std::complex<float> *pDst = &(ref_dst(0, 0, 0, 0, n, s, slc));
                        const std::complex<float> *pSrc = &(ref_src(0, 0, 0, 0, n, s, slc));
                        memcpy(pDst, pSrc, sizeof(std::complex<float>) * RO * E1 * E2 * dstCHA);

                        pDst = &(ref_coil_map_dst(0, 0, 0, 0, n, s, slc));
                        pSrc = &(ref_coil_map(0, 0, 0, 0, n, s, slc));
                        memcpy(pDst, pSrc, sizeof(std::complex<float>) * recon_RO * recon_E1 * recon_E2 * dstCHA);
                    }
                }
            }

            ref_coil_map = ref_coil_map_dst;
        } else {
            ref_dst = ref_src;
        }

    }

    void GenericReconCartesianNonLinearSpirit2DTGadget::perform_nonlinear_spirit_unwrapping(hoNDArray< std::complex<float> >& kspace, 
        hoNDArray< std::complex<float> >& kerIm, hoNDArray< std::complex<float> >& ref2DT, hoNDArray< std::complex<float> >& coilMap2DT, hoNDArray< std::complex<float> >& res, size_t e)
    {
        try
        {
            bool print_iter = this->spirit_print_iter.value();

            size_t RO = kspace.get_size(0);
            size_t E1 = kspace.get_size(1);
            size_t E2 = kspace.get_size(2);
            size_t CHA = kspace.get_size(3);
            size_t N = kspace.get_size(4);
            size_t S = kspace.get_size(5);
            size_t SLC = kspace.get_size(6);

            size_t ref_N = kerIm.get_size(4);
            size_t ref_S = kerIm.get_size(5);

            hoNDArray< std::complex<float> > kspaceLinear(kspace);
            res = kspace;

            // detect whether random sampling is used
            bool use_random_sampling = false;
            std::vector<long long> sampled_step_size;
            long long n, e1;
            for (n=0; n<(long long)N; n++)
            {
                long long prev_sampled_line = -1;
                for (e1=0; e1<(long long)E1; e1++)
                {
                    if(std::abs(kspace(RO/2, e1, 0, 0, 0, 0, 0))>0 && std::abs(kspace(RO/2, e1, 0, CHA-1, 0, 0, 0))>0)
                    {
                        if(prev_sampled_line>0)
                        {
                            sampled_step_size.push_back(e1 - prev_sampled_line);
                        }

                        prev_sampled_line = e1;
                    }
                }
            }

            if(sampled_step_size.size()>4)
            {
                size_t s;
                for (s=2; s<sampled_step_size.size()-1; s++)
                {
                    if(sampled_step_size[s]!=sampled_step_size[s-1])
                    {
                        use_random_sampling = true;
                        break;
                    }
                }
            }

            if(use_random_sampling)
            {
                GDEBUG_STREAM("SPIRIT Non linear, random sampling is detected ... ");
            }

            Gadgetron::GadgetronTimer timer(false);

            // compute linear solution as the initialization
            if(use_random_sampling)
            {
                if (this->perform_timing.value()) timer.start("SPIRIT Non linear, perform linear spirit recon ... ");
                this->perform_spirit_unwrapping(kspace, kerIm, kspaceLinear);
                if (this->perform_timing.value()) timer.stop();
            }
            else
            {
                if (this->perform_timing.value()) timer.start("SPIRIT Non linear, perform linear recon ... ");

                size_t ref2DT_RO = ref2DT.get_size(0);
                size_t ref2DT_E1 = ref2DT.get_size(1);

                // mean over N
                hoNDArray< std::complex<float> > meanKSpace;
                Gadgetron::sum_over_dimension(ref2DT, meanKSpace, 4);

                // if (!debug_folder_full_path_.empty()) { gt_exporter_.export_array_complex(meanKSpace, debug_folder_full_path_ + "spirit_nl_2DT_meanKSpace"); }

                hoNDArray< std::complex<float> > acsSrc(ref2DT_RO, ref2DT_E1, CHA, meanKSpace.begin());
                hoNDArray< std::complex<float> > acsDst(ref2DT_RO, ref2DT_E1, CHA, meanKSpace.begin());

                double grappa_reg_lamda = 0.0005;
                size_t kRO = 5;
                size_t kE1 = 4;

                hoNDArray< std::complex<float> > convKer;
                hoNDArray< std::complex<float> > kIm(RO, E1, CHA, CHA);

                Gadgetron::grappa2d_calib_convolution_kernel(acsSrc, acsDst, (size_t)this->acceFactorE1_[e], grappa_reg_lamda, kRO, kE1, convKer);
                Gadgetron::grappa2d_image_domain_kernel(convKer, RO, E1, kIm);

                // if (!debug_folder_full_path_.empty()) gt_exporter_.export_array_complex(kIm, debug_folder_full_path_ + "spirit_nl_2DT_kIm");

                Gadgetron::hoNDFFT<float>::instance()->ifft2c(kspace, complex_im_recon_buf_);
                // if (!debug_folder_full_path_.empty()) gt_exporter_.export_array_complex(complex_im_recon_buf_, debug_folder_full_path_ + "spirit_nl_2DT_aliasedImage");

//.........这里部分代码省略.........

void grappa2d_unmixing_coeff(const hoNDArray<T>& kerIm, const hoNDArray<T>& coilMap, size_t acceFactorE1, hoNDArray<T>& unmixCoeff, hoNDArray< typename realType<T>::Type >& gFactor)
{
    try
    {
        typedef typename realType<T>::Type value_type;

        size_t RO = kerIm.get_size(0);
        size_t E1 = kerIm.get_size(1);
        size_t srcCHA = kerIm.get_size(2);
        size_t dstCHA = kerIm.get_size(3);

        GADGET_CHECK_THROW(acceFactorE1 >= 1);

        GADGET_CHECK_THROW(coilMap.get_size(0) == RO);
        GADGET_CHECK_THROW(coilMap.get_size(1) == E1);
        GADGET_CHECK_THROW(coilMap.get_size(2) == dstCHA);

        std::vector<size_t> dimUnmixing(3);
        dimUnmixing[0] = RO; dimUnmixing[1] = E1; dimUnmixing[2] = srcCHA;