ProgVolumeHomogenizer Class Reference

#include <volume_homogenizer.h>

Inheritance diagram for ProgVolumeHomogenizer:

Collaboration diagram for ProgVolumeHomogenizer:

Public Member Functions
	ProgVolumeHomogenizer ()
void	readParams ()
void	defineParams ()
void	xmipp2Opencv (const MultidimArray< double > &xmippArray, cv::Mat &opencvMat)
void	convert2Uint8 (cv::Mat opencvDoubleMat, cv::Mat &opencvUintMat)
void	opencv2Xmipp (const cv::Mat &opencvMat, MultidimArray< double > &xmippArray)
void	run ()
void	parameterList (size_t imgSize, int maxParamsTrial, float *paramList)
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	fnVol
FileName	fnRef
FileName	fnSetOfImgIn
FileName	fnSetOfImgOut
int	winSize
double	cutFreq
MetaDataDb	mdPartialParticles
size_t	rank
size_t	Nprocessors
double	pyr_scale
int	levels
int	iterations
bool	addToInput
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 39 of file volume_homogenizer.h.

Constructor & Destructor Documentation

ProgVolumeHomogenizer()

ProgVolumeHomogenizer::ProgVolumeHomogenizer

(

)

Definition at line 32 of file volume_homogenizer.cpp.

{
    rank = 0;
    Nprocessors = 1;
}

Member Function Documentation

convert2Uint8()

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

Definition at line 89 of file volume_homogenizer.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 ProgVolumeHomogenizer::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 52 of file volume_homogenizer.cpp.

{
    addUsageLine("Correcting the set of input images of the input volume with ");
    addUsageLine("respect to the reference map and using optical flow algorithm");
    addParamsLine("  -i <selfile>                   : Selfile with input volume. This is the volume that we aimed to deform its related particles.");
    addParamsLine("                                 : Then we will use its deformed (corrected) images to merge with reference map images.");
    addParamsLine("  -ref <selfile>                 : Selfile with reference volume. This is the volume that we use as the reference in OF algorithm.");
    addParamsLine("                                 : This is the volume that we will be deformed");
    addParamsLine("  -img <selfile>                 : Selfile with input aligned images. These images are related to the input map");
    addParamsLine(" [-o <rootname=\"\">]            : Output fileName");
    addParamsLine("                                 : Output contains the list of corrected (deformed) images and their angles with respect to the reference map");
    addParamsLine("                                 : NOTE: you do not need to define any extension for the output root name. ");
    addParamsLine("                                 : This program will automatically create one .stk and one .xmd");
    addParamsLine(" [--cutFreq <cutFreq=0.5>]       : cut-off frequency to use for low-pass filtering of input and reference volumes. This is digital frequency");
    addParamsLine(" [--pyr_scale <pyr_scale=0.5>]   : parameter, specifying the image scale (<1) to build pyramids for each image.");
    addParamsLine("                                 : pyr_scale=0.5 means a classical pyramid, where each next layer is twice smaller than the previous one.");
    addParamsLine(" [--levels <levels=2>]           : number of pyramid layers including the initial image; levels=1 means");
    addParamsLine("                                 : that no extra layers are created and only the original images are used.");
    addParamsLine(" [--winSize <winSize=50>]        : averaging window size; larger values increase the algorithm robustness to image noise and give more chances");
    addParamsLine("                                 : for fast motion detection, but yield more blurred motion field.");
    addParamsLine(" [--iterations <iterations=10>]  : number of iterations the algorithm does at each pyramid level.");
    addParamsLine(" [--auto]                        : Automatic determination of pyr_scale, levels and iterations");

    addExampleLine("xmipp_volume_homogenizer  -i volume1.vol -ref volume2.vol -img particlesOfVolume1.xmd -o outputParticles");
}

opencv2Xmipp()

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

Definition at line 97 of file volume_homogenizer.cpp.

{
    int h = opencvMat.rows;
    int w = opencvMat.cols;
    xmippArray.initZeros(h, w);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xmippArray)
    DIRECT_A2D_ELEM(xmippArray,i,j) = opencvMat.at<float>(i,j);
}

parameterList()

void ProgVolumeHomogenizer::parameterList	(	size_t	imgSize,
		int	maxParamsTrial,
		float *	paramList
	)

Definition at line 106 of file volume_homogenizer.cpp.

{
    srand (time(NULL));
    int offset = 0;
    for (int i=0; i < maxParamsTrial; i++)
    {
        if (i==0)
        {
            paramList[offset] = pyr_scale;
            paramList[1+offset] = levels;
            paramList[2+offset] = iterations;
        }
        else
        {
            paramList[offset] = (float(0.5*rand()) / RAND_MAX + 0.5); // fraction of image per pyramid 0.5, 2, winSize, 15
            paramList[1+offset] = (rand() % 10+1); // number of pyramids (int)
            paramList[2+offset] = (rand() % 15 + 1); // number of iterations  (int)
        }

        offset +=  numParams;

    }

#ifdef DEBUG
    std::cout << "\n" << std::endl;
    for (int i=0; i < maxParamsTrial*numParams; i++)
    {
        std::cout << paramList[i] << std::endl;
    }
#endif

}

readParams()

void ProgVolumeHomogenizer::readParams ( )

virtual

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 38 of file volume_homogenizer.cpp.

{
    fnVol           = getParam("-i");
    fnRef           = getParam("-ref");
    fnSetOfImgIn    = getParam("-img");
    fnSetOfImgOut   = getParam("-o");
    winSize         = getIntParam("--winSize");
    cutFreq         = getDoubleParam("--cutFreq");
    pyr_scale       = getDoubleParam("--pyr_scale");
    levels          = getIntParam("--levels");
    iterations      = getIntParam("--iterations");
    addToInput      = checkParam("--auto");
}

run()

void ProgVolumeHomogenizer::run ( )

virtual

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

Reimplemented from XmippProgram.

Definition at line 140 of file volume_homogenizer.cpp.

{
    FileName fn_proj, fnIn;
    Image<double> inV, refV;
    Image<double> imgIn, imgCorr;
    MetaDataDb setOfImgIn, setOfImgOut;
    double rot, tilt, psi;
    bool flip;
    Projection projIn, projRef;
    cv::Mat ProjIn, ProjRef, ProjIn8, ProjRef8, ImgIn;
    cv::Mat flow, ImgCorr;
    cv::Mat planes[]={flow, flow};
    FileName stackName = fnSetOfImgOut.removeAllExtensions() + ".stk";
    FileName mdName = fnSetOfImgOut.removeAllExtensions() + ".xmd";
    MDRowSql rowInput;

    inV.read(fnVol);
    inV().setXmippOrigin();

    refV.read(fnRef);
    refV().setXmippOrigin();

    setOfImgIn.read(fnSetOfImgIn);
    setOfImgIn.removeDisabled();

    fn_proj = stackName;
    createEmptyFile(fn_proj, XSIZE(inV()), YSIZE(inV()), 1, setOfImgIn.size(), true, WRITE_OVERWRITE);

    //filtering the input and reference volumes
    FourierFilter Filter;
    Filter.FilterBand=LOWPASS;
    Filter.FilterShape=RAISED_COSINE;
    Filter.w1=cutFreq;
    Filter.generateMask(inV());
    Filter.do_generate_3dmask=true;
    Filter.applyMaskSpace(inV());
    Filter.applyMaskSpace(refV());

    if (addToInput) //Do automatic OF parameter search
    {
        int maxImagsCheck = 5;
        int maxParamsTrial = 100;
        cv::Mat corr;
        double bestCorr = -1e6;
        float* paramList = new float[numParams*maxParamsTrial];

        std::cout << "Estimating optimal parameters for OF alignment" << std::endl;
        init_progress_bar(maxImagsCheck);

        size_t Xdim = 0;
        size_t Ydim = 0;
        size_t Zdim = 0;
        size_t Ndim = 0;

        ApplyGeoParams p;
        p.only_apply_shifts = true;
        imgIn.readApplyGeo(setOfImgIn, 1, p);
        imgIn.getDimensions(Xdim,Ydim,Zdim,Ndim);

        parameterList(Xdim, maxParamsTrial, paramList);

        for (size_t i = 0; i < maxImagsCheck; i++)
        {
            imgIn.readApplyGeo(setOfImgIn, i+1, p);

            setOfImgIn.getRow(rowInput,i+1);
            rowInput.getValue(MDL_ANGLE_ROT, rot);
            rowInput.getValue(MDL_ANGLE_TILT, tilt);
            rowInput.getValue(MDL_ANGLE_PSI, psi);
            if (rowInput.containsLabel(MDL_FLIP))
                rowInput.getValue(MDL_FLIP, flip);

            //Reprojection from input and reference volumes to calculate optical flow (OF)
            projectVolume(inV(), projIn, YSIZE(inV()), XSIZE(inV()), rot, tilt, psi);
            projectVolume(refV(), projRef, YSIZE(refV()), XSIZE(refV()), rot, tilt, psi);

            if (flip)
            {
                projIn.mirrorX();
                projRef.mirrorX();
            }

            //Preparing input data to use for OF algorithm
            xmipp2Opencv(projIn(), ImgIn); //We wamt to compare the projIn() with the projRef() imgIn here is only a holder to compare
            xmipp2Opencv(projIn(), ProjIn);
            xmipp2Opencv(projRef(), ProjRef);

            convert2Uint8(ProjIn,ProjIn8);
            convert2Uint8(ProjRef,ProjRef8);

            int offset = 0;
            for (size_t j = 0; j < maxParamsTrial; ++j)
            {
                cv::calcOpticalFlowFarneback(ProjRef8, ProjIn8, flow, paramList[0+offset], paramList[1+offset],paramList[2+offset], paramList[3+j], 7, 1.5, 0);
                cv::split(flow, planes);

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

                //applying the flow on the input image to use as the corrected one
                cv::remap(ImgIn, ImgCorr, planes[0], planes[1], cv::INTER_CUBIC);
                cv::matchTemplate(ImgCorr, ProjRef, corr, cv::TM_CCOEFF_NORMED);

                if( corr.at<float>(0,0) > bestCorr)
                {
                    bestCorr = corr.at<float>(0,0);

                    pyr_scale = paramList[0+offset];
                    levels = paramList[1+offset];
                    iterations = paramList[2+offset];
                }

                offset += numParams;

            }

            progress_bar(i+1);

        }


            std::cout << "Current parameters for OF: " << std::endl;
            std::cout << "fitness value : " << bestCorr << std::endl;
            std::cout << "Scale : " << pyr_scale << " " <<  " levels : "<< levels <<" WinSize : " << winSize << " Iterations : " << iterations << std::endl;

    }

    //calculating progress time
    size_t maxNImg = setOfImgIn.size();

    if (rank == 0)
        init_progress_bar(maxNImg);

    for (size_t i = 0; i < maxNImg; i++)
    {
        if ((i+1) % Nprocessors == rank)
        {
            ApplyGeoParams p;
            p.only_apply_shifts = true;
            imgIn.readApplyGeo(setOfImgIn, i+1, p);

            setOfImgIn.getRow(rowInput,i+1);
            rowInput.getValue(MDL_ANGLE_ROT, rot);
            rowInput.getValue(MDL_ANGLE_TILT, tilt);
            rowInput.getValue(MDL_ANGLE_PSI, psi);
            if (rowInput.containsLabel(MDL_FLIP))
                rowInput.getValue(MDL_FLIP, flip);

            //Reprojection from input and reference volumes to calculate optical flow (OF)
            projectVolume(inV(), projIn, YSIZE(inV()), XSIZE(inV()), rot, tilt, psi);
            projectVolume(refV(), projRef, YSIZE(refV()), XSIZE(refV()), rot, tilt, psi);

            if (flip)
            {
                projIn.mirrorX();
                projRef.mirrorX();
            }

            //preparing input data to use for OF algorithm
            xmipp2Opencv(imgIn(), ImgIn);
            xmipp2Opencv(projIn(), ProjIn);
            xmipp2Opencv(projRef(), ProjRef);
            convert2Uint8(ProjIn,ProjIn8);
            convert2Uint8(ProjRef,ProjRef8);

            /*OF algorithm/*
            Parameters:
            ----------
            prev         first 8-bit single-channel input image.
            next         second input image of the same size and the same type as prev.
            flow         computed flow image that has the same size as prev and type CV_32FC2.
            pyr_scale    parameter, specifying the image scale (<1) to build pyramids for each image; pyr_scale=0.5 means a classical pyramid, where each next layer is twice smaller than the previous one.
            levels       number of pyramid layers including the initial image; levels=1 means that no extra layers are created and only the original images are used.
            winsize      averaging window size; larger values increase the algorithm robustness to image noise and give more chances for fast motion detection, but yield more blurred motion field.
            iterations   number of iterations the algorithm does at each pyramid level.
            poly_n       size of the pixel neighborhood used to find polynomial expansion in each pixel; larger values mean that the image will be approximated with smoother surfaces, yielding more robust algorithm and more blurred motion field, typically poly_n =5 or 7.
            poly_sigma   standard deviation of the Gaussian that is used to smooth derivatives used as a basis for the polynomial expansion; for poly_n=5, you can set poly_sigma=1.1, for poly_n=7, a good value would be poly_sigma=1.5.
            flags
             */

            //cv::calcOpticalFlowFarneback(ProjRef8, ProjIn8, flow, 0.9, 5, winSize, 15, 7, 1.5, 0);
            //Moshen: cv::calcOpticalFlowFarneback(ProjRef8, ProjIn8, flow, 0.5, 6, winSize, 3, 5, 1.1, 0);
            //JV: cv::calcOpticalFlowFarneback(ProjRef8, ProjIn8, flow, 0.5, 2, winSize, 15, 7, 1.5, 0);


            cv::calcOpticalFlowFarneback(ProjRef8, ProjIn8, flow, pyr_scale, levels, winSize, iterations, 7, 1.5, 0);
            cv::split(flow, planes);

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

            //applying the flow on the input image to use as the corrected one
            cv::remap(ImgIn, ImgCorr, planes[0], planes[1], cv::INTER_CUBIC);

            //preparing output data to use for xmipp
            opencv2Xmipp(ImgCorr, imgCorr());

#ifdef DEBUG
            if (i==5)
            {
                imgIn.write("test_raw.spi"); // Raw image to deform Not filtered or modifed
                projIn.write("test_projIn.spi"); //map projection of conformation to deform (Noise free and filter)
                projRef.write("test_projRef.spi"); //map projection of conformation used as reference (deform to this conformation) (Noise free and filter)

                cv::remap(ProjIn, ImgCorr, planes[0], planes[1], cv::INTER_CUBIC); //apply OF to the map projetion to see something
                opencv2Xmipp(ImgCorr, imgCorr()); //transform to Xmipp
                imgCorr.write("test_ProjIn_Deformed.spi"); //  Image deformed

                Image<double> imgFlowX, imgFlowY;
                opencv2Xmipp(planes[0], imgFlowX());
                opencv2Xmipp(planes[1], imgFlowY());
                imgFlowX.write("test_FlowX.spi");
                imgFlowY.write("test_FlowY.spi");
            }
#endif

            //filling output metaData
            fn_proj.compose(i+1, stackName);

            rowInput.setValue(MDL_IMAGE, fn_proj);
            rowInput.setValue(MDL_SHIFT_X, 0.0);
            rowInput.setValue(MDL_SHIFT_Y, 0.0);
            rowInput.setValue(MDL_SHIFT_Z, 0.0);

            imgCorr.write(fn_proj, i+1,true, WRITE_OVERWRITE);

            mdPartialParticles.addRow(rowInput);
            rowInput.clear();
        }

        if (rank==0)
            progress_bar(i+1);
    }

    synchronize();
    gatherResults();

    if (rank == 0)
        mdPartialParticles.write(mdName);

}

xmipp2Opencv()

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

Definition at line 79 of file volume_homogenizer.cpp.

{
    int h = YSIZE(xmippArray);
    int w = XSIZE(xmippArray);
    opencvMat = cv::Mat::zeros(h, w,CV_32FC1);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xmippArray)
    opencvMat.at<float>(i,j) = (float) DIRECT_A2D_ELEM(xmippArray,i,j);
}

Member Data Documentation

addToInput

bool ProgVolumeHomogenizer::addToInput

Definition at line 58 of file volume_homogenizer.h.

cutFreq

double ProgVolumeHomogenizer::cutFreq

Definition at line 46 of file volume_homogenizer.h.

fnRef

FileName ProgVolumeHomogenizer::fnRef

Definition at line 42 of file volume_homogenizer.h.

fnSetOfImgIn

FileName ProgVolumeHomogenizer::fnSetOfImgIn

Definition at line 42 of file volume_homogenizer.h.

fnSetOfImgOut

FileName ProgVolumeHomogenizer::fnSetOfImgOut

Definition at line 42 of file volume_homogenizer.h.

fnVol

FileName ProgVolumeHomogenizer::fnVol

Definition at line 42 of file volume_homogenizer.h.

iterations

int ProgVolumeHomogenizer::iterations

Definition at line 56 of file volume_homogenizer.h.

levels

int ProgVolumeHomogenizer::levels

Definition at line 54 of file volume_homogenizer.h.

mdPartialParticles

MetaDataDb ProgVolumeHomogenizer::mdPartialParticles

Definition at line 48 of file volume_homogenizer.h.

Nprocessors

size_t ProgVolumeHomogenizer::Nprocessors

Definition at line 50 of file volume_homogenizer.h.

pyr_scale

double ProgVolumeHomogenizer::pyr_scale

Definition at line 52 of file volume_homogenizer.h.

rank

size_t ProgVolumeHomogenizer::rank

Definition at line 50 of file volume_homogenizer.h.

winSize

int ProgVolumeHomogenizer::winSize

Definition at line 44 of file volume_homogenizer.h.

---

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

• xmipp/applications/programs/mpi_volume_homogenizer/volume_homogenizer.h
• xmipp/applications/programs/mpi_volume_homogenizer/volume_homogenizer.cpp