ProgVolDeformSph Class Reference

#include <volume_deform_sph.h>

Inheritance diagram for ProgVolDeformSph:

Collaboration diagram for ProgVolDeformSph:

Public Member Functions
void	defineParams ()
	Define params. More...
void	readParams ()
	Read arguments from command line. More...
void	show ()
	Show. More...
double	distance (double *pclnm)
	Distance. More...
void	run ()
	Run. More...
void	minimizepos (int l2, Matrix1D< double > &steps) const
	Determine the positions to be minimize of a vector containing spherical harmonic coefficients. More...
void	numCoefficients (int l1, int l2, int &nc) const
	Length of coefficients vector. More...
void	fillVectorTerms (int l1, int l2)
	Zernike and SPH coefficients allocation. More...
void	computeStrain ()
	Compute strain. More...
void	writeVector (std::string const &outPath, Matrix1D< double > const &v, bool append) const
	Save vector to file. More...
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	fnVolI
	Volume to deform. More...
FileName	fnVolR
	Reference volume. More...
FileName	fnVolOut
	Output Volume (deformed input volume) More...
FileName	fnRoot
	Root name for several output files. More...
bool	analyzeStrain
	Save the deformation of each voxel for local strain and rotation analysis. More...
bool	optimizeRadius
	Radius optimization. More...
int	L1
	Degree of Zernike polynomials and spherical harmonics. More...
int	L2
std::vector< double >	sigma
	Gaussian width to filter the volumes. More...
std::vector< Image< double > >	volumesI
	Image Vector. More...
std::vector< Image< double > >	volumesR
double	Rmax
	Maximum radius for the transformation. More...
int	vecSize
	Coefficient vector size. More...
Image< double >	VI
	Images. More...
Image< double >	VR
Image< double >	VO
Image< double >	Gx
Image< double >	Gy
Image< double >	Gz
std::vector< double >	absMaxR_vec
	Maxima of reference volumes (in absolute value) More...
double	deformation
double	sumVI
double	sumVD
double	lambda
bool	applyTransformation
bool	saveDeformation
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

Sph Alignment Parameters.

Definition at line 38 of file volume_deform_sph.h.

Member Function Documentation

computeStrain()

void ProgVolDeformSph::computeStrain

(

)

Compute strain.

Definition at line 514 of file volume_deform_sph.cpp.

{
    Image<double> LS, LR;
    LS().initZeros(Gx());
    LR().initZeros(Gx());

    // Gaussian filter of the derivatives
    FourierFilter f;
    f.FilterBand=LOWPASS;
    f.FilterShape=REALGAUSSIAN;
    f.w1=2;
    f.applyMaskSpace(Gx());
    f.applyMaskSpace(Gy());
    f.applyMaskSpace(Gz());

    Gx.write(fnRoot+"_PPPGx.vol");
    Gy.write(fnRoot+"_PPPGy.vol");
    Gz.write(fnRoot+"_PPPGz.vol");

    MultidimArray<double> &mLS=LS();
    MultidimArray<double> &mLR=LR();
    MultidimArray<double> &mGx=Gx();
    MultidimArray<double> &mGy=Gy();
    MultidimArray<double> &mGz=Gz();
    Matrix2D<double> U(3,3), D(3,3), H(3,3);
    std::vector< std::complex<double> > eigs;
    H.initZeros();
    for (int k=STARTINGZ(mLS)+2; k<=FINISHINGZ(mLS)-2; ++k)
    {
        int km1=k-1;
        int kp1=k+1;
        int km2=k-2;
        int kp2=k+2;
        for (int i=STARTINGY(mLS)+2; i<=FINISHINGY(mLS)-2; ++i)
        {
            int im1=i-1;
            int ip1=i+1;
            int im2=i-2;
            int ip2=i+2;
            for (int j=STARTINGX(mLS)+2; j<=FINISHINGX(mLS)-2; ++j)
            {
                int jm1=j-1;
                int jp1=j+1;
                int jm2=j-2;
                int jp2=j+2;
                MAT_ELEM(U,0,0)=Dx(mGx); MAT_ELEM(U,0,1)=Dy(mGx); MAT_ELEM(U,0,2)=Dz(mGx);
                MAT_ELEM(U,1,0)=Dx(mGy); MAT_ELEM(U,1,1)=Dy(mGy); MAT_ELEM(U,1,2)=Dz(mGy);
                MAT_ELEM(U,2,0)=Dx(mGz); MAT_ELEM(U,2,1)=Dy(mGz); MAT_ELEM(U,2,2)=Dz(mGz);

                MAT_ELEM(D,0,0) = MAT_ELEM(U,0,0);
                MAT_ELEM(D,0,1) = MAT_ELEM(D,1,0) = 0.5*(MAT_ELEM(U,0,1)+MAT_ELEM(U,1,0));
                MAT_ELEM(D,0,2) = MAT_ELEM(D,2,0) = 0.5*(MAT_ELEM(U,0,2)+MAT_ELEM(U,2,0));
                MAT_ELEM(D,1,1) = MAT_ELEM(U,1,1);
                MAT_ELEM(D,1,2) = MAT_ELEM(D,2,1) = 0.5*(MAT_ELEM(U,1,2)+MAT_ELEM(U,2,1));
                MAT_ELEM(D,2,2) = MAT_ELEM(U,2,2);

                MAT_ELEM(H,0,1) = 0.5*(MAT_ELEM(U,0,1)-MAT_ELEM(U,1,0));
                MAT_ELEM(H,0,2) = 0.5*(MAT_ELEM(U,0,2)-MAT_ELEM(U,2,0));
                MAT_ELEM(H,1,2) = 0.5*(MAT_ELEM(U,1,2)-MAT_ELEM(U,2,1));
                MAT_ELEM(H,1,0) = -MAT_ELEM(H,0,1);
                MAT_ELEM(H,2,0) = -MAT_ELEM(H,0,2);
                MAT_ELEM(H,2,1) = -MAT_ELEM(H,1,2);

                A3D_ELEM(mLS,k,i,j)=fabs(D.det());
                allEigs(H,eigs);
                for (size_t n=0; n < eigs.size(); n++)
                {
                    double imagabs=fabs(eigs[n].imag());
                    if (imagabs>1e-6)
                    {
                        A3D_ELEM(mLR,k,i,j)=imagabs*180/PI;
                        break;
                    }
                }
            }
        }
        LS.write(fnVolOut.withoutExtension()+"_strain.mrc");
        LR.write(fnVolOut.withoutExtension()+"_rotation.mrc");
    }
}

defineParams()

void ProgVolDeformSph::defineParams ( )

virtual

Define params.

Reimplemented from XmippProgram.

Definition at line 35 of file volume_deform_sph.cpp.

                                    {
    addUsageLine("Compute the deformation that properly fits two volumes using spherical harmonics");
    addParamsLine("   -i <volume>                         : Volume to deform");
    addParamsLine("   -r <volume>                         : Reference volume");
    addParamsLine("  [-o <volume=\"\">]                   : Output volume which is the deformed input volume");
    addParamsLine("  [--oroot <rootname=\"Volumes\">]     : Root name for output files");
    addParamsLine("                                       : By default, the input file is rewritten");
    addParamsLine("  [--sigma <Matrix1D=\"\">]            : Sigma values to filter the volume to perform a multiresolution analysis");
    addParamsLine("  [--analyzeStrain]                    : Save the deformation of each voxel for local strain and rotation analysis");
    addParamsLine("  [--optimizeRadius]                   : Optimize the radius of each spherical harmonic");
    addParamsLine("  [--l1 <l1=3>]                        : Degree Zernike Polynomials=1,2,3,...");
    addParamsLine("  [--l2 <l2=2>]                        : Harmonical depth of the deformation=1,2,3,...");
    addParamsLine("  [--regularization <l=0.00025>]       : Regularization weight");
    addParamsLine("  [--Rmax <r=-1>]                      : Maximum radius for the transformation");
    addParamsLine("   [--thr <N=-1>]                      : Maximal number of the processing CPU threads");
    addExampleLine("xmipp_volume_deform_sph -i vol1.vol -r vol2.vol -o vol1DeformedTo2.vol");
}

distance()

double ProgVolDeformSph::distance

(

double *

pclnm

)

Distance.

Definition at line 256 of file volume_deform_sph.cpp.

{
    if (applyTransformation)
    {
        VO().initZeros(VR());
        VO().setXmippOrigin();
    }
    const MultidimArray<double> &mVR=VR();
    const MultidimArray<double> &mVI=VI();
    const size_t size = m_clnm.size();
    for (size_t i = 0; i < size; ++i) {
        auto &p = m_clnm[i];
        p.x = pclnm[i + 1];
        p.y = pclnm[i + size + 1];
        p.z = pclnm[i + size + size + 1];
    }

    const auto distance_vals = computeDistance();
    deformation=std::sqrt(distance_vals.modg / (double)distance_vals.count);
    sumVD = distance_vals.VD;

    if (applyTransformation)
        VO.write(fnVolOut);

    double massDiff=std::abs(sumVI-sumVD)/sumVI;
    return std::sqrt(distance_vals.diff / (double)distance_vals.count)+lambda*(deformation+massDiff);
}

fillVectorTerms()

void ProgVolDeformSph::fillVectorTerms	(	int	l1,
		int	l2
	)

Zernike and SPH coefficients allocation.

Definition at line 486 of file volume_deform_sph.cpp.

{
    m_zshVals.resize(vecSize);
    int idx = 0;
    for (int h=0; h<=l2; h++)
    {
        int totalSPH = 2*h+1;
        auto aux = (int)(std::floor(totalSPH/2));
        for (int l=h; l<=l1; l+=2)
        {
            for (int m=0; m<totalSPH; m++)
            {
                auto &t = m_zshVals[idx];
                t.l1 = l;
                t.n = h;
                t.l2 = h;
                t.m = m - aux;
                idx++;
            }
        }
    }
}

minimizepos()

void ProgVolDeformSph::minimizepos	(	int	l2,
		Matrix1D< double > &	steps
	)		const

Determine the positions to be minimize of a vector containing spherical harmonic coefficients.

Definition at line 454 of file volume_deform_sph.cpp.

{
    int size = 0;
    numCoefficients(L1,l2,size);
    auto totalSize = (int)(steps.size()/3);
    for (int idx=0; idx<size; idx++)
    {
        VEC_ELEM(steps,idx) = 1;
        VEC_ELEM(steps,idx+totalSize) = 1;
        VEC_ELEM(steps,idx+2*totalSize) = 1;
    }   
}

numCoefficients()

void ProgVolDeformSph::numCoefficients	(	int	l1,
		int	l2,
		int &	nc
	)		const

Length of coefficients vector.

Definition at line 467 of file volume_deform_sph.cpp.

{
    // l1 -> Degree Zernike
    // l2 & h --> Degree SPH
    for (int h=0; h<=l2; h++)
    {
        // For the current SPH degree (h), determine the number of SPH components/equations
        int numSPH = 2*h+1;
        // Finf the total number of radial components with even degree for a given l1 and h
        int count=l1-h+1;
        int numEven=(count>>1)+(count&1 && !(h&1));
        // Total number of components is the number of SPH as many times as Zernike components
        if (h%2 == 0)
            nc += numSPH*numEven;
        else
            nc += numSPH*(l1-h+1-numEven);
    }
}

readParams()

void ProgVolDeformSph::readParams ( )

virtual

Read arguments from command line.

Reimplemented from XmippProgram.

Definition at line 54 of file volume_deform_sph.cpp.

                                  {
    std::string aux;
    fnVolI = getParam("-i");
    fnVolR = getParam("-r");
    L1 = getIntParam("--l1");
    L2 = getIntParam("--l2");
    fnRoot = getParam("--oroot");
    
    aux = getParam("--sigma");
    // Transform string ov values separated by white spaces into substrings stored in a vector
    std::stringstream ss(aux);
    std::istream_iterator<std::string> begin(ss);
    std::istream_iterator<std::string> end;
    std::vector<std::string> vstrings(begin, end);
    sigma.resize(vstrings.size());
    std::transform(vstrings.begin(), vstrings.end(), sigma.begin(), [](const std::string& val)
    {
        return std::stod(val);
    });

    fnVolOut = getParam("-o");
    if (fnVolOut=="")
        fnVolOut=fnVolI;
    analyzeStrain=checkParam("--analyzeStrain");
    optimizeRadius=checkParam("--optimizeRadius");
    lambda = getDoubleParam("--regularization");
    Rmax = getDoubleParam("--Rmax");
    applyTransformation = false;
    int threads = getIntParam("--thr");
    if (0 >= threads) {
        threads = CPU::findCores();
    }
    m_threadPool.resize(threads);
}

run()

void ProgVolDeformSph::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 292 of file volume_deform_sph.cpp.

                           {
    saveDeformation=false;

    VI.read(fnVolI);
    VR.read(fnVolR);
    sumVI = 0.0;
    if (Rmax<0)
        Rmax=XSIZE(VI())/2;

    VI().setXmippOrigin();
    VR().setXmippOrigin();

    // Filter input and reference volumes according to the values of sigma
    FourierFilter filter;
    filter.FilterShape = REALGAUSSIAN;
    filter.FilterBand = LOWPASS;
    filter.generateMask(VI());

    // We need also to normalized the filtered volumes to compare them appropiately
    Image<double> auxI = VI;
    Image<double> auxR = VR;

    MultidimArray<int> bg_mask;
    bg_mask.resizeNoCopy(VI().zdim, VI().ydim, VI().xdim);
    bg_mask.setXmippOrigin();
    normalize_Robust(auxI(), bg_mask, true);
    bg_mask *= 0;
    normalize_Robust(auxR(), bg_mask, true);

    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(auxI())
    {
        if (DIRECT_A3D_ELEM(auxI(),k,i,j) >= 0.0)
            sumVI += DIRECT_A3D_ELEM(auxI(),k,i,j);
    }

    volumesI.push_back(auxI());
    volumesR.push_back(auxR());
    if (sigma.size() > 1 || sigma[0] != 0)
    {
        for (int ids=0; ids<sigma.size(); ids++)
        {
            Image<double> auxI = VI;
            Image<double> auxR = VR;
            filter.w1 = sigma[ids];

            // Filer input vol
            filter.do_generate_3dmask = true;
            filter.applyMaskSpace(auxI());
            bg_mask *= 0;
            normalize_Robust(auxI(), bg_mask, true);
            volumesI.push_back(auxI);
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(auxI())
            {
                if (DIRECT_A3D_ELEM(auxI(),k,i,j) >= 0.0)
                    sumVI += DIRECT_A3D_ELEM(auxI(),k,i,j);
            }

            // Filter ref vol
            filter.applyMaskSpace(auxR());
            bg_mask *= 0;
            normalize_Robust(auxR(), bg_mask, true);
            volumesR.push_back(auxR);
        }
    }

    Matrix1D<double> steps;
    Matrix1D<double> x;
    vecSize = 0;
    numCoefficients(L1,L2,vecSize);
    size_t totalSize = 3*vecSize;
    fillVectorTerms(L1,L2);
    m_clnm.resize(vecSize);
    x.initZeros(totalSize);
    for (int h=0;h<=L2;h++)
    {
        steps.clear();
        steps.initZeros(totalSize);
        minimizepos(h,steps);

        std::cout<<std::endl;
        std::cout<<"-------------------------- Basis Degrees: ("<<L1<<","<<h<<") --------------------------"<<std::endl;
#ifdef NEVERDEFINED
        for (int d=0;d<VEC_XSIZE(x);d++)
        {
            if (x(d)==0)
            {
                steps(d) = 1;
            }
            else if (d>=VEC_XSIZE(steps)+nh(h)-nh(h+1)&d<VEC_XSIZE(steps)+nh(h+1))
            {
                steps(d) = 1;
            }
            else
            {
                steps(d) = 0;
            }
            //std::cout<<steps(d)<<" ";
        }
        //std::cout<<std::endl;
#endif
        int iter;
        double fitness;
        powellOptimizer(x, 1, (int)(totalSize), &volDeformSphGoal, this,
                        0.01, fitness, iter, steps, true);

        std::cout<<std::endl;
        std::cout << "Deformation " << deformation << std::endl;
        std::ofstream deformFile;
        deformFile.open (fnRoot+"_deformation.txt");
        deformFile << deformation;
        deformFile.close();

// #define DEBUG
#ifdef DEBUG
    Image<double> save;
    save() = VI();
    save.write(fnRoot+"_PPPIdeformed.vol");
    save()-=VR();
    save.write(fnRoot+"_PPPdiff.vol");
    save()=VR();
    save.write(fnRoot+"_PPPR.vol");
    std::cout << "Error=" << deformation << std::endl;
    std::cout << "Press any key\n";
    char c; std::cin >> c;
#endif

    }
    applyTransformation=true;
    Matrix1D<double> degrees;
    degrees.initZeros(3);
    VEC_ELEM(degrees,0) = L1;
    VEC_ELEM(degrees,1) = L2;
    VEC_ELEM(degrees,2) = Rmax;
    writeVector(fnRoot+"_clnm.txt", degrees, false);
    writeVector(fnRoot+"_clnm.txt", x, true);
    if (analyzeStrain)
    {
        saveDeformation=true;
        Gx().initZeros(VR());
        Gy().initZeros(VR());
        Gz().initZeros(VR());
        Gx().setXmippOrigin();
        Gy().setXmippOrigin();
        Gz().setXmippOrigin();
    }

    distance(x.adaptForNumericalRecipes()); // To save the output volume

#ifdef DEBUG
        Image<double> save;
        save() = Gx();
        save.write(fnRoot+"_PPPGx.vol");
        save() = Gy();
        save.write(fnRoot+"_PPPGy.vol");
        save() = Gz();
        save.write(fnRoot+"_PPPGz.vol");
#endif

    if (analyzeStrain)
        computeStrain();
}

show()

void ProgVolDeformSph::show

(

)

Show.

Definition at line 90 of file volume_deform_sph.cpp.

                            {
    if (verbose==0)
        return;
    std::cout
            << "Volume to deform:     " << fnVolI         << std::endl
            << "Reference volume:     " << fnVolR         << std::endl
            << "Output volume:        " << fnVolOut       << std::endl
            << "Zernike Degree:       " << L1             << std::endl
            << "SH Degree:            " << L2             << std::endl
            << "Save deformation:     " << analyzeStrain  << std::endl
            << "Regularization:       " << lambda         << std::endl
    ;

}

writeVector()

void ProgVolDeformSph::writeVector	(	std::string const &	outPath,
		Matrix1D< double > const &	v,
		bool	append
	)		const

Save vector to file.

Definition at line 595 of file volume_deform_sph.cpp.

{
    std::ofstream outFile;
    if (append)
        outFile.open(outPath, std::ios_base::app);
    else
        outFile.open(outPath);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
        outFile << VEC_ELEM(v,i) << " ";
    outFile << std::endl;
}

Member Data Documentation

absMaxR_vec

std::vector<double> ProgVolDeformSph::absMaxR_vec

Maxima of reference volumes (in absolute value)

Definition at line 86 of file volume_deform_sph.h.

analyzeStrain

bool ProgVolDeformSph::analyzeStrain

Save the deformation of each voxel for local strain and rotation analysis.

Definition at line 54 of file volume_deform_sph.h.

applyTransformation

bool ProgVolDeformSph::applyTransformation

Definition at line 97 of file volume_deform_sph.h.

deformation

double ProgVolDeformSph::deformation

Definition at line 89 of file volume_deform_sph.h.

fnRoot

FileName ProgVolDeformSph::fnRoot

Root name for several output files.

Definition at line 51 of file volume_deform_sph.h.

fnVolI

FileName ProgVolDeformSph::fnVolI

Volume to deform.

Definition at line 42 of file volume_deform_sph.h.

fnVolOut

FileName ProgVolDeformSph::fnVolOut

Output Volume (deformed input volume)

Definition at line 48 of file volume_deform_sph.h.

fnVolR

FileName ProgVolDeformSph::fnVolR

Reference volume.

Definition at line 45 of file volume_deform_sph.h.

Gx

Image<double> ProgVolDeformSph::Gx

Definition at line 81 of file volume_deform_sph.h.

Gy

Image<double> ProgVolDeformSph::Gy

Definition at line 82 of file volume_deform_sph.h.

Gz

Image<double> ProgVolDeformSph::Gz

Definition at line 83 of file volume_deform_sph.h.

L1

int ProgVolDeformSph::L1

Degree of Zernike polynomials and spherical harmonics.

Definition at line 60 of file volume_deform_sph.h.

L2

int ProgVolDeformSph::L2

Definition at line 61 of file volume_deform_sph.h.

lambda

double ProgVolDeformSph::lambda

Definition at line 94 of file volume_deform_sph.h.

optimizeRadius

bool ProgVolDeformSph::optimizeRadius

Radius optimization.

Definition at line 57 of file volume_deform_sph.h.

Rmax

double ProgVolDeformSph::Rmax

Maximum radius for the transformation.

Definition at line 71 of file volume_deform_sph.h.

saveDeformation

bool ProgVolDeformSph::saveDeformation

Definition at line 100 of file volume_deform_sph.h.

sigma

std::vector<double> ProgVolDeformSph::sigma

Gaussian width to filter the volumes.

Definition at line 64 of file volume_deform_sph.h.

sumVD

double ProgVolDeformSph::sumVD

Definition at line 91 of file volume_deform_sph.h.

sumVI

double ProgVolDeformSph::sumVI

Definition at line 90 of file volume_deform_sph.h.

vecSize

int ProgVolDeformSph::vecSize

Coefficient vector size.

Definition at line 75 of file volume_deform_sph.h.

VI

Image<double> ProgVolDeformSph::VI

Images.

Definition at line 78 of file volume_deform_sph.h.

VO

Image<double> ProgVolDeformSph::VO

Definition at line 80 of file volume_deform_sph.h.

volumesI

std::vector<Image<double> > ProgVolDeformSph::volumesI

Image Vector.

Definition at line 67 of file volume_deform_sph.h.

volumesR

std::vector<Image<double> > ProgVolDeformSph::volumesR

Definition at line 68 of file volume_deform_sph.h.

VR

Image<double> ProgVolDeformSph::VR

Definition at line 79 of file volume_deform_sph.h.

---

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

• xmipp/libraries/reconstruction/volume_deform_sph.h
• xmipp/libraries/reconstruction/volume_deform_sph.cpp