ProgReconsADMM Class Reference

#include <reconstruct_ADMM.h>

Inheritance diagram for ProgReconsADMM:

Collaboration diagram for ProgReconsADMM:

Public Member Functions
	ProgReconsADMM ()
void	defineParams ()
void	readParams ()
void	show ()
void	run ()
void	produceSideInfo ()
void	project (double rot, double tilt, double psi, MultidimArray< double > &P, bool adjoint=false, double weight=1.)
void	project (const Matrix1D< double > &r1, const Matrix1D< double > &r2, MultidimArray< double > &P, bool adjoint=false, double weight=1.)
void	constructHtb ()
void	computeHtKH (MultidimArray< double > &kernelV)
void	addRegularizationTerms ()
void	applyKernel3D (MultidimArray< double > &x, MultidimArray< double > &AtAx)
void	applyLFilter (MultidimArray< std::complex< double > > &fourierL, bool adjoint=false)
void	applyLtFilter (MultidimArray< std::complex< double > > &fourierL, MultidimArray< double > &u, MultidimArray< double > &d)
void	applyConjugateGradient ()
void	doPOCSProjection ()
void	updateUD ()
void	produceVolume ()
virtual void	shareVolume (MultidimArray< double > &V)
virtual void	synchronize ()
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	fnIn
FileName	fnRoot
FileName	fnFirst
FileName	fnHtb
FileName	fnHtKH
String	kernelShape
double	a
double	alpha
double	mu
double	lambda
double	lambda1
double	Ti
double	Tp
bool	useWeights
bool	useCTF
int	Ncgiter
int	Nadmmiter
bool	positivity
bool	applyMask
double	Ts
Mask	mask
String	symmetry
bool	saveIntermediate
size_t	Nprocs
size_t	rank
AdmmKernel	kernel
Image< double >	CHtb
Image< double >	Ck
Image< double >	Vk
MetaDataVec	mdIn
MultidimArray< std::complex< double > >	fourierKernelV
MultidimArray< double >	paddedx
FourierTransformer	transformerPaddedx
FourierTransformer	transformerL
MultidimArray< double >	ux
MultidimArray< double >	uy
MultidimArray< double >	uz
MultidimArray< double >	dx
MultidimArray< double >	dy
MultidimArray< double >	dz
MultidimArray< double >	ud
MultidimArray< std::complex< double > >	fourierLx
MultidimArray< std::complex< double > >	fourierLy
MultidimArray< std::complex< double > >	fourierLz
SymList	SL
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 69 of file reconstruct_ADMM.h.

Constructor & Destructor Documentation

ProgReconsADMM()

ProgReconsADMM::ProgReconsADMM

(

)

Definition at line 39 of file reconstruct_ADMM.cpp.

                              : XmippProgram()
{
    rank=0;
    Nprocs=1;
}

Member Function Documentation

addRegularizationTerms()

void ProgReconsADMM::addRegularizationTerms

(

)

Add regularization to the kernel

Definition at line 467 of file reconstruct_ADMM.cpp.

{
    MultidimArray<double> L;
    L.initZeros(2*ZSIZE(CHtb())-1,2*YSIZE(CHtb())-1,2*XSIZE(CHtb())-1);
    L.setXmippOrigin();
    FourierTransformer transformer;
    transformer.setReal(L);

    if (mu>0)
    {
        addGradientTerm(mu, kernel,L,transformer,fourierKernelV, 'x');
        addGradientTerm(mu, kernel,L,transformer,fourierKernelV, 'y');
        addGradientTerm(mu, kernel,L,transformer,fourierKernelV, 'z');
    }

    if (lambda1>0)
    {
        L.initZeros();
        L(0,0,0)=lambda1;
        transformer.FourierTransform();
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(fourierKernelV)
            DIRECT_MULTIDIM_ELEM(fourierKernelV,n)+=DIRECT_MULTIDIM_ELEM(transformer.fFourier,n);
    }
}

applyConjugateGradient()

void ProgReconsADMM::applyConjugateGradient

(

)

Apply conjugate gradient

Definition at line 543 of file reconstruct_ADMM.cpp.

{
    // Conjugate gradient is solving Ax=b, when A is symmetric (=> positive semidefinite) and we apply it
    // to the problem
    // (H^T K H+mu L^TL + lambda I) x = H^Tb + mu L^T(u-d)
    // Being H the projection operator
    //       K the CTF operator
    //       L the gradient operator

    // Compute H^tb+mu*L^t(u-d)
    MultidimArray<double> r;
    applyLtFilter(fourierLx,ux,dx);
    r=ud;
    applyLtFilter(fourierLy,uy,dy);
    r+=ud;
    applyLtFilter(fourierLz,uz,dz);
    r+=ud;
    r*=mu;

    r+=CHtb();

    // Apply A^tA to the current estimate of the reconstruction
    MultidimArray<double> AtAVk;
    applyKernel3D(Ck(),AtAVk);

    // Compute first residual. This is the negative gradient of ||Ax-b||^2
    r-=AtAVk;
    double d=r.sum2();

    // Search direction
    MultidimArray<double> p, AtAp;
    p=r;

    // Perform CG iterations
    MultidimArray<double> &mCk=Ck();

    for (int iter=0; iter<Ncgiter; ++iter)
    {
        applyKernel3D(p,AtAp);
        double alpha=d/p.dotProduct(AtAp);

        // Update residual and current estimate
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(r)
        {
            DIRECT_MULTIDIM_ELEM(r,n)  -=alpha*DIRECT_MULTIDIM_ELEM(AtAp,n);
            DIRECT_MULTIDIM_ELEM(mCk,n)+=alpha*DIRECT_MULTIDIM_ELEM(p,n);
        }
        double newd=r.sum2();
        double beta=newd/d;
        d=newd;

        // Update search direction
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(p)
            DIRECT_MULTIDIM_ELEM(p,n)  = DIRECT_MULTIDIM_ELEM(r,n)+beta*DIRECT_MULTIDIM_ELEM(p,n);

//      Image<double> save;
//      save()=Vk();
//      save.write("PPPVk.vol");
//      std::cout << "Press any key" << std::endl;
//      char c; std::cin >> c;
    }
}

applyKernel3D()

void ProgReconsADMM::applyKernel3D	(	MultidimArray< double > &	x,
		MultidimArray< double > &	AtAx
	)

Apply kernel 3D

Definition at line 492 of file reconstruct_ADMM.cpp.

{
    paddedx.initZeros(2*ZSIZE(x)-1,2*YSIZE(x)-1,2*XSIZE(x)-1);
    paddedx.setXmippOrigin();

    // Copy Vk into paddedx
    for (int k=STARTINGZ(x); k<=FINISHINGZ(x); ++k)
        for (int i=STARTINGY(x); i<=FINISHINGY(x); ++i)
            memcpy(&A3D_ELEM(paddedx,k,i,STARTINGX(x)),&A3D_ELEM(x,k,i,STARTINGX(x)),XSIZE(x)*sizeof(double));

    // Compute Fourier transform of paddedx
    transformerPaddedx.setReal(paddedx);
    transformerPaddedx.FourierTransform();

    // Apply kernel
    double K=MULTIDIM_SIZE(paddedx);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(transformerPaddedx.fFourier)
    {
        DIRECT_MULTIDIM_ELEM(transformerPaddedx.fFourier,n)*=DIRECT_MULTIDIM_ELEM(fourierKernelV,n);
        DIRECT_MULTIDIM_ELEM(transformerPaddedx.fFourier,n)*=K;
    }

    // Inverse Fourier transform
    transformerPaddedx.inverseFourierTransform();
    CenterFFT(paddedx,false);

    // Crop central region
    AtAx.resize(x);
    for (int k=STARTINGZ(x); k<=FINISHINGZ(x); ++k)
        for (int i=STARTINGY(x); i<=FINISHINGY(x); ++i)
            memcpy(&A3D_ELEM(AtAx,k,i,STARTINGX(x)),&A3D_ELEM(paddedx,k,i,STARTINGX(x)),XSIZE(x)*sizeof(double));
}

applyLFilter()

void ProgReconsADMM::applyLFilter	(	MultidimArray< std::complex< double > > &	fourierL,
		bool	adjoint = false
	)

Apply Lt filter

Definition at line 525 of file reconstruct_ADMM.cpp.

{
    transformerL.FourierTransform();
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(fourierL)
        DIRECT_MULTIDIM_ELEM(transformerL.fFourier,n)*=DIRECT_MULTIDIM_ELEM(fourierL,n);
    transformerL.inverseFourierTransform();
    CenterFFT(ud,false);
    if (adjoint)
        ud*=-1;
}

applyLtFilter()

void ProgReconsADMM::applyLtFilter	(	MultidimArray< std::complex< double > > &	fourierL,
		MultidimArray< double > &	u,
		MultidimArray< double > &	d
	)

Definition at line 536 of file reconstruct_ADMM.cpp.

{
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(ud)
        DIRECT_MULTIDIM_ELEM(ud,n)=DIRECT_MULTIDIM_ELEM(u,n)-DIRECT_MULTIDIM_ELEM(d,n);
    applyLFilter(fourierL,true);
}

computeHtKH()

void ProgReconsADMM::computeHtKH

(

MultidimArray< double > &

kernelV

)

Compute H^t*K*H. H is the projection operator, K is the CTF operator

Definition at line 367 of file reconstruct_ADMM.cpp.

{
    kernelV.initZeros(2*ZSIZE(CHtb())-1,2*YSIZE(CHtb())-1,2*XSIZE(CHtb())-1);
    kernelV.setXmippOrigin();

    if (rank==0)
    {
        std::cerr << "Calculating H'KH ...\n";
        init_progress_bar(mdIn.size());
    }
    size_t i=0;
    double rot, tilt, psi, weight=1;
    bool hasWeight=mdIn.containsLabel(MDL_WEIGHT) && useWeights;
    bool hasCTF=(mdIn.containsLabel(MDL_CTF_MODEL) || mdIn.containsLabel(MDL_CTF_DEFOCUSU)) && useCTF;
    CTFDescription ctf;
    MultidimArray<double> kernelAutocorr;
    if (!hasCTF)
        kernel.getKernelAutocorrelation(kernelAutocorr);
    Matrix2D<double> E;
    Matrix1D<double> r1(3), r2(3);
    for (size_t objId : mdIn.ids())
    {
        if ((i+1)%Nprocs==rank)
        {
            // COSS: Read also MIRROR
            mdIn.getValue(MDL_ANGLE_ROT,rot,objId);
            mdIn.getValue(MDL_ANGLE_TILT,tilt,objId);
            mdIn.getValue(MDL_ANGLE_PSI,psi,objId);
            if (hasWeight)
                mdIn.getValue(MDL_WEIGHT,weight,objId);
            if (hasCTF)
            {
                ctf.readFromMetadataRow(mdIn,objId);
                ctf.produceSideInfo();
                kernel.applyCTFToKernelAutocorrelation(ctf,Ts,kernelAutocorr);
            }
            kernelAutocorr.setXmippOrigin();

            // Update kernel
            Euler_angles2matrix(rot,tilt,psi,E,false);
            E.getRow(0,r1);
            E.getRow(1,r2);
            double iStep=1.0/kernel.autocorrStep;
            for (auto k=(kernelV.zinit); k<=(kernelV.zinit + kernelV.zdim - 1); ++k)
            {
                double r1_z=k*ZZ(r1);
                double r2_z=k*ZZ(r2);
                for (auto i=(kernelV.yinit); i<=(kernelV.yinit + kernelV.ydim - 1); ++i)
                {
                    double r1_yz=i*YY(r1)+r1_z;
                    double r2_yz=i*YY(r2)+r2_z;
                    for (auto j=(kernelV.xinit); j<=(kernelV.xinit + kernelV.xdim - 1); ++j)
                    {
                        double r1_xyz=j*XX(r1)+r1_yz;
                        double r2_xyz=j*XX(r2)+r2_yz;
                        A3D_ELEM(kernelV,k,i,j)+=weight*kernelAutocorr.interpolatedElement2D(r1_xyz*iStep,r2_xyz*iStep);
                    }
                }
            }
        }

        i++;
        if (i%100==0 && rank==0)
            progress_bar(i);
    }
    if (rank==0)
        progress_bar(mdIn.size());

#ifndef SYMMETRIZE_PROJECTIONS
    MultidimArray<double> kernelVsym;
    symmetrizeVolume(SL, kernelV, kernelVsym, false, false, true);
    kernelV=kernelVsym;
#endif

    shareVolume(kernelV);
}

constructHtb()

void ProgReconsADMM::constructHtb

(

)

H^t*b

Definition at line 243 of file reconstruct_ADMM.cpp.

{
    size_t xdim, ydim, zdim, ndim;
    getImageSize(mdIn,xdim,ydim,zdim,ndim);
    CHtb().initZeros(xdim,xdim,xdim);
    CHtb().setXmippOrigin();

    Image<double> I;
    double rot, tilt, psi;
    size_t i=0;
    if (rank==0)
    {
        std::cerr << "Performing Htb ...\n";
        init_progress_bar(mdIn.size());
    }
    double weight=1.;
    ApplyGeoParams geoParams;
    geoParams.only_apply_shifts=true;
    geoParams.wrap=xmipp_transformation::DONT_WRAP;
    for (size_t objId : mdIn.ids())
    {
        if ((i+1)%Nprocs==rank)
        {
            I.readApplyGeo(mdIn,objId,geoParams);
            I().setXmippOrigin();
            mdIn.getValue(MDL_ANGLE_ROT,rot,objId);
            mdIn.getValue(MDL_ANGLE_TILT,tilt,objId);
            mdIn.getValue(MDL_ANGLE_PSI,psi,objId);
            if (useWeights && mdIn.containsLabel(MDL_WEIGHT))
                mdIn.getValue(MDL_WEIGHT,weight,objId);

            project(rot,tilt,psi,I(),true,weight);
        }
        i++;
        if (i%100==0 && rank==0)
            progress_bar(i);
    }
    if (rank==0)
        progress_bar(mdIn.size());

    // Symmetrize Htb
#ifndef SYMMETRIZE_PROJECTIONS
    MultidimArray<double> CHtbsym;
    symmetrizeVolume(SL, CHtb(), CHtbsym, false, false, true);
    CHtb=CHtbsym;
#endif
    shareVolume(CHtb());
}

defineParams()

void ProgReconsADMM::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 45 of file reconstruct_ADMM.cpp.

{
    addUsageLine("Reconstruct with Alternative Direction Method of Multipliers");
    addParamsLine(" -i <metadata>: Input metadata with images and angles");
    addParamsLine(" [--oroot <root=reconstruct_admm>]: Rootname for output files");
    addParamsLine(" [--Htb <volume=\"\">]: Filename of the Htb volume");
    addParamsLine(" [--HtKH <volume=\"\">]: Filename of the HtKH volume");
    addParamsLine(" [--firstVolume <volume=\"\">]: Filename of the first initial guess");
    addParamsLine(" [--kernel <shape=KaiserBessel>]: Kernel shape");
    addParamsLine("    where <shape>");
    addParamsLine("      KaiserBessel: Kaiser-Bessel function as kernel");
    addParamsLine(" [--downsamplingV <Ti=1>]: Downsampling factor for volume");
    addParamsLine(" [--downsamplingI <Tp=1>]: Downsampling factor for projections");
    addParamsLine(" [--dontUseWeights]: Do not use weights if available in the input metadata");
    addParamsLine(" [--dontUseCTF]: Do not use CTF if available in the input metadata");
    addParamsLine("               : If the CTF information is used, the algorithm assumes that the input images are phase flipped");
    addParamsLine(" [--sampling <Ts=1>]: The sampling rate is only used to correct for the CTF");
    addParamsLine(" [--mu <mu=1e-4>]: Augmented Lagrange penalty");
    addParamsLine(" [--lambda <lambda=1e-5>]: Total variation parameter");
    addParamsLine(" [--lambda1 <lambda1=1e-7>]: Tikhonov parameter");
    addParamsLine(" [--cgiter <N=3>]: Conjugate Gradient iterations");
    addParamsLine(" [--admmiter <N=30>]: ADMM iterations");
    addParamsLine(" [--positivity]: Positivity constraint");
    addParamsLine(" [--sym <s=c1>]: Symmetry constraint");
    addParamsLine(" [--saveIntermediate]: Save Htb and HtKH volumes for posterior calls");

    mask.defineParams(this,INT_MASK);
}

doPOCSProjection()

void ProgReconsADMM::doPOCSProjection

(

)

POCS projection

Definition at line 606 of file reconstruct_ADMM.cpp.

{
    if (applyMask)
        mask.apply_mask(Ck(),Ck());

    if (positivity)
    {
        Ck().threshold("below",0.,0.);
#ifdef NEVERDEFINED
        MultidimArray<double> smallKernel3D(generate_mask2*a+1,2*a+1,2*a+1);
        smallKernel3D.setXmippOrigin();
        kernel.computeKernel3D(smallKernel3D);
        smallKernel3D/=sqrt(smallKernel3D.sum2());

        produceVolume();
        MultidimArray<double> &mVk=Vk();
        MultidimArray<double> &mCk=Ck();
        Ck.write("PPPCkBefore.vol");
        int k0=STARTINGZ(mVk)+a;
        int i0=STARTINGY(mVk)+a;
        int j0=STARTINGX(mVk)+a;
        int kF=FINISHINGZ(mVk)-a;
        int iF=FINISHINGY(mVk)-a;
        int jF=FINISHINGX(mVk)-a;
        FOR_ALL_ELEMENTS_IN_ARRAY3D(mVk)
            if (k<k0 || k>kF || i<i0 || i>iF || j<j0 || j>jF)
                A3D_ELEM(mCk,k,i,j)=0.;
            else if (A3D_ELEM(mVk,k,i,j)<0)
            {
                // Calculate dot product
                double dot=0.;
                for (int kk=-a; kk<=a; ++kk)
                    for (int ii=-a; ii<=a; ++ii)
                        for (int jj=-a; jj<=a; ++jj)
                            dot+=A3D_ELEM(mCk,k+kk,i+ii,j+jj)*A3D_ELEM(smallKernel3D,kk,ii,jj);

                // Update C
                for (int kk=-a; kk<=a; ++kk)
                    for (int ii=-a; ii<=a; ++ii)
                        for (int jj=-a; jj<=a; ++jj)
                            A3D_ELEM(mCk,k+kk,i+ii,j+jj)-=dot*A3D_ELEM(smallKernel3D,kk,ii,jj);
            }
        Ck.write("PPPCkAfter.vol");
#endif
    }
}

produceSideInfo()

void ProgReconsADMM::produceSideInfo

(

)

Definition at line 106 of file reconstruct_ADMM.cpp.

{
    // Read input images
    mdIn.read(fnIn);

    // Symmetrize input images
#ifdef SYMMETRIZE_PROJECTIONS
    Matrix2D<double> Lsym(3,3), Rsym(3,3);
    MetaDataVec mdSym;
    double rot, tilt, psi, newrot, newtilt, newpsi;
    for (auto& row : mdIn)
    {
        row.getValue(MDL_ANGLE_ROT,rot);
        row.getValue(MDL_ANGLE_TILT,tilt);
        row.getValue(MDL_ANGLE_PSI,psi);
        mdSym.addRow(dynamic_cast<MDRowVec&>(row));
        for (int isym = 0; isym < SL.symsNo(); isym++)
        {
            SL.getMatrices(isym, Lsym, Rsym, false);
            Euler_apply_transf(Lsym,Rsym,rot,tilt,psi,newrot,newtilt,newpsi);
            row.setValue(MDL_ANGLE_ROT,newrot);
            row.setValue(MDL_ANGLE_TILT,newtilt);
            row.setValue(MDL_ANGLE_PSI,newpsi);
            mdSym.addRow(dynamic_cast<MDRowVec&>(row));
        }
    }
    mdIn=mdSym;
    mdSym.clear();
#endif

    // Prepare kernel
    if (kernelShape=="KaiserBessel")
        kernel.initializeKernel(alpha,a,0.0001);
    kernel.convolveKernelWithItself(0.005);

    // Get Htb reconstruction
    if (fnHtb=="")
    {
        constructHtb();
        if (saveIntermediate && rank==0)
            CHtb.write(fnRoot+"_Htb.vol");
    }
    else
    {
        CHtb.read(fnHtb);
        CHtb().setXmippOrigin();
    }

    // Adjust regularization weights
    //double Nimgs=mdIn.size();
    //double CHtb_norm=sqrt(CHtb().sum2()/MULTIDIM_SIZE(CHtb()));
    //COSS mu*=Nimgs/XSIZE(CHtb())*CHtb_norm;
    //COSS lambda*=Nimgs/XSIZE(CHtb())*CHtb_norm;
    //COSS lambda1*=Nimgs/XSIZE(CHtb());

    // Get first volume
    if (fnFirst=="")
        Ck().initZeros(CHtb());
    else
    {
        Ck.read(fnFirst);
        Ck().setXmippOrigin();
    }

    // Compute H'*K*H+mu*L^T*L+lambda_1 I
    Image<double> kernelV;
    if (fnHtKH=="")
    {
        computeHtKH(kernelV());
        if (saveIntermediate && rank==0)
            kernelV.write(fnRoot+"_HtKH.vol");
    }
    else
    {
        kernelV.read(fnHtKH);
        kernelV().setXmippOrigin();
    }
    FourierTransformer transformer;
    transformer.FourierTransform(kernelV(),fourierKernelV,true);

    // Add regularization in Fourier space
    addRegularizationTerms();

    // Resize u and d volumes
    ux.initZeros(CHtb());
    ux.setXmippOrigin();
    uy=ux;
    uz=ux;
    dx=ux;
    dy=ux;
    dz=ux;

    // Prepare Lt filters
    MultidimArray<double> L;
    L.resize(ux);
    kernel.computeGradient(L,'x');
    transformer.FourierTransform(L,fourierLx);
    kernel.computeGradient(L,'y');
    transformer.FourierTransform(L,fourierLy);
    kernel.computeGradient(L,'z');
    transformer.FourierTransform(L,fourierLz);

    ud=ux;
    transformerL.setReal(ud);

    // Prepare mask
    if (applyMask)
        mask.generate_mask(CHtb());
    synchronize();
}

produceVolume()

void ProgReconsADMM::produceVolume

(

)

Convert from coefficients to voxel volume

Definition at line 679 of file reconstruct_ADMM.cpp.

{
    MultidimArray<double> kernel3D;
    kernel3D.resize(Ck());
    kernel.computeKernel3D(kernel3D);
    convolutionFFT(Ck(),kernel3D,Vk());
    // Vk()*=kernel3D.sum();
}

project() [1/2]

void ProgReconsADMM::project	(	double	rot,
		double	tilt,
		double	psi,
		MultidimArray< double > &	P,
		bool	adjoint = false,
		double	weight = 1.
	)

Project the volume V onto P using rot,tilt,psi

Definition at line 292 of file reconstruct_ADMM.cpp.

{
    Matrix2D<double> E;
    Euler_angles2matrix(rot,tilt,psi,E,false);
    Matrix1D<double> r1(3), r2(3);
    E.getRow(0,r1);
    E.getRow(1,r2);
    project(r1,r2,P,adjoint,weight);
}

project() [2/2]

void ProgReconsADMM::project	(	const Matrix1D< double > &	r1,
		const Matrix1D< double > &	r2,
		MultidimArray< double > &	P,
		bool	adjoint = false,
		double	weight = 1.
	)

Project the volume V onto P using r1 and r2 as the coordinate system

Definition at line 302 of file reconstruct_ADMM.cpp.

{
    const MultidimArray<double> &mV=CHtb();
    if (!adjoint)
    {
        P.initZeros(std::ceil(XSIZE(mV)/Tp),std::ceil(YSIZE(mV)/Tp));
        P.setXmippOrigin();
    }

    // Tomographic projection
    for (int k=STARTINGZ(mV); k<=FINISHINGZ(mV); ++k) {
        // initial computation
        double rzn  = k;
        double rzn1= rzn*VEC_ELEM(r1,2);
        double rzn2= rzn*VEC_ELEM(r2,2);
        for (int i=STARTINGY(mV); i<=FINISHINGY(mV); ++i) {
            // initial computation
            double ryn  = i;
            double ryn1= ryn*VEC_ELEM(r1,1);
            double ryn2= ryn*VEC_ELEM(r2,1);
            double sx0=rzn1+ryn1;
            double sy0=rzn2+ryn2;

            for (int j=STARTINGX(mV); j<=FINISHINGX(mV); ++j) {
                double rxn  = j;
                double rxn1= rxn*VEC_ELEM(r1,0);
                double rxn2= rxn*VEC_ELEM(r2,0);

                double sx = Ti*(sx0+rxn1);
                double sy = Ti*(sy0+rxn2);

                int sxmin = std::floor(sx-kernel.supp);
                int sxmax = std::ceil(sx+kernel.supp);
                int symin = std::floor(sy-kernel.supp);
                int symax = std::ceil(sy+kernel.supp);
                if (sxmin<STARTINGX(P))
                    sxmin = STARTINGX(P);
                if (sxmax>FINISHINGX(P))
                    sxmax = FINISHINGX(P);
                if (symin<STARTINGY(P))
                    symin = STARTINGY(P);
                if (symax>FINISHINGY(P))
                    symax = FINISHINGY(P);

                if (adjoint)
                    for (int ii=symin; ii<=symax; ii++) {
                        double u=(ii-sy)*Tp;
                        for (int jj=sxmin; jj<=sxmax; jj++) {
                            double v=(jj-sx)*Tp;
                            A3D_ELEM(mV,k,i,j)+=weight*A2D_ELEM(P,ii,jj)*kernel.projectionValueAt(u,v);
                        }
                    }
                else
                    for (int ii=symin; ii<=symax; ii++) {
                        double u=(ii-sy)*Tp;
                        for (int jj=sxmin; jj<=sxmax; jj++) {
                            double v=(jj-sx)*Tp;
                            A2D_ELEM(P,ii,jj)+=weight*A3D_ELEM(mV,k,i,j)*kernel.projectionValueAt(u,v);
                        }
                    }
            }
        }
    }
}

readParams()

void ProgReconsADMM::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 74 of file reconstruct_ADMM.cpp.

{
    fnIn=getParam("-i");
    fnRoot=getParam("--oroot");
    fnFirst=getParam("--firstVolume");
    fnHtb=getParam("--Htb");
    fnHtKH=getParam("--HtKH");
    kernelShape=getParam("--kernel");
    if (kernelShape=="KaiserBessel")
    {
        a=4; //getDoubleParam("--kernel",1);
        alpha=19; // getDoubleParam("--kernel",2);
    }
    Ti=getDoubleParam("--downsamplingV");
    Tp=getDoubleParam("--downsamplingI");
    useWeights=!checkParam("--dontUseWeights");
    useCTF=!checkParam("--dontUseCTF");
    Ts=getDoubleParam("--sampling");
    mu=getDoubleParam("--mu");
    lambda=getDoubleParam("--lambda");
    lambda1=getDoubleParam("--lambda1");
    Ncgiter=getIntParam("--cgiter");
    Nadmmiter=getIntParam("--admmiter");
    positivity=checkParam("--positivity");
    saveIntermediate=checkParam("--saveIntermediate");

    applyMask=checkParam("--mask");
    if (applyMask)
        mask.readParams(this);
    SL.readSymmetryFile(getParam("--sym"));
}

run()

void ProgReconsADMM::run ( )

virtual

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

Reimplemented from XmippProgram.

Definition at line 222 of file reconstruct_ADMM.cpp.

{
    produceSideInfo();
    if (rank==0)
    {
        std::cout << "Running ADMM iterations ..." << std::endl;
        init_progress_bar(Nadmmiter);
        for (int iter=0; iter<Nadmmiter; ++iter)
        {
            applyConjugateGradient();
            doPOCSProjection();
            updateUD();
            progress_bar(iter);
        }
        progress_bar(Nadmmiter);
        produceVolume();
        Vk.write(fnRoot+".vol");
    }
    synchronize();
}

shareVolume()

virtual void ProgReconsADMM::shareVolume ( MultidimArray< double > &  V )

inlinevirtual

Reimplemented in MpiProgReconstructADMM.

Definition at line 145 of file reconstruct_ADMM.h.

{}

show()

void ProgReconsADMM::show

(

)

Definition at line 217 of file reconstruct_ADMM.cpp.

{

}

synchronize()

virtual void ProgReconsADMM::synchronize ( )

inlinevirtual

Reimplemented in MpiProgReconstructADMM.

Definition at line 148 of file reconstruct_ADMM.h.

{}

updateUD()

void ProgReconsADMM::updateUD

(

)

Update u and d

Definition at line 653 of file reconstruct_ADMM.cpp.

{
    double threshold;
    if (mu>0)
        threshold=lambda/mu;
    else
        threshold=0;
//  std::cout << "lambda=" << lambda << std::endl;
//  std::cout << "mu=" << mu << std::endl;
//  std::cout << "Threshold=" << threshold << std::endl;
    MultidimArray<double> LCk;

    // Update gradient and Lagrange parameters
    ud=Ck(); applyLFilter(fourierLx); LCk=ud;
    ux=LCk; ux+=dx; ux.threshold("soft",threshold);
    dx+=LCk; dx-=ux;

    ud=Ck(); applyLFilter(fourierLy); LCk=ud;
    uy=LCk; uy+=dy; uy.threshold("soft",threshold);
    dy+=LCk; dy-=uy;

    ud=Ck(); applyLFilter(fourierLz); LCk=ud;
    uz=LCk; uz+=dz; uz.threshold("soft",threshold);
    dz+=LCk; dz-=uz;
}

Member Data Documentation

a

double ProgReconsADMM::a

Definition at line 78 of file reconstruct_ADMM.h.

alpha

double ProgReconsADMM::alpha

Definition at line 79 of file reconstruct_ADMM.h.

applyMask

bool ProgReconsADMM::applyMask

Definition at line 91 of file reconstruct_ADMM.h.

CHtb

Image<double> ProgReconsADMM::CHtb

Definition at line 152 of file reconstruct_ADMM.h.

Ck

Image<double> ProgReconsADMM::Ck

Definition at line 153 of file reconstruct_ADMM.h.

dx

MultidimArray<double> ProgReconsADMM::dx

Definition at line 159 of file reconstruct_ADMM.h.

dy

MultidimArray<double> ProgReconsADMM::dy

Definition at line 159 of file reconstruct_ADMM.h.

dz

MultidimArray<double> ProgReconsADMM::dz

Definition at line 159 of file reconstruct_ADMM.h.

fnFirst

FileName ProgReconsADMM::fnFirst

Definition at line 74 of file reconstruct_ADMM.h.

fnHtb

FileName ProgReconsADMM::fnHtb

Definition at line 75 of file reconstruct_ADMM.h.

fnHtKH

FileName ProgReconsADMM::fnHtKH

Definition at line 76 of file reconstruct_ADMM.h.

fnIn

FileName ProgReconsADMM::fnIn

Definition at line 72 of file reconstruct_ADMM.h.

fnRoot

FileName ProgReconsADMM::fnRoot

Definition at line 73 of file reconstruct_ADMM.h.

fourierKernelV

MultidimArray<std::complex<double> > ProgReconsADMM::fourierKernelV

Definition at line 155 of file reconstruct_ADMM.h.

fourierLx

MultidimArray< std::complex<double> > ProgReconsADMM::fourierLx

Definition at line 160 of file reconstruct_ADMM.h.

fourierLy

MultidimArray< std::complex<double> > ProgReconsADMM::fourierLy

Definition at line 160 of file reconstruct_ADMM.h.

fourierLz

MultidimArray< std::complex<double> > ProgReconsADMM::fourierLz

Definition at line 160 of file reconstruct_ADMM.h.

kernel

AdmmKernel ProgReconsADMM::kernel

Definition at line 151 of file reconstruct_ADMM.h.

kernelShape

String ProgReconsADMM::kernelShape

Definition at line 77 of file reconstruct_ADMM.h.

lambda

double ProgReconsADMM::lambda

Definition at line 81 of file reconstruct_ADMM.h.

lambda1

double ProgReconsADMM::lambda1

Definition at line 82 of file reconstruct_ADMM.h.

mask

Mask ProgReconsADMM::mask

Definition at line 93 of file reconstruct_ADMM.h.

mdIn

MetaDataVec ProgReconsADMM::mdIn

Definition at line 154 of file reconstruct_ADMM.h.

mu

double ProgReconsADMM::mu

Definition at line 80 of file reconstruct_ADMM.h.

Nadmmiter

int ProgReconsADMM::Nadmmiter

Definition at line 89 of file reconstruct_ADMM.h.

Ncgiter

int ProgReconsADMM::Ncgiter

Definition at line 88 of file reconstruct_ADMM.h.

Nprocs

size_t ProgReconsADMM::Nprocs

Definition at line 96 of file reconstruct_ADMM.h.

paddedx

MultidimArray<double> ProgReconsADMM::paddedx

Definition at line 156 of file reconstruct_ADMM.h.

positivity

bool ProgReconsADMM::positivity

Definition at line 90 of file reconstruct_ADMM.h.

rank

size_t ProgReconsADMM::rank

Definition at line 97 of file reconstruct_ADMM.h.

saveIntermediate

bool ProgReconsADMM::saveIntermediate

Definition at line 95 of file reconstruct_ADMM.h.

SL

SymList ProgReconsADMM::SL

Definition at line 161 of file reconstruct_ADMM.h.

symmetry

String ProgReconsADMM::symmetry

Definition at line 94 of file reconstruct_ADMM.h.

Ti

double ProgReconsADMM::Ti

Definition at line 84 of file reconstruct_ADMM.h.

Tp

double ProgReconsADMM::Tp

Definition at line 85 of file reconstruct_ADMM.h.

transformerL

FourierTransformer ProgReconsADMM::transformerL

Definition at line 157 of file reconstruct_ADMM.h.

transformerPaddedx

FourierTransformer ProgReconsADMM::transformerPaddedx

Definition at line 157 of file reconstruct_ADMM.h.

Ts

double ProgReconsADMM::Ts

Definition at line 92 of file reconstruct_ADMM.h.

ud

MultidimArray<double> ProgReconsADMM::ud

Definition at line 159 of file reconstruct_ADMM.h.

useCTF

bool ProgReconsADMM::useCTF

Definition at line 87 of file reconstruct_ADMM.h.

useWeights

bool ProgReconsADMM::useWeights

Definition at line 86 of file reconstruct_ADMM.h.

ux

MultidimArray<double> ProgReconsADMM::ux

Definition at line 159 of file reconstruct_ADMM.h.

uy

MultidimArray<double> ProgReconsADMM::uy

Definition at line 159 of file reconstruct_ADMM.h.

uz

MultidimArray<double> ProgReconsADMM::uz

Definition at line 159 of file reconstruct_ADMM.h.

Vk

Image<double> ProgReconsADMM::Vk

Definition at line 153 of file reconstruct_ADMM.h.

---

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

• xmipp/legacy/libraries/reconstruction/reconstruct_ADMM.h
• xmipp/legacy/libraries/reconstruction/reconstruct_ADMM.cpp