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

//#define DEBUG
void ProgSSNR::run()
{
    show();
    produceSideInfo();

    Matrix2D<double> output;
    if (!radial_avg)
    {
        if (!generate_VSSNR)
            estimateSSNR(1, output);
        else
            estimateSSNR(2, output);
        if (fn_out == "")
        	fn_out=fn_S.insertBeforeExtension("_SSNR").removeLastExtension().addExtension("xmd");
    }
    else
    {
        radialAverage(output);
        if (fn_out == "")
        	fn_out=fn_VSSNR.insertBeforeExtension("_radial_avg").removeLastExtension().addExtension("xmd");
    }
#ifdef DEBUG
    output.write(fn_out);
#endif
    MetaData MD;
    for (size_t i=1; i<MAT_YSIZE(output); ++i)
    {
    	size_t id=MD.addObject();
    	MD.setValue(MDL_RESOLUTION_FREQ,output(i,1),id);
  		MD.setValue(MDL_RESOLUTION_SSNR,output(i,2),id);
    	MD.setValue(MDL_RESOLUTION_FREQREAL,1.0/output(i,1),id);
    }
    MD.write(fn_out);
}

    void run()
    {
    	std::string extension=fnStack.getExtension();
    	if (extension=="mrc")
    		fnStack+=":mrcs";
    	else if (extension=="ali")
    		fnStack+=":mrcs";

    	std::ifstream fhAngles;
    	fhAngles.open(fnAngles.c_str());
    	if (!fhAngles)
    		REPORT_ERROR(ERR_IO_NOTOPEN,fnAngles);

    	Image<double> stack;
    	stack.read(fnStack, HEADER);
    	size_t xdim, ydim, zdim, ndim;
    	stack.getDimensions(xdim, ydim, zdim, ndim);

    	MetaData MD;
    	FileName fnSlice;
    	for (unsigned long slice = FIRST_IMAGE; slice <= ndim; slice++)
    	{
    		double angle;
    		fhAngles >> angle;
    		fnSlice.compose((int)slice,fnStack);
    		size_t objId=MD.addObject();
    		MD.setValue(MDL_IMAGE,fnSlice,objId);
    		MD.setValue(MDL_ANGLE_TILT,angle,objId);
    		MD.setValue(MDL_ENABLED,1,objId);
    	}
    	MD.write(fnOut);
    }

// Compute fitness =========================================================
double ObjFunc_nma_alignment::eval(Vector X, int *nerror) {
	int dim = global_nma_prog->numberOfModes;

	for (int i = 0; i < dim; i++) {
		global_nma_prog->trial(i) = X[i];
	}

	int pyramidLevelDisc = 1;
	int pyramidLevelCont = (global_nma_prog->currentStage == 1) ? 1 : 0;

	FileName fnRandom = global_nma_prog->createDeformedPDB(pyramidLevelCont);
	const char * randStr = fnRandom.c_str();

	if (global_nma_prog->currentStage == 1) {
		global_nma_prog->performCompleteSearch(fnRandom, pyramidLevelDisc);
	} else {
		double rot, tilt, psi, xshift, yshift;
		MetaData DF;

		rot = global_nma_prog->bestStage1(
				VEC_XSIZE(global_nma_prog->bestStage1) - 5);
		tilt = global_nma_prog->bestStage1(
				VEC_XSIZE(global_nma_prog->bestStage1) - 4);
		psi = global_nma_prog->bestStage1(
				VEC_XSIZE(global_nma_prog->bestStage1) - 3);
		xshift = global_nma_prog->bestStage1(
				VEC_XSIZE(global_nma_prog->bestStage1) - 2);
		yshift = global_nma_prog->bestStage1(
				VEC_XSIZE(global_nma_prog->bestStage1) - 1);

		size_t objId = DF.addObject();
		FileName fnDown = formatString("%s_downimg.xmp", randStr);
		DF.setValue(MDL_IMAGE, fnDown, objId);
		DF.setValue(MDL_ENABLED, 1, objId);
		DF.setValue(MDL_ANGLE_ROT, rot, objId);
		DF.setValue(MDL_ANGLE_TILT, tilt, objId);
		DF.setValue(MDL_ANGLE_PSI, psi, objId);
		DF.setValue(MDL_SHIFT_X, xshift, objId);
		DF.setValue(MDL_SHIFT_Y, yshift, objId);

		DF.write(formatString("%s_angledisc.xmd", randStr));
		copyImage(global_nma_prog->currentImgName.c_str(), fnDown.c_str());
	}
	double fitness = global_nma_prog->performContinuousAssignment(fnRandom,
			pyramidLevelCont);

	runSystem("rm", formatString("-rf %s* &", randStr));

	global_nma_prog->updateBestFit(fitness, dim);
	return fitness;
}

void ProgNmaAlignment::writeImageParameters(const FileName &fnImg) {
	MetaData md;
	size_t objId = md.addObject();
	md.setValue(MDL_IMAGE, fnImg, objId);
	md.setValue(MDL_ENABLED, 1, objId);
	md.setValue(MDL_ANGLE_ROT, parameters(0), objId);
	md.setValue(MDL_ANGLE_TILT, parameters(1), objId);
	md.setValue(MDL_ANGLE_PSI, parameters(2), objId);
	md.setValue(MDL_SHIFT_X, parameters(3), objId);
	md.setValue(MDL_SHIFT_Y, parameters(4), objId);

	int dim = numberOfModes;
	std::vector<double> vectortemp;

	for (int j = 5; j < 5 + dim; j++) {
		vectortemp.push_back(parameters(j));
	}

	md.setValue(MDL_NMA, vectortemp, objId);
	md.setValue(MDL_COST, parameters(5 + dim), objId);

	md.append(fnOutDir+"/nmaDone.xmd");
}

/* Save coordinates to disk ------------------------------------------------ */
void Micrograph::write_coordinates(int label, double minCost,
                                   const FileName &_fn_coords)
{
    std::ofstream fh;
    if (_fn_coords != "")
        fn_coords = _fn_coords;

    MetaData MD;
    MD.setComment((std::string) "Selected Coordinates for file " + fn_coords);
    int imax = coords.size();
    size_t id;
    for (int i = 0; i < imax; i++)
    {
        if (coords[i].valid && coords[i].cost > minCost
                && coords[i].label == label)
        {
            id = MD.addObject();
            MD.setValue(MDL_XCOOR, coords[i].X, id);
            MD.setValue(MDL_YCOOR, coords[i].Y, id);
        }
    }
    MD.write(fn_coords);
}

void runThread(ThreadArgument &thArg)
{
    int thread_id = thArg.thread_id;
    ProgXrayImport * ptrProg= (ProgXrayImport *)thArg.workClass;

    MetaData localMD;
    Image<double> Iaux;
    FileName fnImgIn, fnImgOut;
    size_t first = 0, last = 0;
    MultidimArray<char> mask;

    while (ptrProg->td->getTasks(first, last))
    {
        for (size_t i=first; i<=last; i++)
        {
            ptrProg->inMD.getValue(MDL_IMAGE, fnImgIn, ptrProg->objIds[i]);


            MDRow rowGeo;
            ptrProg->readGeoInfo(fnImgIn, rowGeo);
//            ptrProg->readAndCrop(fnImgIn, Iaux, ptrProg->cropSizeX, ptrProg->cropSizeY);

            Iaux.read(fnImgIn);
            Iaux().selfWindow(ptrProg->cropSizeYi,ptrProg->cropSizeXi,
            		(int)(YSIZE(Iaux())-ptrProg->cropSizeYe-1),(int)(XSIZE(Iaux())-ptrProg->cropSizeXe-1));

            Iaux().resetOrigin();




            if (XSIZE(ptrProg->IavgDark())!=0)
            {
                Iaux()-=ptrProg->IavgDark();
                forcePositive(Iaux());
            }


            double currentBeam = 1;
            double expTime = 1;
            double slitWidth = 1;

            if ( ptrProg->dSource == ptrProg->MISTRAL )
            {
                size_t idx = fnImgIn.getPrefixNumber();
                currentBeam = dMi(ptrProg->cBeamArray, idx-1);
                expTime = dMi(ptrProg->expTimeArray, idx-1);
                slitWidth = dMi(ptrProg->slitWidthArray, idx-1);
            }
            else
                ptrProg->readCorrectionInfo(fnImgIn, currentBeam, expTime, slitWidth);

            Iaux() *= 1.0/(currentBeam*expTime*slitWidth);
            if (XSIZE(ptrProg->IavgFlat())!=0)
                Iaux()/=ptrProg->IavgFlat();

            // Assign median filter to zero valued pixels to avoid -inf when applying log10
            Iaux().equal(0,mask);
            mask.resizeNoCopy(Iaux());

            if (XSIZE(ptrProg->bpMask()) != 0)
                mask += ptrProg->bpMask();

            boundMedianFilter(Iaux(), mask);

            if (ptrProg->logFix)
            {
                Iaux().selfLog();
                if (ptrProg->selfAttFix)
                    Iaux() *= -1.;
            }

            fnImgOut.compose(i+1, ptrProg->fnOut);


            size_t objId = localMD.addObject();
            localMD.setValue(MDL_IMAGE,fnImgOut,objId);
            localMD.setRow(rowGeo, objId); //
            //            localMD.setValue(MDL_ANGLE_TILT,Iaux.tilt(),objId);
            Iaux.write(fnImgOut);
            if (thread_id==0)
                progress_bar(i);
        }
    }
    //Lock for update the total counter
    ptrProg->mutex.lock();
    ptrProg->outMD.unionAll(localMD);
    ptrProg->mutex.unlock();
}

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

/* Produce all images ------------------------------------------------------ */
void Micrograph::produce_all_images(int label, double minCost,
                                    const FileName &fn_rootIn, const FileName &fn_image, double ang,
                                    double tilt, double psi, bool rmStack)
{
    MetaData SF;
    Image<double> I;
    Micrograph *M;

    // Set Source image
    if (fn_image == "")
        M = this;
    else
    {
        M = new Micrograph;
        M->open_micrograph(fn_image/*, swapbyte*/);
        M->set_window_size(X_window_size, Y_window_size);
        M->set_transmitance_flag(compute_transmitance);
        M->set_inverse_flag(compute_inverse);
    }

    // Set scale for particles
    int MXdim, MYdim, thisXdim, thisYdim;
    M->size(MXdim, MYdim);
    this->size(thisXdim, thisYdim);
    double scaleX = (double) MXdim / thisXdim;
    double scaleY = (double) MYdim / thisYdim;

    // Compute max and minimum if compute_transmitance
    // or compute_inverse flags are ON
    double Dmax=0., Dmin=0.;
    if (compute_transmitance || compute_inverse)
    {
        (*this).computeDoubleMinMax(Dmin, Dmax);

        if (compute_transmitance)
        {
            if (Dmin > 1)
                Dmin = log10(Dmin);
            if (Dmax > 1)
                Dmax = log10(Dmax);
        }
    }
    // Scissor all particles
    if (ang != 0)
        std::cout << "Angle from Y axis to tilt axis " << ang << std::endl
        << "   applying appropriate rotation\n";
    int nmax = ParticleNo();
    FileName fn_aux;
    FileName _ext = fn_rootIn.getFileFormat();
    FileName fn_out;
    FileName fn_root = fn_rootIn.removeFileFormat().removeLastExtension();
    if (fn_rootIn.hasStackExtension())
        fn_out=fn_root.addExtension(_ext);
    else
    	fn_out=fn_rootIn.addExtension("stk");

    if (rmStack)
        fn_out.deleteFile();
    size_t ii = 0;
    size_t id;
    for (int n = 0; n < nmax; n++)
        if (coords[n].valid && coords[n].cost > minCost && coords[n].label == label)
        {
            fn_aux.compose(++ii, fn_out);
            id = SF.addObject();
            // If the ctfRow was set, copy the info to images metadata
            if (ctfRow.containsLabel(MDL_CTF_DEFOCUSU))
                SF.setRow(ctfRow, id);
            SF.setValue(MDL_IMAGE, fn_aux, id);
            SF.setValue(MDL_MICROGRAPH, M->fn_micrograph, id);
            SF.setValue(MDL_XCOOR, coords[n].X, id);
            SF.setValue(MDL_YCOOR, coords[n].Y, id);
            bool t = M->scissor(coords[n], I(), Dmin, Dmax, scaleX, scaleY);
            if (!t)
            {
                std::cout << "Particle " << fn_aux
                << " is very near the border, "
                << "corresponding image is set to blank\n";
                SF.setValue(MDL_ENABLED, -1, id);
            }
            else
                SF.setValue(MDL_ENABLED, 1, id);
            //  if (ang!=0) I().rotate(-ang);
            I.write(fn_out, ii, true, WRITE_APPEND);
        }
    SF.write(fn_out.withoutExtension() + ".xmd");


    // Free source image??
    if (fn_image != "")
    {
        M->close_micrograph();
        delete M;
    }
}

//.........这里部分代码省略.........
                    readMat(flowXFileName.c_str(), flowx);
                    readMat(flowYFileName.c_str(), flowy);
                    d_flowx.upload(flowx);
                    d_flowy.upload(flowy);
                    d_calc.flags=cv::OPTFLOW_USE_INITIAL_FLOW;
                    d_calc(d_avgcurr, d_preimg, d_flowx, d_flowy);
                }
                d_flowx.download(planes[0]);
                d_flowy.download(planes[1]);
                d_avgcurr.release();
                d_preimg.release();
                d_flowx.release();
                d_flowy.release();
#else

                if (cnt==2)
                    calcOpticalFlowFarneback(avgcurr8, preimg8, flow, 0.5, 6, winSize, 1, 5, 1.1, 0);
                else
                {
                    flowFileName=foname.removeLastExtension()+formatString("flow%d%d.txt",div*2,flowCounter);
                    readMat(flowFileName.c_str(), flow);
                    calcOpticalFlowFarneback(avgcurr8, preimg8, flow, 0.5, 6, winSize, 1, 5, 1.1, cv::OPTFLOW_USE_INITIAL_FLOW);
                }
                split(flow, planes);

#endif
                // Save the flows if we are in the last step
                if (div==groupSize)
                {
                    if (i > 0)
                    {
                        std_dev2(planes,flowxPre,flowyPre,meanStdev);
                        size_t id=MD.addObject();
                        MD.setValue(MDL_OPTICALFLOW_MEANX, double(meanStdev(0)), id);
                        MD.setValue(MDL_OPTICALFLOW_MEANY, double(meanStdev(2)), id);
                        MD.setValue(MDL_OPTICALFLOW_STDX, double(meanStdev(1)), id);
                        MD.setValue(MDL_OPTICALFLOW_STDY, double(meanStdev(3)), id);
                        MD.write(motionInfFile, MD_APPEND);
                    }
                    planes[0].copyTo(flowxPre);
                    planes[1].copyTo(flowyPre);
                }
                else
                {
#ifdef GPU
                    flowXFileName=foname.removeLastExtension()+formatString("flowx%d%d.txt",div,i+1);
                    flowYFileName=foname.removeLastExtension()+formatString("flowy%d%d.txt",div,i+1);
                    saveMat(flowXFileName.c_str(), planes[0]);
                    saveMat(flowYFileName.c_str(), planes[1]);
#else

                    flowFileName=foname.removeLastExtension()+formatString("flow%d%d.txt",div,i+1);
                    saveMat(flowFileName.c_str(), flow);
#endif

                    if ((i+1)%2==0)
                        flowCounter++;
                }
                for( int row = 0; row < planes[0].rows; row++ )
                    for( int col = 0; col < planes[0].cols; col++ )
                    {
                        planes[0].at<float>(row,col) += (float)col;
                        planes[1].at<float>(row,col) += (float)row;
                    }
                cv::remap(preimg, dest, planes[0], planes[1], cv::INTER_CUBIC);
                if (div==1 && saveCorrMovie)