Inheritance diagram for ProgOpticalAligment:

Collaboration diagram for ProgOpticalAligment:

Public Member Functions
void	defineParams ()

void	readParams ()

void	saveMat (const FileName &fnOut, const cv::Mat &M)

void	readMat (const FileName &fnIn, cv::Mat &M)

void	xmipp2Opencv (const MultidimArray< float > &xmippArray, cv::Mat &opencvMat)

void	opencv2Xmipp (const cv::Mat &opencvMat, MultidimArray< float > &xmippArray)

void	convert2Uint8 (cv::Mat opencvDoubleMat, cv::Mat &opencvUintMat)

void	applyWindow (MultidimArray< double > &data)

int	computeAvg (size_t begin, size_t end, cv::Mat &cvAvgImg)

void	evaluateDisplacements (const cv::Mat flowCurrentGroup, const cv::Mat flowPreviousGroup, Matrix1D< double > &meanStdDev)

void	removeFlows (int level)

void	produceSideInfo ()

void	run ()

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Public Attributes
FileName	fnMovie

FileName	fnOut

FileName	fnGain

FileName	fnDark

FileName	fnMovieOut

FileName	fnMicOut

FileName	fnMovieUncOut

FileName	fnMicUncOut

FileName	fnMicInitial

MetaDataVec	movie

int	winSize

int	gpuDevice

int	nfirst

int	nlast

int	numberOfFrames

int	finalGroupSize

bool	globalShiftCorr

bool	inMemory

double	dose0

double	doseStep

double	accelerationVoltage

double	sampling

int	compensationMode

int	xLTcorner

int	yLTcorner

int	xDRcorner

int	yDRcorner

FileName	fnmovieRoot

FileName	fnTempStack

FileName	fnBaseOut

size_t	Xdim

size_t	Ydim

size_t	Zdim

size_t	Ndim

MultidimArray< float >	tempStack

Image< float >	tempAvg

FourierTransformer	transformer

ProgMovieFilterDose *	filterDose

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Detailed Description

Definition at line 62 of file movie_optical_alignment_cpu.cpp.

Member Function Documentation

◆ applyWindow()

void ProgOpticalAligment::applyWindow ( MultidimArray< double > & data )

inline

Definition at line 224 of file movie_optical_alignment_cpu.cpp.

     {
         if (yDRcorner!=-1)
             data.selfWindow(yLTcorner, xLTcorner, yDRcorner, xDRcorner);
     }

◆ computeAvg()

int ProgOpticalAligment::computeAvg	(	size_t	begin,
		size_t	end,
		cv::Mat &	cvAvgImg
	)

inline

Definition at line 231 of file movie_optical_alignment_cpu.cpp.

     {
         Image<float> frame;
         end=std::min(end,movie.size()-1);
         end=std::min(end,(size_t)nlast);
         if (begin>end)
             return 0;
 //      std::cout << "Computing average: " << begin << " " << end << std::endl;
 
         for (size_t i=begin;i<=end;i++)
         {
             int actualIdx=i-nfirst;
             if (inMemory)
                 frame().aliasImageInStack(tempStack,actualIdx);
             else
             {
 //                std::cout << "Reading " << fnTempStack << " " << actualIdx+1 << std::endl;
                 frame.read(fnTempStack,DATA,actualIdx+1);
             }
 
             if (i==begin)
                 tempAvg()=frame();
             else
                 tempAvg()+=frame();
         }
         tempAvg()/=float(end-begin+1);
 //        std::cout << "Avg stats: "; tempAvg().printStats(); std::cout << std::endl;
         xmipp2Opencv(tempAvg(), cvAvgImg);
         return end-begin+1;
     }

◆ convert2Uint8()

void ProgOpticalAligment::convert2Uint8	(	cv::Mat	opencvDoubleMat,
		cv::Mat &	opencvUintMat
	)

inline

Definition at line 216 of file movie_optical_alignment_cpu.cpp.

     {
         double min,max;
         cv::minMaxLoc(opencvDoubleMat, &min, &max);
         opencvDoubleMat.convertTo(opencvUintMat, CV_8U, 255.0/(max - min), -min * 255.0/(max - min));
     }

◆ defineParams()

void ProgOpticalAligment::defineParams ( )

inlinevirtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 95 of file movie_optical_alignment_cpu.cpp.

     {
         addUsageLine ("Align movies using optical flow");
         addParamsLine("    -i <inMoviewFnName>           : input movie File Name");
         addParamsLine("   [-o <fn=\"out.xmd\">]           : Metadata with the shifts of each frame.");
         addParamsLine("   [--oavg <fn=\"\">]              : Give the name of a micrograph to generate an aligned micrograph");
         addParamsLine("   [--oavgInitial <fn=\"\">]       : Use this option to save the initial average micrograph ");
         addParamsLine("                                   : before applying any alignment. ");
         addParamsLine("   [--oUnc <fnMic=\"\"> <fnMovie=\"\">] : Give the name of a micrograph and movie to generate an aligned and dose uncompensated micrograph");
         addParamsLine("   [--cropULCorner <x=0> <y=0>]    : crop up left corner (unit=px, index starts at 0)");
         addParamsLine("   [--cropDRCorner <x=-1> <y=-1>]  : crop down right corner (unit=px, index starts at 0), -1 -> no crop");
         addParamsLine("   [--frameRange <n0=-1> <nF=-1>]  : First and last frame to align, frame numbers start at 0");
         addParamsLine("   [--bin <s=-1>]               : Binning factor, it may be any floating number");
         addParamsLine("   [--winSize <int=150>]        : window size for optical flow algorithm");
         addParamsLine("   [--groupSize <int=1>]        : the depth of pyramid for optical flow algorithm");
         addParamsLine("   [--outMovie <fn=\"\">]       : save corrected, dose compensated stack");
         addParamsLine("   [--dark <fn=\"\">]           : Dark correction image");
         addParamsLine("   [--gain <fn=\"\">]           : Gain correction image (we will multiply by it)");
         addParamsLine("   [--inmemory]                 : Do not write a temporary file with the ");
         addParamsLine("   [--doseCorrection <dosePerFrame=0> <Ts=1> <kV=200> <previousDose=0> <mode=\"pre\">] : Set dosePerFrame to 0 if you do not want to correct by the dose");
         addParamsLine("                                : Dose in e/A^2, Sampling rate (Ts) in A. Valid modes are pre and post");
         addParamsLine("                                : Pre compensates before aligning and post after aligning");
 
 #ifdef GPU
         addParamsLine("   [--gpu <int=0>]              : GPU device to be used");
 #endif
     }

◆ evaluateDisplacements()

void ProgOpticalAligment::evaluateDisplacements	(	const cv::Mat *	flowCurrentGroup,
		const cv::Mat *	flowPreviousGroup,
		Matrix1D< double > &	meanStdDev
	)

inline

Definition at line 262 of file movie_optical_alignment_cpu.cpp.

     {
         double sumX=0, sumY=0;
         double absSumX=0, absSumY=0;
         double sqSumX=0, sqSumY=0;
         double valSubtract;
         int h=flowCurrentGroup[0].rows;
         int w=flowCurrentGroup[0].cols;
         for(int i=0;i<h;i++)
             for(int j=0;j<w;j++)
             {
                 valSubtract=flowCurrentGroup[0].at<float>(i,j)-flowPreviousGroup[0].at<float>(i,j);
                 sumX+=valSubtract;
                 absSumX+=abs(valSubtract);
                 sqSumX+=valSubtract*valSubtract;
                 valSubtract=flowCurrentGroup[1].at<float>(i,j)-flowPreviousGroup[1].at<float>(i,j);
                 sumY+=valSubtract;
                 absSumY+=abs(valSubtract);
                 sqSumY+=valSubtract*valSubtract;
             }
         double n=h*w;
         double avgX=sumX/n;
         double avgY=sumY/n;
         meanStdDev(0)=absSumX/n;
         meanStdDev(1)=sqrt(std::max(sqSumX/n-avgX*avgX,0.0));
         meanStdDev(2)=absSumY/n;
         meanStdDev(3)=sqrt(std::max(sqSumY/n-avgY*avgY,0.0));
     }

◆ opencv2Xmipp()

void ProgOpticalAligment::opencv2Xmipp	(	const cv::Mat &	opencvMat,
		MultidimArray< float > &	xmippArray
	)

inline

Definition at line 206 of file movie_optical_alignment_cpu.cpp.

     {
         int h = opencvMat.rows;
         int w = opencvMat.cols;
         xmippArray.resizeNoCopy(h, w);
         for (int row=0; row<h; ++row)
             memcpy(&DIRECT_A2D_ELEM(xmippArray,row,0),&opencvMat.at<float>(row,0),XSIZE(xmippArray)*sizeof(float));
     }

◆ produceSideInfo()

void ProgOpticalAligment::produceSideInfo ( )

inline

Definition at line 305 of file movie_optical_alignment_cpu.cpp.

     {
         fnmovieRoot=fnMovie.getDir();
         fnBaseOut=fnOut.removeDirectories();
 
         // If input is an stack create a metadata.
         if (fnMovie.isMetaData())
         {
             movie.read(fnMovie);
             getImageSize(movie, Xdim, Ydim, Zdim, Ndim);
             Ndim=movie.size();
         }
         else
         {
             ImageGeneric movieStack;
             movieStack.read(fnMovie,HEADER);
             movieStack.getDimensions(Xdim,Ydim,Zdim,Ndim);
             if (fnMovie.getExtension()=="mrc" and Ndim ==1)
                 Ndim = Zdim;
             size_t id;
             FileName fn;
             for (size_t i=0;i<Ndim;i++)
             {
                 id = movie.addObject();
                 fn.compose(i+FIRST_IMAGE,fnMovie);
                 movie.setValue(MDL_IMAGE, fn, id);
             }
         }
 
         // Get corners
         if (yDRcorner!=-1)
         {
             Xdim = xDRcorner - xLTcorner +1 ;
             Ydim = yDRcorner - yLTcorner +1 ;
         }
         if (Zdim!=1)
             REPORT_ERROR(ERR_ARG_INCORRECT,"This program is meant to align 2D frames, not 3D");
 
         // Dark and gain images
         Image<double> dark, igain;
         if (fnDark!="")
         {
             dark.read(fnDark);
             applyWindow(dark());
         }
 
         if (fnGain!="")
         {
             igain.read(fnGain);
             applyWindow(igain());
             // Pablo with coss: Consistency with other alignment methods
             // Do not divide by the gain, as the gain is expected to be inverted already
             double avg = igain().computeAvg();
             if (std::isinf(avg) || std::isnan(avg))
                 REPORT_ERROR(ERR_ARG_INCORRECT,
                              "The input gain image is incorrect, it contains infinite or nan");
         }
 
         // Count the number of frames to process
         numberOfFrames=Ndim;
         if (nfirst<0)
             nfirst=0;
         if (nlast<0)
             nlast=numberOfFrames-1;
         numberOfFrames=nlast-nfirst+1;
 
         // Initialize the stack for the output movie
         if (!fnMovieOut.isEmpty())
             createEmptyFile(fnMovieOut, Xdim, Ydim, 1, numberOfFrames, true, WRITE_REPLACE);
         if (!fnMovieUncOut.isEmpty())
             createEmptyFile(fnMovieUncOut, Xdim, Ydim, 1, numberOfFrames, true, WRITE_REPLACE);
 
         // Prepare stack
         if (inMemory)
             tempStack.resizeNoCopy(numberOfFrames,1,Ydim,Xdim);
         else
         {
             fnTempStack=fnmovieRoot+"tmpMovie.stk";
             createEmptyFile(fnTempStack,Xdim,Ydim,1,numberOfFrames,true);
         }
 
         int currentFrameInIdx=0,currentFrameOutIdx=0;
         FileName fnFrame;
         Image<double> frameImage, translatedImage;
         MultidimArray<float> Ifloat;
         Matrix1D<double> shift(2);
         filterDose=new ProgMovieFilterDose(accelerationVoltage);
         filterDose->pixel_size=sampling;
         MultidimArray< std::complex<double> > FFTI;
         for (size_t objId : movie.ids())
         {
             if (currentFrameInIdx>=nfirst && currentFrameInIdx<=nlast)
             {
                 movie.getValue(MDL_IMAGE, fnFrame, objId);
                 frameImage.read(fnFrame);
                 applyWindow(frameImage());
 
                 if (XSIZE(dark())>0)
                     frameImage()-=dark();
                 if (XSIZE(igain())>0)
                     frameImage()*=igain();
 
                 if (movie.containsLabel(MDL_SHIFT_X))
                 {
                     movie.getValue(MDL_SHIFT_X, XX(shift), objId);
                     movie.getValue(MDL_SHIFT_Y, YY(shift), objId);
                     if (fabs(XX(shift))>0 || fabs(YY(shift))>0)
                     {
 //                      std::cout << "Translating " << fnFrame << " by " << shift.transpose() << std::endl;
                         translate(xmipp_transformation::BSPLINE3, translatedImage(), frameImage(), shift, xmipp_transformation::WRAP);
                         frameImage()=translatedImage();
                     }
                 }
 
                 if (fnMovieUncOut!="")
                     frameImage.write(fnMovieUncOut, currentFrameOutIdx+1, true, WRITE_REPLACE);
 
                 if (compensationMode==PRECOMPENSATION)
                 {
                     transformer.FourierTransform(frameImage(), FFTI, false);
                     filterDose->applyDoseFilterToImage(YSIZE(frameImage()), XSIZE(frameImage()), FFTI,
                                                       dose0+currentFrameInIdx*doseStep, dose0+(currentFrameInIdx+1)*doseStep);
                     transformer.inverseFourierTransform();
                 }
 //                std::cout << "Reading " << fnFrame << " "; frameImage().printStats(); std::cout << std::endl;
                 if (inMemory)
                 {
                     typeCast(frameImage(),Ifloat);
                     memcpy(&DIRECT_NZYX_ELEM(tempStack,currentFrameOutIdx,0,0,0),&Ifloat(0,0),MULTIDIM_SIZE(Ifloat)*sizeof(float));
                 }
                 else
                 {
                     frameImage.write(fnTempStack, currentFrameOutIdx+1, true, WRITE_REPLACE);
 //                  std::cout << "Writing " << currentFrameOutIdx+1 << "@" << fnTempStack << std::endl;
                 }
                 currentFrameOutIdx++;
             }
             currentFrameInIdx++;
         }
     }

◆ readMat()

void ProgOpticalAligment::readMat	(	const FileName &	fnIn,
		cv::Mat &	M
	)

inline

Definition at line 183 of file movie_optical_alignment_cpu.cpp.

     {
         std::ifstream ifp(fnIn.c_str(), std::ios::in | std::ios::binary);
         int rows, cols;
         ifp.read(reinterpret_cast<char*>(&rows), sizeof(int));
         ifp.read(reinterpret_cast<char*>(&cols), sizeof(int));
         M.create(rows, cols, CV_32FC1);
         for (int row=0; row<rows; ++row)
                 ifp.read(reinterpret_cast<char*>(&(M.at<float>(row,0))), cols * sizeof(float));
         ifp.close();
     }

◆ readParams()

void ProgOpticalAligment::readParams ( )

inlinevirtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

Definition at line 123 of file movie_optical_alignment_cpu.cpp.

     {
         fnMovie = getParam("-i");
         fnOut = getParam("-o");
         fnGain = getParam("--gain");
         fnDark = getParam("--dark");
         fnMicOut = getParam("--oavg");
         fnMicInitial = getParam("--oavgInitial");
         fnMovieOut = getParam("--outMovie");
         if (checkParam("--oUnc"))
         {
             fnMicUncOut = getParam("--oUnc",0);
             fnMovieUncOut = getParam("--oUnc",1);
             if (fnMicUncOut!="" && fnMovieUncOut=="")
                 REPORT_ERROR(ERR_ARG_MISSING,"Writing the uncompensated images require both movie and micrograph");
         }
         finalGroupSize = getIntParam("--groupSize");
         nfirst = getIntParam("--frameRange",0);
         nlast = getIntParam("--frameRange",1);
         winSize   = getIntParam("--winSize");
         xLTcorner= getIntParam("--cropULCorner",0);
         yLTcorner= getIntParam("--cropULCorner",1);
         xDRcorner = getIntParam("--cropDRCorner",0);
         yDRcorner = getIntParam("--cropDRCorner",1);
         inMemory = checkParam("--inmemory");
         doseStep = getDoubleParam("--doseCorrection",0);
         sampling = getDoubleParam("--doseCorrection",1);
         accelerationVoltage = getDoubleParam("--doseCorrection",2);
         dose0 = getDoubleParam("--doseCorrection",3);
         String aux;
         aux=getParam("--doseCorrection",4);
         if (doseStep==0)
             compensationMode=NOCOMPENSATION;
         else
         {
             if (aux=="pre")
                 compensationMode=PRECOMPENSATION;
             else
                 compensationMode=POSTCOMPENSATION;
         }
 #ifdef GPU
         gpuDevice = getIntParam("--gpu");
 #endif
     }

◆ removeFlows()

void ProgOpticalAligment::removeFlows ( int level )

inline

Definition at line 291 of file movie_optical_alignment_cpu.cpp.

     {
         int g=0;
         FileName fnAuxX=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowx_%d_%d.flow",level,g);
         FileName fnAuxY=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowy_%d_%d.flow",level,g++);
         while (fileExists(fnAuxX))
         {
             deleteFile(fnAuxX);
             deleteFile(fnAuxY);
             fnAuxX=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowx_%d_%d.flow",level,g);
             fnAuxY=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowy_%d_%d.flow",level,g++);
         }
     }

◆ run()

void ProgOpticalAligment::run ( )

inlinevirtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

Definition at line 446 of file movie_optical_alignment_cpu.cpp.

     {
         produceSideInfo();
 
         Matrix1D<double> meanStdev;
         Image<float> undeformedGroupAverage, uncompensatedMic;
 
 #ifdef GPU
 #ifndef CV_VERSION_EPOCH // version 3 or newer
         // Matrices required in GPU part
         cv::cuda::GpuMat d_currentReference8, d_currentGroupAverage8;
 #else
     cv::gpu::GpuMat d_flowx, d_flowy, d_currentReference8, d_currentGroupAverage8;
 #endif 
 #else
         cv::Mat flow;
 #endif // GPU
 
         // Matrices required by Opencv
         cv::Mat cvCurrentReference, cvNewReference, cvCurrentGroupAverage, cvUndeformedGroupAverage, cvCurrentGroupAverage8, cvCurrentReference8;
         cv::Mat flowxCurrentGroup, flowyCurrentGroup, flowxPreviousGroup, flowyPreviousGroup;
         cv::Mat flowCurrentGroup[]={flowxCurrentGroup, flowyCurrentGroup};
         cv::Mat flowPreviousGroup[]={flowxPreviousGroup, flowyPreviousGroup};
 
         meanStdev.initZeros(4);
 #ifdef GPU
 #ifndef CV_VERSION_EPOCH // version 3 or newer
         cv::cuda::setDevice(gpuDevice);
     cv::Ptr<cv::cuda::FarnebackOpticalFlow> d_calc = cv::cuda::FarnebackOpticalFlow::create(
         6, // numLevels
         0.5, // pyrScale
         true, // fastPyramids
         winSize, // winSize
         1, // numIters
         5, // polyN
         1.1, // polySigma
         0); // flags
 #else
     cv::gpu::setDevice(gpuDevice);
         // Object for optical flow
         cv::gpu::FarnebackOpticalFlow d_calc;
         // 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
 #endif // GPU
 
         computeAvg(nfirst, nlast, cvCurrentReference);
         if (!fnMicInitial.isEmpty())
             tempAvg.write(fnMicInitial);
 
         cout << "Frames " << nfirst << " to " << nlast << " under processing ..." << std::endl;
 
         int numberOfGroups=2;
         int levelNum=int(ceil(log(double(numberOfFrames)/finalGroupSize)/log(2.0))), levelCounter=0;
         MetaDataVec MDout; // To save plot information
         while (levelCounter<levelNum)
         {
             convert2Uint8(cvCurrentReference,cvCurrentReference8);
     #ifdef GPU
             d_currentReference8.upload(cvCurrentReference8);
     #endif
             int currentGroupSize=int(ceil((float)numberOfFrames/numberOfGroups));
 
             bool lastLevel = levelCounter==levelNum-1;
             // avgStep to hold the sum of aligned frames of each group at each step
             cvNewReference=cv::Mat::zeros(Ydim, Xdim, CV_32FC1);
 
             cout << "Level " << levelCounter << "/" << levelNum-1
                  << " of the pyramid is under processing" << std::endl;
             cout << "   Group size: " << currentGroupSize << ", Number of groups: " << numberOfGroups << std::endl;
 
             // Compute time for each level
             clock_t tStart = clock();
 
             for (int currentGroup=0; currentGroup<numberOfGroups; currentGroup++)
             {
                 int NimgsInAvg=computeAvg(currentGroup*currentGroupSize+nfirst, (currentGroup+1)*currentGroupSize+nfirst-1, cvCurrentGroupAverage);
                 if (NimgsInAvg==0)
                     continue;
                 convert2Uint8(cvCurrentGroupAverage,cvCurrentGroupAverage8);
 
                 if (numberOfGroups>2)
                 {
                     FileName flowXFileName=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowx_%d_%d.flow",levelCounter-1,currentGroup/2);
                     readMat(flowXFileName.c_str(), flowCurrentGroup[0]);
 //                    std::cout << "Reading flow " << flowXFileName << std::endl;
 
                     FileName flowYFileName=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowy_%d_%d.flow",levelCounter-1,currentGroup/2); // ***
                     readMat(flowYFileName.c_str(), flowCurrentGroup[1]);
                 }
 
 #ifdef GPU
                 d_currentGroupAverage8.upload(cvCurrentGroupAverage8);
 #ifndef CV_VERSION_EPOCH // version 3 or newer
         cv::cuda::GpuMat d_flow(cvCurrentGroupAverage8.size(), CV_32FC2);
         cv::cuda::GpuMat planes[2];
         if (numberOfGroups>2)
                 {
                 cv::cuda::split(d_flow, planes);
                     planes[0].upload(flowCurrentGroup[0]);
                     planes[1].upload(flowCurrentGroup[1]);
                     d_calc->setFlags(cv::OPTFLOW_USE_INITIAL_FLOW);
                 }
                 d_calc->calc(d_currentReference8, d_currentGroupAverage8, d_flow);
 
                 cv::cuda::split(d_flow, planes);
                 planes[0].download(flowCurrentGroup[0]);
                 planes[1].download(flowCurrentGroup[1]);
 #else
                 if (numberOfGroups>2)
                 {
                     d_flowx.upload(flowCurrentGroup[0]);
                     d_flowy.upload(flowCurrentGroup[1]);
                     d_calc.flags=cv::OPTFLOW_USE_INITIAL_FLOW;
                 }
                 d_calc(d_currentReference8, d_currentGroupAverage8, d_flowx, d_flowy);
                 d_flowx.download(flowCurrentGroup[0]);
                 d_flowy.download(flowCurrentGroup[1]);
                 d_currentGroupAverage8.release();
                 d_flowx.release();
                 d_flowy.release();
 #endif
 #else // CPU version
                 int ofFlags=0;
                 // Check if we should use the flows from the previous steps
                 if (numberOfGroups==2)
                 {
                     flowxCurrentGroup=cv::Mat::zeros(Ydim, Xdim, CV_32FC1);
                     flowyCurrentGroup=cv::Mat::zeros(Ydim, Xdim, CV_32FC1);
                     flow=cv::Mat::zeros(Ydim, Xdim, CV_32FC1);
                 }
                 else
                     ofFlags=cv::OPTFLOW_USE_INITIAL_FLOW;
 
                 merge(flowCurrentGroup,2,flow);
 
                 calcOpticalFlowFarneback(cvCurrentReference8, cvCurrentGroupAverage8, flow, 0.5, 6, winSize, 1, 5, 1.1, ofFlags);
                 split(flow, flowCurrentGroup);
 #endif
 
                 // Save the flows if we are in the last step
                 if (lastLevel)
                 {
                     if (currentGroup > 0)
                     {
                         evaluateDisplacements(flowCurrentGroup,flowPreviousGroup,meanStdev);
                         size_t id=MDout.addObject();
                         MDout.setValue(MDL_OPTICALFLOW_MEANX, double(meanStdev(0)), id);
                         MDout.setValue(MDL_OPTICALFLOW_MEANY, double(meanStdev(2)), id);
                         MDout.setValue(MDL_OPTICALFLOW_STDX, double(meanStdev(1)), id);
                         MDout.setValue(MDL_OPTICALFLOW_STDY, double(meanStdev(3)), id);
                     }
                     flowCurrentGroup[0].copyTo(flowPreviousGroup[0]);
                     flowCurrentGroup[1].copyTo(flowPreviousGroup[1]);
                 }
                 else
                 {
                     FileName flowXFileName=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowx_%d_%d.flow",levelCounter,currentGroup);
                     FileName flowYFileName=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowy_%d_%d.flow",levelCounter,currentGroup);
                     saveMat(flowXFileName.c_str(), flowCurrentGroup[0]);
                     saveMat(flowYFileName.c_str(), flowCurrentGroup[1]);
 //                  flowXFileName=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowx_%d_%d.mrc",levelCounter,currentGroup);
 //                  flowYFileName=fnmovieRoot+fnBaseOut.removeLastExtension()+formatString("flowy_%d_%d.mrc",levelCounter,currentGroup);
 //                  std::cout << "Saving flows " << flowXFileName << " " << flowYFileName << std::endl;
 //                  Image<float> save;
 //                  opencv2Xmipp(flowCurrentGroup[0],save());
 //                  save.write(flowXFileName);
 //                  opencv2Xmipp(flowCurrentGroup[1],save());
 //                  save.write(flowYFileName);
                 }
 
                 // Update cvNewReference
                 for (float row = 0; row < flowCurrentGroup[0].rows; row++ )
                     for (float col = 0; col < flowCurrentGroup[0].cols; col++ )
                     {
                         flowCurrentGroup[0].at<float>(row,col) += col;
                         flowCurrentGroup[1].at<float>(row,col) += row;
                     }
                 cv::remap(cvCurrentGroupAverage, cvUndeformedGroupAverage, flowCurrentGroup[0], flowCurrentGroup[1], cv::INTER_CUBIC);
 //                Image<float> save;
 //                opencv2Xmipp(cvUndeformedGroupAverage,save());
 //                std::cout << "cvUndeformedGroupAverage stats "; save().printStats(); std::cout << std::endl;
 //                opencv2Xmipp(flowCurrentGroup[0],save());
 //                std::cout << "flowCurrentGroup[0] stats "; save().printStats(); std::cout << std::endl;
 //                opencv2Xmipp(flowCurrentGroup[1],save());
 //                std::cout << "flowCurrentGroup[1] stats "; save().printStats(); std::cout << std::endl;
                 if (lastLevel)
                 {
                     if (compensationMode==POSTCOMPENSATION)
                     {
                         MultidimArray<float> If;
                         MultidimArray<double> Id;
                         MultidimArray<std::complex<double> > fId;
                         opencv2Xmipp(cvUndeformedGroupAverage,If);
                         typeCast(If,Id);
                         transformer.FourierTransform(Id, fId, false);
                         filterDose->applyDoseFilterToImage(YSIZE(Id), XSIZE(Id), fId,
                                                           dose0+currentGroup*doseStep, dose0+(currentGroup+1)*doseStep);
                         transformer.inverseFourierTransform();
                         typeCast(Id,If);
                         xmipp2Opencv(If,cvUndeformedGroupAverage);
                     }
                     if (!fnMovieOut.isEmpty())
                     {
                         opencv2Xmipp(cvUndeformedGroupAverage,undeformedGroupAverage());
                         // std::cout << "Writing frame " << fnMovieOut << " " << currentGroup+1 << std::endl;
                         undeformedGroupAverage.write(fnMovieOut, currentGroup+1, true, WRITE_REPLACE);
                     }
                     if (!fnMovieUncOut.isEmpty())
                     {
                         Image<float> frame;
                         cv::Mat cvFrame, cvUndeformedFrame;
                         frame.read(formatString("%d@%s",currentGroup+1,fnMovieUncOut.c_str()));
                         xmipp2Opencv(frame(), cvFrame);
                         cv::remap(cvFrame, cvUndeformedFrame, flowCurrentGroup[0], flowCurrentGroup[1], cv::INTER_CUBIC);
                         opencv2Xmipp(cvUndeformedFrame,frame());
                         frame.write(fnMovieUncOut, currentGroup+1, true, WRITE_REPLACE);
 
                         if (currentGroup==0)
                             uncompensatedMic()=frame();
                         else
                             uncompensatedMic()+=frame();
                     }
                 }
                 cvNewReference+=cvUndeformedGroupAverage;
             }
 #ifdef GPU
             d_currentReference8.release();
 #endif
             cvCurrentReference=cvNewReference;
             cvCurrentReference*=1.0/numberOfGroups;
 //            Image<float> save;
 //            opencv2Xmipp(cvCurrentReference,save());
 //            std::cout << "cvCurrentReference stats "; save().printStats(); std::cout << std::endl;
             printf("Processing time: %.2fs\n", (double)(clock() - tStart)/CLOCKS_PER_SEC);
             numberOfGroups=std::min(2*numberOfGroups,numberOfFrames);
             levelCounter++;
             if (levelCounter>1)
                 removeFlows(levelCounter-2);
         }
         MDout.write(fnOut, MD_APPEND);
         removeFlows(levelNum-1);
         removeFlows(levelNum);
 
         if (!fnMicOut.isEmpty())
         {
             opencv2Xmipp(cvCurrentReference, tempAvg());
             tempAvg.write(fnMicOut);
         }
 
         if (!fnMicUncOut.isEmpty())
         {
             uncompensatedMic.write(fnMicUncOut);
 //          std::cout << "Writing " << fnMicUncOut << std::endl;
         }
 
         if (!inMemory)
             deleteFile(fnTempStack);
     }

◆ saveMat()

void ProgOpticalAligment::saveMat	(	const FileName &	fnOut,
		const cv::Mat &	M
	)

inline

Definition at line 169 of file movie_optical_alignment_cpu.cpp.

     {
         std::ofstream ofp(fnOut.c_str(), std::ios::out | std::ios::binary);
         ofp.write(reinterpret_cast<const char*>(&M.rows), sizeof(int));
         ofp.write(reinterpret_cast<const char*>(&M.cols), sizeof(int));
 //  ofp << M.rows;
 //  ofp << M.cols;
     for (int row=0; row<M.rows; ++row) {
                 ofp.write(reinterpret_cast<const char*>(&M.at<float>(row,0)), M.cols * sizeof(float));
     }
         ofp.close();
     }

◆ xmipp2Opencv()

void ProgOpticalAligment::xmipp2Opencv	(	const MultidimArray< float > &	xmippArray,
		cv::Mat &	opencvMat
	)

inline

Definition at line 196 of file movie_optical_alignment_cpu.cpp.

     {
         int h = YSIZE(xmippArray);
         int w = XSIZE(xmippArray);
         opencvMat = cv::Mat::zeros(h, w,CV_32FC1);
         for (int row=0; row<h; ++row)
             memcpy(&opencvMat.at<float>(row,0),&DIRECT_A2D_ELEM(xmippArray,row,0),XSIZE(xmippArray)*sizeof(float));
     }