#include <xmipp_gpu_correlation.h>

Inheritance diagram for ProgGpuCorrelation:

Collaboration diagram for ProgGpuCorrelation:

Public Member Functions
void	readParams ()
	Read argument from command line. More...

void	show ()
	Show. More...

void	defineParams ()
	Define parameters. More...

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
MetaDataVec	SF

MetaDataVec	SFexp

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 36 of file xmipp_gpu_correlation.h.

Member Function Documentation

◆ defineParams()

void ProgGpuCorrelation::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippProgram.

Definition at line 536 of file xmipp_gpu_correlation.cpp.

 {
 
     addParamsLine("   -i_ref  <md_ref_file>                : Metadata file with input reference images");
     addParamsLine("   -i_exp  <md_exp_file>                : Metadata file with input experimental images");
     addParamsLine("   -o      <md_out>                     : Output metadata file");
     addParamsLine("   [--classify <md_classes_out=\"output_classes.xmd\">]         : To generate the aligned output images and write the associated metadata");
     addParamsLine("   [--keep_best <N=2>]                  : To keep N aligned images with the highest correlation");
     addParamsLine("   [--significance <alpha=0.2>]         : To use significance with the indicated value");
     addParamsLine("   [--odir <outputDir=\".\">]           : Output directory to save the aligned images");
     addParamsLine("   [--maxShift <s=10>]                  : Maximum shift allowed (+-this amount)");
     addParamsLine("   [--simplifiedMd <b=false>]           : To generate a simplified metadata with only the maximum weight image stores");
     addParamsLine("   [--sizePad <pad=100>]    ");
     addParamsLine("   [--device <dev=0>]                   : GPU device to use. 0th by default");
     addUsageLine("Computes the correlation between a set of experimental images with respect "
                  "to a set of reference images with CUDA in GPU");
 
 }

◆ readParams()

void ProgGpuCorrelation::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippProgram.

Definition at line 494 of file xmipp_gpu_correlation.cpp.

 {
 
     fn_ref = getParam("-i_ref");
     fn_exp = getParam("-i_exp");
     fn_out = getParam("-o");
     generate_out = checkParam("--classify");
     fn_classes_out = getParam("--classify");
     significance = checkParam("--significance");
     simplifiedMd = checkParam("--simplifiedMd");
     if(significance){
         alpha=getDoubleParam("--significance");
         keepN=false;
     }
     if(checkParam("--keep_best") && !significance){
         keepN=true;
         n_keep=getIntParam("--keep_best");
     }
     if(!keepN && !significance){
         keepN=true;
         n_keep=getIntParam("--keep_best");
     }
     fnDir = getParam("--odir");
     maxShift = getIntParam("--maxShift");
     sizePad = getIntParam("--sizePad");
     int device = getIntParam("--device");
     gpu = GPU(device);
 
 }

◆ run()

void ProgGpuCorrelation::run ( )

virtual

processImage

Reimplemented from XmippProgram.

Definition at line 1337 of file xmipp_gpu_correlation.cpp.

 {
 
     //Setting the cuda device to use
     gpu.set();
 
     //PROJECTION IMAGES
     size_t Xdim, Ydim, Zdim, Ndim;
     SF.read(fn_ref,NULL);
     size_t mdInSize = SF.size();
     getImageSize(SF, Xdim, Ydim, Zdim, Ndim);
 
 
     //EXPERIMENTAL IMAGES
     SFexp.read(fn_exp,NULL);
     size_t mdExpSize = SFexp.size();
 
     // Generate mask
     Mask mask;
     mask.type = BINARY_CIRCULAR_MASK;
     mask.mode = INNER_MASK;
     auto rad = (size_t)std::min(Xdim*0.48, Ydim*0.48);
 
     int number = rad;
     auto *out = new int[5];
 
     while(true){
         if (number%2!=0){
             number--;
             continue;
         }
         for (int z=0; z<5; z++)
             out[z]=0;
         primeFactors(number, out);
         if ((out[0]!=0 || out[1]!=0 || out[2]!=0 || out[3]!=0) && out[4]==0){
             rad = number;
             break;
         }
         else
             number--;
     }
 
     mask.R1 = rad;
     mask.resize(Ydim,Xdim);
     mask.get_binary_mask().setXmippOrigin();
     mask.generate_mask();
     int maskCount = mask.get_binary_mask().sum();
 
     //AJ check the size of the data to avoid exceed the GPU memory
     float memory[3]={0, 0, 0}; //total, free, used
     cuda_check_gpu_memory(memory);
 
     int maxGridSize[3];
     cuda_check_gpu_properties(maxGridSize);
 
 
     //AJ check_gpu_memory to know how many images we can copy in the gpu memory
     float limit=0.4; //0.877; 1.3;
     int available_images_proj = mdExpSize; //mdInSize
     int available1 = mdExpSize;
     int available2 = mdExpSize;
     if(Xdim*Ydim*mdExpSize*4*100/memory[1]>limit){ //mdInSize
         available1 = floor(memory[1]*(limit/100)/(Xdim*Ydim*4));
     }
     if(Xdim*2*Ydim*2*mdExpSize>maxGridSize[0]){ //mdInSize
         available2 = floor((round(maxGridSize[0]*0.9))/(Xdim*Ydim*2*2));
     }
     if (available1<available2)
         available_images_proj = available1;
     else
         available_images_proj = available2;
 
 
     //matrix with all the best transformations in CPU
     auto *matrixTransCpu = new MultidimArray<float> [mdInSize]; //mdExpSize
     for(int i=0; i<mdInSize; i++) //mdExpSize
         matrixTransCpu[i].coreAllocate(1, mdExpSize, 3, 3); //mdInSize
     auto *matrixTransCpu_mirror = new MultidimArray<float> [mdInSize]; //mdExpSize
     for(int i=0; i<mdInSize; i++) //mdExpSize
         matrixTransCpu_mirror[i].coreAllocate(1, mdExpSize, 3, 3); //mdInSize
 
     //correlation matrix
     MultidimArray<float> matrixCorrCpu(1, 1, mdInSize, mdExpSize); //mdExpSize, mdInSize
     MultidimArray<float> matrixCorrCpu_mirror(1, 1, mdInSize, mdExpSize); //mdExpSize, mdInSize
 
     //Aux vectors with maximum values of correlation in RT and TR steps
     float *max_vector_rt;
     float *max_vector_tr;
     float *max_vector_rt_mirror;
     float *max_vector_tr_mirror;
 
     //Transformation matrix in GPU and CPU
     TransformMatrix<float> transMat_tr;
     TransformMatrix<float> transMat_rt;
     TransformMatrix<float> transMat_tr_mirror;
     TransformMatrix<float> transMat_rt_mirror;
 
     TransformMatrix<float> resultTR;
     TransformMatrix<float> resultRT;
 
     int firstIdx=0;
     bool finish=false;
 
     mycufftHandle myhandlePadded_tr, myhandleMask_tr, myhandlePolar_tr, myhandleAux_tr, myhandlePaddedB_tr, myhandleMaskB_tr, myhandlePolarB_tr, myhandleAuxB_tr;
     mycufftHandle myhandlePadded_rt, myhandleMask_rt, myhandlePolar_rt, myhandleAux_rt, myhandlePaddedB_rt, myhandleMaskB_rt, myhandlePolarB_rt, myhandleAuxB_rt;
     mycufftHandle ifftcb;
 
     myStreamHandle myStreamTR, myStreamRT;
     myStreamCreate(myStreamTR);
     myStreamCreate(myStreamRT);
 
 
     GpuCorrelationAux d_referenceAux;
 
     size_t pad_Xdim=2*Xdim-1;
     size_t pad_Ydim=2*Ydim-1;
 
     number = pad_Xdim;
     while(true){
         if (number%2!=0){
             number++;
             continue;
         }
         for (int z=0; z<5; z++)
             out[z]=0;
         primeFactors(number, out);
         if ((out[0]!=0 || out[1]!=0 || out[2]!=0 || out[3]!=0) && out[4]==0){
             pad_Xdim = number;
             break;
         }
         else
             number++;
     }
 
     pad_Ydim = pad_Xdim;
     d_referenceAux.XdimOrig=Xdim;
     d_referenceAux.YdimOrig=Ydim;
     d_referenceAux.Xdim=pad_Xdim;
     d_referenceAux.Ydim=pad_Ydim;
     d_referenceAux.XdimPolar=360;
     d_referenceAux.YdimPolar=(size_t)mask.R1;
 
 
     StructuresAux myStructureAux_tr, myStructureAux_rt;
 
     auto iter = SFexp.ids().begin();
 
     GpuMultidimArrayAtCpu<float> original_image_stack;
 
     //Loop over the reference images
     size_t totalWork=mdInSize*mdExpSize;
     size_t workDone=0;
     init_progress_bar(totalWork);
     size_t lastProgressShown=0;
     while(!finish){
 
         original_image_stack.resize(Xdim,Ydim,1,available_images_proj);
 
         //Aux vectors with maximum values of correlation in RT and TR steps
         cpuMalloc((void**)&max_vector_tr, sizeof(float)*available_images_proj);
         cpuMalloc((void**)&max_vector_rt, sizeof(float)*available_images_proj);
         cpuMalloc((void**)&max_vector_tr_mirror, sizeof(float)*available_images_proj);
         cpuMalloc((void**)&max_vector_rt_mirror, sizeof(float)*available_images_proj);
 
 
         //Transformation matrix in GPU and CPU
         transMat_tr.resize(myStreamTR, available_images_proj);
         transMat_rt.resize(myStreamRT, available_images_proj);
         transMat_tr_mirror.resize(myStreamTR, available_images_proj);
         transMat_rt_mirror.resize(myStreamRT, available_images_proj);
 
         resultTR.resize(myStreamTR, available_images_proj);
         resultRT.resize(myStreamRT, available_images_proj);
 
         //TODO allocate memory with care
         myStructureAux_tr.padded_image_gpu.resize(pad_Xdim, pad_Ydim, 1, available_images_proj);
         myStructureAux_tr.padded_image2_gpu.resize(pad_Xdim, pad_Ydim, 1, available_images_proj);
         myStructureAux_tr.padded_mask_gpu.resize(pad_Xdim, pad_Ydim, 1, 1);
         myStructureAux_tr.polar_gpu.resize(d_referenceAux.XdimPolar,d_referenceAux.YdimPolar,1,available_images_proj);
         myStructureAux_tr.polar2_gpu.resize(d_referenceAux.XdimPolar,d_referenceAux.YdimPolar,1,available_images_proj);
 
         myStructureAux_rt.padded_image_gpu.resize(pad_Xdim, pad_Ydim, 1, available_images_proj);
         myStructureAux_rt.padded_image2_gpu.resize(pad_Xdim, pad_Ydim, 1, available_images_proj);
         myStructureAux_rt.padded_mask_gpu.resize(pad_Xdim, pad_Ydim, 1, 1);
         myStructureAux_rt.polar_gpu.resize(d_referenceAux.XdimPolar,d_referenceAux.YdimPolar,1,available_images_proj);
         myStructureAux_rt.polar2_gpu.resize(d_referenceAux.XdimPolar,d_referenceAux.YdimPolar,1,available_images_proj);
 
         //SF
         preprocess_images_reference(SFexp, firstIdx, available_images_proj, mask, d_referenceAux,
                 myhandlePadded_tr, myhandleMask_tr, myhandlePolar_tr,  myStructureAux_tr, iter, myStreamTR);
 
         d_referenceAux.maskCount=maskCount;
         d_referenceAux.produceSideInfo(myhandlePaddedB_tr, myhandleMaskB_tr, myStructureAux_tr, myStreamTR);
 
         //AJ calling a cudaDeviceSyncrhonize to be sure that these images are loaded in gpu memory
         // and available for all the streams
         waitGpu(myStreamTR, true);
 
         //EXPERIMENTAL IMAGES PART
         size_t expIndex = 0;
         FileName fnImgExp;
         auto iterExp = SF.begin();
 
         GpuCorrelationAux d_experimentalAuxTR, d_experimentalAuxRT;
         d_experimentalAuxTR.XdimOrig=d_referenceAux.XdimOrig;
         d_experimentalAuxTR.YdimOrig=d_referenceAux.YdimOrig;
         d_experimentalAuxTR.Xdim=d_referenceAux.Xdim;
         d_experimentalAuxTR.Ydim=d_referenceAux.Ydim;
         d_experimentalAuxTR.XdimPolar=d_referenceAux.XdimPolar;
         d_experimentalAuxTR.YdimPolar=d_referenceAux.YdimPolar;
 
         d_experimentalAuxRT.XdimOrig=d_referenceAux.XdimOrig;
         d_experimentalAuxRT.YdimOrig=d_referenceAux.YdimOrig;
         d_experimentalAuxRT.Xdim=d_referenceAux.Xdim;
         d_experimentalAuxRT.Ydim=d_referenceAux.Ydim;
         d_experimentalAuxRT.XdimPolar=d_referenceAux.XdimPolar;
         d_experimentalAuxRT.YdimPolar=d_referenceAux.YdimPolar;
 
         //TODO: here we can use threads to carry out the alignment of different images in different threads
         size_t n=0;
         int available_images_exp = mdInSize; //mdExpSize
         while(available_images_exp && (*iterExp).id()!=0){
 
             transMat_tr.initialize(myStreamTR);
             transMat_rt.initialize(myStreamRT);
             transMat_tr_mirror.initialize(myStreamTR);
             transMat_rt_mirror.initialize(myStreamRT);
 
             for(int i=0; i<available_images_proj; i++){
                 max_vector_tr[i]=-1;
                 max_vector_rt[i]=-1;
                 max_vector_tr_mirror[i]=-1;
                 max_vector_rt_mirror[i]=-1;
             }
 
             available_images_exp--;
 
             MDRow& rowExp = *iterExp;
             rowExp.getValue(MDL_IMAGE, fnImgExp);
             //std::cerr << expIndex << ". Image: " << fnImgExp << std::endl;
 
             //AJ calling the function to align the images
             bool mirror=false;
             //SFexp
             align_experimental_image(fnImgExp, d_referenceAux, d_experimentalAuxTR, d_experimentalAuxRT, transMat_tr, transMat_rt,
                     max_vector_tr, max_vector_rt, SF, available_images_proj, mirror, maxShift,
                     myhandlePadded_tr, myhandleMask_tr, myhandlePolar_tr, myhandlePaddedB_tr, myhandleMaskB_tr, myhandlePolarB_tr,
                     myhandlePadded_rt, myhandleMask_rt, myhandlePolar_rt, myhandlePaddedB_rt, myhandleMaskB_rt, myhandlePolarB_rt,
                     myStructureAux_tr, myStructureAux_rt, myStreamTR, myStreamRT,
                     resultTR, resultRT, original_image_stack, ifftcb);
 
 
             mirror=true;
             //SFexp
             align_experimental_image(fnImgExp, d_referenceAux, d_experimentalAuxTR, d_experimentalAuxRT, transMat_tr_mirror, transMat_rt_mirror,
                             max_vector_tr_mirror, max_vector_rt_mirror, SF, available_images_proj, mirror, maxShift,
                             myhandlePadded_tr, myhandleMask_tr, myhandlePolar_tr, myhandlePaddedB_tr, myhandleMaskB_tr, myhandlePolarB_tr,
                             myhandlePadded_rt, myhandleMask_rt, myhandlePolar_rt, myhandlePaddedB_rt, myhandleMaskB_rt, myhandlePolarB_rt,
                             myStructureAux_tr, myStructureAux_rt, myStreamTR, myStreamRT,
                             resultTR, resultRT, original_image_stack, ifftcb);
 
             //AJ to check the best transformation among all the evaluated
             transMat_tr.copyMatrixToCpu(myStreamTR);
             transMat_tr_mirror.copyMatrixToCpu(myStreamRT);
             transMat_rt.copyMatrixToCpu(myStreamTR);
             transMat_rt_mirror.copyMatrixToCpu(myStreamRT);
 
             waitGpu(myStreamTR, false);
             waitGpu(myStreamRT, false);
 
             MultidimArray<float> out2(3,3);
             for(int i=0; i<available_images_proj; i++){
                 if(max_vector_tr[i]>max_vector_rt[i]){
                     memcpy(MULTIDIM_ARRAY(out2), &transMat_tr.h_data[i*9], 9*sizeof(float));
                     matrixTransCpu[n].setSlice(firstIdx+i, out2);
                     A2D_ELEM(matrixCorrCpu, n, firstIdx+i) = max_vector_tr[i];
                 }else{
                     memcpy(MULTIDIM_ARRAY(out2), &transMat_rt.h_data[i*9], 9*sizeof(float));
                     matrixTransCpu[n].setSlice(firstIdx+i, out2);
                     A2D_ELEM(matrixCorrCpu, n, firstIdx+i) = max_vector_rt[i];
                 }
                 //mirror image
                 if(max_vector_tr_mirror[i]>max_vector_rt_mirror[i]){
                     memcpy(MULTIDIM_ARRAY(out2), &transMat_tr_mirror.h_data[i*9], 9*sizeof(float));
                     matrixTransCpu_mirror[n].setSlice(firstIdx+i, out2);
                     A2D_ELEM(matrixCorrCpu_mirror, n, firstIdx+i) = max_vector_tr_mirror[i];
                 }else{
                     memcpy(MULTIDIM_ARRAY(out2), &transMat_rt_mirror.h_data[i*9], 9*sizeof(float));
                     matrixTransCpu_mirror[n].setSlice(firstIdx+i, out2);
                     A2D_ELEM(matrixCorrCpu_mirror, n, firstIdx+i) = max_vector_rt_mirror[i];
                 }
             }
 
             if(iterExp != SF.end())
                 ++iterExp;
 
             n++;
             workDone+=available_images_proj;
             if (size_t(workDone/100)>lastProgressShown)
             {
                 progress_bar(workDone);
                 lastProgressShown=size_t(workDone/100);
             }
         }//end while experimental images
 
         firstIdx +=available_images_proj;
         int aux;
         aux=available_images_proj;
         if(firstIdx+available_images_proj > mdExpSize){ //mdInSize
             aux=available_images_proj;
             available_images_proj=mdExpSize-firstIdx; //mdInSize
         }
         if(firstIdx==mdExpSize){ //mdInSize
             finish=true;
         }
         if(aux!=available_images_proj){
             myhandlePadded_tr.clear();
             myhandleMask_tr.clear();
             myhandlePolar_tr.clear();
             myhandlePaddedB_tr.clear();
             myhandleMaskB_tr.clear();
             myhandlePolarB_tr.clear();
 
             myhandlePadded_rt.clear();
             myhandleMask_rt.clear();
             myhandlePolar_rt.clear();
             myhandlePaddedB_rt.clear();
             myhandleMaskB_rt.clear();
             myhandlePolarB_rt.clear();
         }
 
 
     }//End loop over the reference images while(!finish)
     progress_bar(totalWork);
 
     myhandlePadded_tr.clear();
     myhandleMask_tr.clear();
     myhandlePolar_tr.clear();
     myhandlePaddedB_tr.clear();
     myhandleMaskB_tr.clear();
     myhandlePolarB_tr.clear();
 
     myhandlePadded_rt.clear();
     myhandleMask_rt.clear();
     myhandlePolar_rt.clear();
     myhandlePaddedB_rt.clear();
     myhandleMaskB_rt.clear();
     myhandlePolarB_rt.clear();
 
     MultidimArray<float> weights(1,1,mdExpSize,2*mdInSize);
     MultidimArray<float> weightsMax;
     MultidimArray<float> corrTotalRow(1,1,mdExpSize, 2*mdInSize);
     int Nref;
     if(keepN){
         Nref=n_keep;
     }else if(significance){
         Nref=round(corrTotalRow.xdim*alpha);
         if(Nref==0)
             Nref=1;
     }
 
 
     calculate_weights(matrixCorrCpu, matrixCorrCpu_mirror, corrTotalRow, weights, Nref, mdExpSize, mdInSize, weightsMax, simplifiedMd,
             matrixTransCpu, matrixTransCpu_mirror, maxShift);
 
     std::cerr << "Creating output metadatas..." << std::endl;
 
     generate_metadata(SF, SFexp, fnDir, fn_out, mdExpSize, mdInSize, weights, corrTotalRow, matrixTransCpu,
             matrixTransCpu_mirror, maxShift, weightsMax, simplifiedMd, Nref);
 
     if(generate_out)
         generate_output_classes(SF, SFexp, fnDir, mdExpSize, mdInSize, weights, matrixTransCpu,
                 matrixTransCpu_mirror, maxShift, fn_classes_out, weightsMax, simplifiedMd, Nref);
 
     //Free memory in CPU
     for(int i=0; i<mdInSize; i++) //mdExpSize
         matrixTransCpu[i].coreDeallocate();
     delete []matrixTransCpu;
     for(int i=0; i<mdInSize; i++) //mdExpSize
         matrixTransCpu_mirror[i].coreDeallocate();
     delete []matrixTransCpu_mirror;
 
     cpuFree(max_vector_tr);
     cpuFree(max_vector_rt);
     cpuFree(max_vector_tr_mirror);
     cpuFree(max_vector_rt_mirror);
 
 
 
 }

◆ show()

void ProgGpuCorrelation::show ( )

Show.

Definition at line 526 of file xmipp_gpu_correlation.cpp.

 {
     std::cout
     << "Input projected:                " << fn_ref    << std::endl
     << "Input experimental:             " << fn_exp    << std::endl
     << "Generate output images (y/n):   " << generate_out    << std::endl
     ;
 }

Member Data Documentation

◆ SF

MetaDataVec ProgGpuCorrelation::SF

Definition at line 53 of file xmipp_gpu_correlation.h.

◆ SFexp

MetaDataVec ProgGpuCorrelation::SFexp

Definition at line 53 of file xmipp_gpu_correlation.h.

The documentation for this class was generated from the following files:

xmipp/legacy/libraries/reconstruction_adapt_cuda/xmipp_gpu_correlation.h
xmipp/legacy/libraries/reconstruction_adapt_cuda/xmipp_gpu_correlation.cpp

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ defineParams()

◆ readParams()

◆ run()

◆ show()

Member Data Documentation

◆ SF

◆ SFexp