C++ MetaData::getValue方法代码示例

/* Read coordinates from disk ---------------------------------------------- */
void Micrograph::read_coordinates(int label, const FileName &_fn_coords)
{
    std::ifstream fh;
    int line_no = 0;
    std::string line;

    fn_coords = _fn_coords;

    MetaData MD;
    MD.read(fn_coords);
    line_no = MD.size();

    // Resize coordinate list and read
    coords.reserve(line_no);
    struct Particle_coords aux;
    aux.valid = true;
    aux.label = label;
    aux.cost = 1;

    FOR_ALL_OBJECTS_IN_METADATA(MD)
    {
        MD.getValue(MDL_XCOOR, aux.X, __iter.objId); //aux.X=x;
        MD.getValue(MDL_YCOOR, aux.Y, __iter.objId); //aux.Y=y;
        coords.push_back(aux);
    }
}

void ProgValidationNonTilt::obtainSumW(const MetaData & tempMd, double & sum_W, std::vector<double> & sum_u, std::vector<double> & H, const double factor)
{
    double a;
    double rot,tilt,w;
    double x,y,z;
    double xx,yy,zz;
    double w2;
    double tempW;
    double W;
    double sumW;

    sumW = 0;
    FOR_ALL_OBJECTS_IN_METADATA(tempMd)
    {
        tempMd.getValue(MDL_ANGLE_ROT,rot,__iter.objId);
        tempMd.getValue(MDL_ANGLE_TILT,tilt,__iter.objId);
        //tempMd.getValue(MDL_WEIGHT,w,__iter.objId);
        tempMd.getValue(MDL_MAXCC,w,__iter.objId);
        x = sin(tilt*PI/180.)*cos(rot*PI/180.);
        y = sin(tilt*PI/180.)*sin(rot*PI/180.);
        z = std::abs(cos(tilt*PI/180.));

        tempW = 1e3;
        _FOR_ALL_OBJECTS_IN_METADATA2(tempMd)
        {
            tempMd.getValue(MDL_ANGLE_ROT,rot,__iter2.objId);
            tempMd.getValue(MDL_ANGLE_TILT,tilt,__iter2.objId);
            //tempMd.getValue(MDL_WEIGHT,w2,__iter2.objId);
            tempMd.getValue(MDL_MAXCC,w2,__iter2.objId);
            xx = sin(tilt*PI/180.)*cos(rot*PI/180.);
            yy = sin(tilt*PI/180.)*sin(rot*PI/180.);
            zz = std::abs(cos(tilt*PI/180.));
            a = std::abs(std::acos(x*xx+y*yy+z*zz));

            if ( (a<tempW) && (a != 0))
            {
                W = a*std::exp(std::abs(w-w2))*std::exp(-(w+w2));
                tempW = a;
            }
        }
        sumW +=  W;
    }

    sum_W = sumW;

    for (size_t n=0; n<sum_u.size(); n++)
    {
        H[n] = sumW/(sumW+factor*sum_u.at(n));
    }
}

    void run()
    {
        // Get angles ==============================================================
        MetaData angles;
        angles.read(fnIn);
        size_t AngleNo = angles.size();
        if (AngleNo == 0 || !angles.containsLabel(MDL_ANGLE_ROT))
            REPORT_ERROR(ERR_MD_BADLABEL, "Input file doesn't contain angular information");

        double maxWeight = -99.e99;
        MultidimArray<double> weight;
        weight.initZeros(AngleNo);
        if (angles.containsLabel(MDL_WEIGHT))
        {
            // Find maximum weight
            int i=0;
            FOR_ALL_OBJECTS_IN_METADATA(angles)
            {
                double w;
                angles.getValue(MDL_WEIGHT,w,__iter.objId);
                DIRECT_A1D_ELEM(weight,i++)=w;
                maxWeight=XMIPP_MAX(w,maxWeight);
            }
        }

void ProgValidationNonTilt::obtainSumU_2(const MetaData & mdGallery, const MetaData & tempMd,std::vector<double> & sum_u,std::vector<double> & H0)
{

	const size_t tempMdSz= tempMd.size();
	const size_t numGallery= mdGallery.size();

    double xRan,yRan,zRan;
    size_t indx;
    double sumWRan;
    double * xRanArray = new double[tempMdSz];
    double * yRanArray = new double[tempMdSz];
    double * zRanArray  = new double[tempMdSz];
    std::vector<double> weightV;
    double a;

    double rot,tilt,w;
    bool mirror;

    for (size_t n=0; n<sum_u.size(); n++)
    {
        sumWRan = 0;

        for (size_t nS=0; nS<tempMd.size(); nS++)
        {
			indx = (size_t) (double( std::rand())*numGallery)/RAND_MAX;

			while ( (indx ==0) || (indx > numGallery) )
				indx = (size_t) (double( std::rand())*numGallery )/RAND_MAX;

        	mdGallery.getValue(MDL_ANGLE_ROT,rot,indx);
        	mdGallery.getValue(MDL_ANGLE_TILT,tilt,indx);
        	mdGallery.getValue(MDL_FLIP,mirror,indx);

            if (mirror == 1)
            	tilt = tilt + 180;

            xRan = sin(tilt*PI/180.)*cos(rot*PI/180.);
            yRan = sin(tilt*PI/180.)*sin(rot*PI/180.);
            zRan = std::abs(cos(tilt*PI/180.));

            xRanArray[nS] = xRan;
            yRanArray[nS] = yRan;
            zRanArray[nS] = zRan;
            //tempMd.getColumnValues(MDL_WEIGHT, weightV);
            tempMd.getColumnValues(MDL_MAXCC, weightV);
            std::random_shuffle(weightV.begin(), weightV.end());
        }


        sumWRan = 0;
        double WRan=0;
        double tempWRan, tempW1, tempW2, temp;
        for (size_t nS1=0; nS1<tempMd.size(); nS1++)
        {
            tempWRan = 1e3;
            for (size_t nS2=0; nS2<tempMd.size(); nS2++)
            {
            	temp = xRanArray[nS1]*xRanArray[nS2]+yRanArray[nS1]*yRanArray[nS2]+zRanArray[nS1]*zRanArray[nS2];
                if (temp < 1)
                	a = std::abs(std::acos(temp));
                else
                	a = 0;

                if ( (a<tempWRan) && (a > 0.00001) && (temp<1) )
                {
                    tempWRan = a;
                    tempW2 = weightV[nS2];
                    tempW1 = weightV[nS1];
                    WRan = a*std::exp(std::abs(tempW1-tempW2))*std::exp(-(tempW1+tempW2));
                    if (WRan == 0)
                    	WRan = a;
                }
            }
            sumWRan += WRan;
        }

        if (sumWRan == 0)
        	sumWRan = 0.075*tempMd.size();

        sum_u.at(n)=sumWRan;

    }

    size_t idx = 0;
    while (idx < sum_u.size())
    {
        std::random_shuffle(sum_u.begin(), sum_u.end());

        if(sum_u.at(0) != sum_u.at(1))
        {
            H0[idx] = sum_u.at(0)/(sum_u.at(0)+sum_u.at(1));
            idx += 1;
        }
    }

    delete xRanArray;
    delete yRanArray;
    delete zRanArray;

}

void ProgValidationNonTilt::obtainSumW(const MetaData & tempMd, double & sum_W, std::vector<double> & sum_u, std::vector<double> & H, const double factor)
{
    double a;
    double rot,tilt,w;
    double x,y,z;
    double xx,yy,zz;
    bool mirror;
    double w2;
    double tempW;
    double W=0;
    double temp;
    double sumW;

    sumW = 0;
    FOR_ALL_OBJECTS_IN_METADATA(tempMd)
    {
        tempMd.getValue(MDL_ANGLE_ROT,rot,__iter.objId);
        tempMd.getValue(MDL_ANGLE_TILT,tilt,__iter.objId);
        tempMd.getValue(MDL_FLIP,mirror,__iter.objId);
        if (mirror == 1)
        	tilt = tilt + 180;

        //tempMd.getValue(MDL_WEIGHT,w,__iter.objId);
        tempMd.getValue(MDL_MAXCC,w,__iter.objId);
        x = sin(tilt*PI/180.)*cos(rot*PI/180.);
        y = sin(tilt*PI/180.)*sin(rot*PI/180.);
        z = std::abs(cos(tilt*PI/180.));

        tempW = 1e3;
        _FOR_ALL_OBJECTS_IN_METADATA2(tempMd)
        {
            tempMd.getValue(MDL_ANGLE_ROT,rot,__iter2.objId);
            tempMd.getValue(MDL_ANGLE_TILT,tilt,__iter2.objId);
            //tempMd.getValue(MDL_WEIGHT,w2,__iter2.objId);
            tempMd.getValue(MDL_MAXCC,w2,__iter2.objId);
            tempMd.getValue(MDL_FLIP,mirror,__iter2.objId);
            if (mirror == 1)
            	tilt = tilt + 180;

            xx = sin(tilt*PI/180.)*cos(rot*PI/180.);
            yy = sin(tilt*PI/180.)*sin(rot*PI/180.);
            zz = std::abs(cos(tilt*PI/180.));
            temp = x*xx+y*yy+z*zz;

            if (temp < 1)
            	a = std::abs(std::acos(temp));
            else
            	a = 0;

            if ( (a<tempW) && (a > 0.00001) && (temp<1 ))
            {
                W = a*std::exp(std::abs(w-w2))*std::exp(-(w+w2));
                tempW = a;

                if (W == 0)
                	W = a;
            }
        }
        sumW +=  W;
    }

    //Here the problem is that maybe the changes are only in psi angle. We give the best solution assuming an angular sampling of 0.5º
    if (sumW == 0)
    	sumW = 0.075*tempMd.size();

    sum_W = sumW;
    for (size_t n=0; n<sum_u.size(); n++)
    {
        H[n] = sumW/(sumW+factor*sum_u.at(n));
    }
}

//.........这里部分代码省略.........
            for (size_t i = tIni; i <= tEnd; ++i)
            {
                objId = inMD.addObject();
                inMD.setValue(MDL_IMAGE, fnInput + formatString("/img%d.spe", i), objId);
            }
            break;
        }
    case GENERIC:
        {
            // Get Darkfield
            std::cerr << "Getting darkfield from "+fnInput << " ..." << std::endl;
            getDarkfield(fnInput, IavgDark);
            if (XSIZE(IavgDark())!=0)
                IavgDark.write(fnRoot+"_darkfield.xmp");


            std::vector<FileName> listDir;
            fnInput.getFiles(listDir);
            size_t objId;

            for (size_t i = 0; i < listDir.size(); ++i)
            {
                if (!listDir[i].hasImageExtension())
                    continue;
                objId = inMD.addObject();
                inMD.setValue(MDL_IMAGE, fnInput+"/"+listDir[i], objId);
            }
        }
        break;
    }

    inMD.findObjects(objIds);
    size_t nIm = inMD.size();

    // Create empty output stack file

    getImageInfo(inMD, imgInfo);


    /* Get the flatfield:: We get the FF after the image list because we need the image size to adapt the FF
     * in case they were already cropped.
     */
    if (!fnFlat.empty())
    {
        std::cout << "Getting flatfield from "+fnFlat << " ..." << std::endl;
        getFlatfield(fnFlat,IavgFlat);
        if ( XSIZE(IavgFlat()) != 0 )
        {
            FileName ffName = fnRoot+"_flatfield_avg.xmp";
            IavgFlat.write(ffName);
            fMD.setValue(MDL_IMAGE, ffName, fMD.addObject());
        }
    }

    createEmptyFile(fnOut, imgInfo.adim.xdim-cropSizeXi-cropSizeXe, imgInfo.adim.ydim-cropSizeYi-cropSizeYe, 1, nIm);

    // Process images
    td = new ThreadTaskDistributor(nIm, XMIPP_MAX(1, nIm/30));
    tm = new ThreadManager(thrNum, this);
    std::cerr << "Getting data from " << fnInput << " ...\n";
    init_progress_bar(nIm);
    tm->run(runThread);
    progress_bar(nIm);

    // Write Metadata and angles
    MetaData MDSorted;
    MDSorted.sort(outMD,MDL_ANGLE_TILT);
    MDSorted.write("[email protected]"+fnRoot + ".xmd");
    if ( fMD.size() > 0 )
        fMD.write("[email protected]"+fnRoot + ".xmd", MD_APPEND);

    // We also reference initial and final images at 0 degrees for Mistral tomograms
    if ( dSource == MISTRAL )
    {
        fMD.clear();
        FileName degree0Fn = "NXtomo/instrument/sample/0_degrees_initial_image";
        if ( H5File.checkDataset(degree0Fn.c_str()))
            fMD.setValue(MDL_IMAGE, degree0Fn + "@" + fnInput, fMD.addObject());
        degree0Fn = "NXtomo/instrument/sample/0_degrees_final_image";
        if ( H5File.checkDataset(degree0Fn.c_str()))
            fMD.setValue(MDL_IMAGE, degree0Fn + "@" + fnInput, fMD.addObject());
        if ( fMD.size() > 0 )
            fMD.write("[email protected]"+fnRoot + ".xmd", MD_APPEND);
    }

    // Write tlt file for IMOD
    std::ofstream fhTlt;
    fhTlt.open((fnRoot+".tlt").c_str());
    if (!fhTlt)
        REPORT_ERROR(ERR_IO_NOWRITE,fnRoot+".tlt");
    FOR_ALL_OBJECTS_IN_METADATA(MDSorted)
    {
        double tilt;
        MDSorted.getValue(MDL_ANGLE_TILT,tilt,__iter.objId);
        fhTlt << tilt << std::endl;
    }
    fhTlt.close();
    delete td;
    delete tm;
}

// Perform complete search =================================================
void ProgNmaAlignment::performCompleteSearch(const FileName &fnRandom,
		int pyramidLevel) const {
	String program;
	String arguments;
	const char * randStr = fnRandom.c_str();

	// Reduce the image
	FileName fnDown = formatString("%s_downimg.xmp", fnRandom.c_str());
	if (pyramidLevel != 0) {
		Image<double> I;
		I.read(currentImgName);
		selfPyramidReduce(BSPLINE3, I(), pyramidLevel);
		I.write(fnDown);
	} else
		copyImage(currentImgName.c_str(), fnDown.c_str());

	mkdir((fnRandom+"_ref").c_str(), S_IRWXU);

	double angSampling=2*RAD2DEG(atan(1.0/((double) imgSize / pow(2.0, (double) pyramidLevel+1))));
	angSampling=std::max(angSampling,discrAngStep);
	program = "xmipp_angular_project_library";
	arguments = formatString(
			"-i %s_deformedPDB.vol -o %s_ref/ref.stk --sampling_rate %f -v 0",
			randStr, randStr, angSampling);
	if (projMatch)
		arguments +=formatString(
						" --compute_neighbors --angular_distance -1 --experimental_images %s_downimg.xmp", randStr);

	runSystem(program, arguments, false);

	String refSelStr = formatString("%s_ref/ref.doc", randStr);

	if (fnmask != "") {
		program = "xmipp_transform_mask";
		arguments = formatString("-i %s --mask binary_file %s", refSelStr.c_str(),
				fnmask.c_str());
		runSystem(program, arguments, false);
	}

	// Perform alignment
	String fnOut=formatString("%s_angledisc.xmd",randStr);
	if (!projMatch) {
		program = "xmipp_angular_discrete_assign";
		arguments = formatString(
						"-i %s_downimg.xmp --ref %s -o %s --psi_step 5 --max_shift_change %d --search5D -v 0",
						randStr, refSelStr.c_str(), fnOut.c_str(), (int)round((double) imgSize / (10.0 * pow(2.0, (double) pyramidLevel))));
	} else {
		String refStkStr = formatString("%s_ref/ref.stk", randStr);
		program = "xmipp_angular_projection_matching";
		arguments =	formatString(
				        "-i %s_downimg.xmp --ref %s -o %s --search5d_step 1 --max_shift %d -v 0",
				        randStr, refStkStr.c_str(), fnOut.c_str(), (int)round((double) imgSize / (10.0 * pow(2.0, (double) pyramidLevel))));
	}
	runSystem(program, arguments, false);
	if (projMatch)
	{
		MetaData MD;
		MD.read(fnOut);
		bool flip;
		size_t id=MD.firstObject();
		MD.getValue(MDL_FLIP,flip,id);
		if (flip)
		{
			// This is because continuous assignment does not understand flips

			double shiftX, rot, tilt, psi, newrot, newtilt, newpsi;
			// Change sign in shiftX
			MD.getValue(MDL_SHIFT_X,shiftX,id);
			MD.setValue(MDL_SHIFT_X,-shiftX,id);

			// Change Euler angles
			MD.getValue(MDL_ANGLE_ROT,rot,id);
			MD.getValue(MDL_ANGLE_TILT,tilt,id);
			MD.getValue(MDL_ANGLE_PSI,psi,id);
			Euler_mirrorY(rot,tilt,psi,newrot,newtilt,newpsi);
			MD.setValue(MDL_ANGLE_ROT,newrot,id);
			MD.setValue(MDL_ANGLE_TILT,newtilt,id);
			MD.setValue(MDL_ANGLE_PSI,newpsi,id);
			MD.write(fnOut);
		}
	}
}

//majorAxis and minorAxis is the estimated particle size in px
void ProgSortByStatistics::processInprocessInputPrepareSPTH(MetaData &SF, bool trained)
{
    //#define DEBUG
    PCAMahalanobisAnalyzer tempPcaAnalyzer0;
    PCAMahalanobisAnalyzer tempPcaAnalyzer1;
    PCAMahalanobisAnalyzer tempPcaAnalyzer2;
    PCAMahalanobisAnalyzer tempPcaAnalyzer3;
    PCAMahalanobisAnalyzer tempPcaAnalyzer4;

    //Morphology
    tempPcaAnalyzer0.clear();
    //Signal to noise ratio
    tempPcaAnalyzer1.clear();
    tempPcaAnalyzer2.clear();
    tempPcaAnalyzer3.clear();
    //Histogram analysis, to detect black points and saturated parts
    tempPcaAnalyzer4.clear();

    double sign = 1;//;-1;
    int numNorm = 3;
    int numDescriptors0=numNorm;
    int numDescriptors2=4;
    int numDescriptors3=11;
    int numDescriptors4 = 10;

    MultidimArray<float> v0(numDescriptors0);
    MultidimArray<float> v2(numDescriptors2);
    MultidimArray<float> v3(numDescriptors3);
    MultidimArray<float> v4(numDescriptors4);

    if (verbose>0)
    {
        std::cout << " Sorting particle set by new xmipp method..." << std::endl;
    }

    int nr_imgs = SF.size();
    if (verbose>0)
        init_progress_bar(nr_imgs);

    int c = XMIPP_MAX(1, nr_imgs / 60);
    int imgno = 0, imgnoPCA=0;

    bool thereIsEnable=SF.containsLabel(MDL_ENABLED);
    bool first=true;

    // We assume that at least there is one particle
    size_t Xdim, Ydim, Zdim, Ndim;
    getImageSize(SF,Xdim,Ydim,Zdim,Ndim);

    //Initialization:
    MultidimArray<double> nI, modI, tempI, tempM, ROI;
    MultidimArray<bool> mask;
    nI.resizeNoCopy(Ydim,Xdim);
    modI.resizeNoCopy(Ydim,Xdim);
    tempI.resizeNoCopy(Ydim,Xdim);
    tempM.resizeNoCopy(Ydim,Xdim);
    mask.resizeNoCopy(Ydim,Xdim);
    mask.initConstant(true);

    MultidimArray<double> autoCorr(2*Ydim,2*Xdim);
    MultidimArray<double> smallAutoCorr;

    Histogram1D hist;
    Matrix2D<double> U,V,temp;
    Matrix1D<double> D;

    MultidimArray<int> radial_count;
    MultidimArray<double> radial_avg;
    Matrix1D<int> center(2);
    MultidimArray<int> distance;
    int dim;
    center.initZeros();

    v0.initZeros(numDescriptors0);
    v2.initZeros(numDescriptors2);
    v3.initZeros(numDescriptors3);
    v4.initZeros(numDescriptors4);

    ROI.resizeNoCopy(Ydim,Xdim);
    ROI.setXmippOrigin();
    FOR_ALL_ELEMENTS_IN_ARRAY2D(ROI)
    {
        double temp = std::sqrt(i*i+j*j);
        if ( temp < (Xdim/2))
            A2D_ELEM(ROI,i,j)= 1;
        else
            A2D_ELEM(ROI,i,j)= 0;
    }

    Image<double> img;
    FourierTransformer transformer(FFTW_BACKWARD);

    FOR_ALL_OBJECTS_IN_METADATA(SF)
    {
        if (thereIsEnable)
        {
            int enabled;
            SF.getValue(MDL_ENABLED,enabled,__iter.objId);
            if ( (enabled==-1)  )
//.........这里部分代码省略.........

void ProgSortByStatistics::processInputPrepare(MetaData &SF)
{
    PCAMahalanobisAnalyzer tempPcaAnalyzer;
    tempPcaAnalyzer.clear();

    Image<double> img;
    MultidimArray<double> img2;
    MultidimArray<int> radial_count;
    MultidimArray<double> radial_avg;
    Matrix1D<int> center(2);
    center.initZeros();

    if (verbose>0)
        std::cout << " Processing training set ..." << std::endl;

    int nr_imgs = SF.size();
    if (verbose>0)
        init_progress_bar(nr_imgs);
    int c = XMIPP_MAX(1, nr_imgs / 60);
    int imgno = 0, imgnoPCA=0;
    MultidimArray<float> v;
    MultidimArray<int> distance;
    int dim;

    bool thereIsEnable=SF.containsLabel(MDL_ENABLED);
    bool first=true;
    FOR_ALL_OBJECTS_IN_METADATA(SF)
    {
        if (thereIsEnable)
        {
            int enabled;
            SF.getValue(MDL_ENABLED,enabled,__iter.objId);
            if (enabled==-1)
                continue;
        }
        img.readApplyGeo(SF,__iter.objId);
        if (targetXdim!=-1 && targetXdim!=XSIZE(img()))
        	selfScaleToSize(LINEAR,img(),targetXdim,targetXdim,1);
        MultidimArray<double> &mI=img();
        mI.setXmippOrigin();
        mI.statisticsAdjust(0,1);

        // Overall statistics
        Histogram1D hist;
        compute_hist(mI,hist,-4,4,31);

        // Radial profile
        img2.resizeNoCopy(mI);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img2)
        {
            double val=DIRECT_MULTIDIM_ELEM(mI,n);
            DIRECT_MULTIDIM_ELEM(img2,n)=val*val;
        }
        if (first)
        {
            radialAveragePrecomputeDistance(img2, center, distance, dim);
            first=false;
        }
        fastRadialAverage(img2, distance, dim, radial_avg, radial_count);

        // Build vector
        v.initZeros(XSIZE(hist)+XSIZE(img2)/2);
        int idx=0;
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(hist)
        v(idx++)=(float)DIRECT_A1D_ELEM(hist,i);
        for (size_t i=0; i<XSIZE(img2)/2; i++)
            v(idx++)=(float)DIRECT_A1D_ELEM(radial_avg,i);

        tempPcaAnalyzer.addVector(v);

        if (imgno % c == 0 && verbose>0)
            progress_bar(imgno);
        imgno++;
        imgnoPCA++;
    }
    if (verbose>0)
        progress_bar(nr_imgs);

    MultidimArray<double> vavg,vstddev;
    tempPcaAnalyzer.computeStatistics(vavg,vstddev);
    tempPcaAnalyzer.evaluateZScore(2,20,false);
    pcaAnalyzer.insert(pcaAnalyzer.begin(), tempPcaAnalyzer);
}

    int main2()
    {

        MultidimArray<double> preImg, avgCurr, mappedImg;
        MultidimArray<double> outputMovie;
        Matrix1D<double> meanStdev;
        ImageGeneric movieStack;
        Image<double> II;
        MetaData MD; // To save plot information
        FileName motionInfFile, flowFileName, flowXFileName, flowYFileName;
        ArrayDim aDim;

        // For measuring times (both for whole process and for each level of the pyramid)
        clock_t tStart, tStart2;

#ifdef GPU
        // Matrix that we required in GPU part
        GpuMat d_flowx, d_flowy, d_dest;
        GpuMat d_avgcurr, d_preimg;
#endif

        // Matrix required by Opencv
        cv::Mat flow, dest, flowx, flowy;
        cv::Mat flowxPre, flowyPre;
        cv::Mat avgcurr, avgstep, preimg, preimg8, avgcurr8;
        cv::Mat planes[]={flowxPre, flowyPre};

        int imagenum, cnt=2, div=0, flowCounter;
        int h, w, levelNum, levelCounter=1;

        motionInfFile=foname.replaceExtension("xmd");
        std::string extension=fname.getExtension();
        if (extension=="mrc")
            fname+=":mrcs";
        movieStack.read(fname,HEADER);
        movieStack.getDimensions(aDim);
        imagenum = aDim.ndim;
        h = aDim.ydim;
        w = aDim.xdim;
        if (darkImageCorr)
        {
            II.read(darkRefFilename);
            darkImage=II();
        }
        if (gainImageCorr)
        {
            II.read(gianRefFilename);
            gainImage=II();
        }
        meanStdev.initZeros(4);
        //avgcurr=cv::Mat::zeros(h, w,CV_32FC1);
        avgCurr.initZeros(h, w);
        flowxPre=cv::Mat::zeros(h, w,CV_32FC1);
        flowyPre=cv::Mat::zeros(h, w,CV_32FC1);
#ifdef GPU

        // Object for optical flow
        FarnebackOpticalFlow d_calc;
        setDevice(gpuDevice);

        // Initialize the parameters for optical flow structure
        d_calc.numLevels=6;
        d_calc.pyrScale=0.5;
        d_calc.fastPyramids=true;
        d_calc.winSize=winSize;
        d_calc.numIters=1;
        d_calc.polyN=5;
        d_calc.polySigma=1.1;
        d_calc.flags=0;
#endif
        // Initialize the stack for the output movie
        if (saveCorrMovie)
            outputMovie.initZeros(imagenum, 1, h, w);
        // Correct for global motion from a cross-correlation based algorithms
        if (globalShiftCorr)
        {
            Matrix1D<double> shiftMatrix(2);
            shiftVector.reserve(imagenum);
            shiftMD.read(globalShiftFilename);
            FOR_ALL_OBJECTS_IN_METADATA(shiftMD)
            {
                shiftMD.getValue(MDL_SHIFT_X, XX(shiftMatrix), __iter.objId);
                shiftMD.getValue(MDL_SHIFT_Y, YY(shiftMatrix), __iter.objId);
                shiftVector.push_back(shiftMatrix);
            }
        }
        tStart2=clock();
        // Compute the average of the whole stack
        fstFrame++; // Just to adapt to Li algorithm
        lstFrame++; // Just to adapt to Li algorithm
        if (lstFrame>=imagenum || lstFrame==1)
            lstFrame=imagenum;
        imagenum=lstFrame-fstFrame+1;
        levelNum=sqrt(double(imagenum));
        computeAvg(fname, fstFrame, lstFrame, avgCurr);
        // if the user want to save the PSD
        if (doAverage)
        {
            II()=avgCurr;
            II.write(foname);
//.........这里部分代码省略.........