#include <angular_sph_alignment.h>

Inheritance diagram for ProgAngularSphAlignment:

Collaboration diagram for ProgAngularSphAlignment:

Public Member Functions
	ProgAngularSphAlignment ()
	Empty constructor. More...

	~ProgAngularSphAlignment ()
	Destructor. More...

void	readParams ()
	Read argument from command line. More...

void	show ()
	Show. More...

void	defineParams ()
	Define parameters. More...

void	preProcess ()

void	processImage (const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut)

void	numCoefficients (int l1, int l2, int &nc) const
	Length of coefficients vector. More...

void	minimizepos (int l2, Matrix1D< double > &steps) const
	Determine the positions to be minimize of a vector containing spherical harmonic coefficients. More...

void	fillVectorTerms (int l1, int l2)
	Zernike and SPH coefficients allocation. More...

void	deformVol (MultidimArray< double > &mVD, const MultidimArray< double > &mV, double &def, double rot, double tilt, double psi)
	Deform a volumen using Zernike-Spherical harmonic basis. More...

void	applyCTFImage (double const &deltaDefocusU, double const &deltaDefocusV, double const &deltaDefocusAngle)

void	updateCTFImage (double defocusU, double defocusV, double angle)

double	tranformImageSph (double *pclnm, double rot, double tilt, double psi, double deltaDefocusU, double deltaDefocusV, double deltaDefocusAngle)

virtual void	finishProcessing ()

virtual void	writeImageParameters (const FileName &fnImg)

Public Member Functions inherited from XmippMetadataProgram
MetaData *	getInputMd ()

MetaDataVec &	getOutputMd ()

	XmippMetadataProgram ()
	Empty constructor. More...

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

virtual void	init ()

virtual void	setup (MetaData *md, const FileName &o="", const FileName &oroot="", bool applyGeo=false, MDLabel label=MDL_IMAGE)

virtual	~XmippMetadataProgram ()

void	setMode (WriteModeMetaData _mode)

void	setupRowOut (const FileName &fnImgIn, const MDRow &rowIn, const FileName &fnImgOut, MDRow &rowOut) const
	Prepare rowout. More...

virtual void	wait ()
	Wait for the distributor to finish. More...

virtual void	checkPoint ()
	For very long programs, it may be needed to write checkpoints. More...

virtual void	run ()
	Run over all images. 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)

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	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	fnVolR

FileName	fnMaskR

FileName	fnOutDir
	Output directory. More...

int	L1

int	L2

Matrix1D< int >	vL1

Matrix1D< int >	vN

Matrix1D< int >	vL2

Matrix1D< int >	vM

double	maxShift

double	maxAngularChange

double	maxResol

double	Ts

int	Rmax

bool	optimizeAlignment

bool	optimizeDeformation

bool	optimizeDefocus

bool	ignoreCTF

bool	optimizeRadius

bool	phaseFlipped

double	lambda

int	RmaxDef

Matrix1D< double >	p

int	flagEnabled

bool	resume

MultidimArray< int >	mask2D

MultidimArray< int >	V_mask

size_t	Xdim

Image< double >	V

Image< double >	Vdeformed

Image< double >	I

Image< double >	Ip

Image< double >	Ifiltered

Image< double >	Ifilteredp

Image< double >	P

FourierFilter	filter

Matrix2D< double >	A

double	old_rot

double	old_tilt

double	old_psi

double	old_shiftX

double	old_shiftY

bool	old_flip

bool	hasCTF

double	old_defocusU

double	old_defocusV

double	old_defocusAngle

double	currentDefocusU

double	currentDefocusV

double	currentAngle

CTFDescription	ctf

FourierFilter	FilterCTF

int	vecSize

Matrix1D< double >	clnm

Matrix1D< double >	steps_cp

double	totalDeformation

double	sumV

double	sumVd

bool	showOptimization

double	correlation

Public Attributes inherited from XmippMetadataProgram
FileName	fn_in
	Filenames of input and output Metadata. More...

FileName	fn_out

FileName	baseName

FileName	pathBaseName

FileName	oextBaseName

bool	apply_geo
	Apply geo. More...

size_t	ndimOut
	Output dimensions. More...

size_t	zdimOut

size_t	ydimOut

size_t	xdimOut

DataType	datatypeOut

size_t	mdInSize
	Number of input elements. More...

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Protected Member Functions
virtual void	createWorkFiles ()

Protected Member Functions inherited from XmippMetadataProgram
virtual void	initComments ()

virtual void	postProcess ()

virtual bool	getImageToProcess (size_t &objId, size_t &objIndex)

void	show () const override

virtual void	startProcessing ()

virtual void	writeOutput ()

virtual void	showProgress ()

virtual void	defineLabelParam ()

Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Member Functions inherited from Rerunable
	Rerunable (const FileName &fn)

virtual void	createWorkFiles (bool resume, MetaData *md)

const FileName &	getFileName () const

void	setFileName (const FileName &fn)

Additional Inherited Members
Protected Attributes inherited from XmippMetadataProgram
WriteModeMetaData	mode
	Metadata writing mode: OVERWRITE, APPEND. More...

FileName	oext
	Output extension and root. More...

FileName	oroot

MDLabel	image_label
	MDLabel to be used to read/write images, usually will be MDL_IMAGE. More...

bool	produces_an_output
	Indicate that a unique final output is produced. More...

bool	produces_a_metadata
	Indicate that the unique final output file is a Metadata. More...

bool	each_image_produces_an_output
	Indicate that an output is produced for each image in the input. More...

bool	allow_apply_geo

bool	decompose_stacks
	Input Metadata will treat a stack file as a set of images instead of a unique file. More...

bool	delete_output_stack
	Delete previous output stack file prior to process images. More...

bool	get_image_info
	Get the input image file dimensions to further operations. More...

bool	save_metadata_stack
	Save the associated output metadata when output file is a stack. More...

bool	track_origin
	Include the original input image filename in the output stack. More...

bool	keep_input_columns
	Keep input metadata columns. More...

bool	remove_disabled
	Remove disabled images from the input selfile. More...

bool	allow_time_bar
	Show process time bar. More...

bool	input_is_metadata
	Input is a metadata. More...

bool	single_image
	Input is a single image. More...

bool	input_is_stack
	Input is a stack. More...

bool	output_is_stack
	Output is a stack. More...

bool	create_empty_stackfile

bool	delete_mdIn

size_t	time_bar_step
	Some time bar related counters. More...

size_t	time_bar_size

size_t	time_bar_done

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

Predict Continuous Parameters.

Definition at line 42 of file angular_sph_alignment.h.

Constructor & Destructor Documentation

◆ ProgAngularSphAlignment()

ProgAngularSphAlignment::ProgAngularSphAlignment ( )

Empty constructor.

Definition at line 35 of file angular_sph_alignment.cpp.

                                                  : Rerunable("")
 {
     resume = false;
     produces_a_metadata = true;
     each_image_produces_an_output = false;
     showOptimization = false;
 }

◆ ~ProgAngularSphAlignment()

ProgAngularSphAlignment::~ProgAngularSphAlignment ( )

default

Destructor.

Member Function Documentation

◆ applyCTFImage()

void ProgAngularSphAlignment::applyCTFImage	(	double const &	deltaDefocusU,
		double const &	deltaDefocusV,
		double const &	deltaDefocusAngle
	)

Definition at line 506 of file angular_sph_alignment.cpp.

                                                                              {
     double defocusU=old_defocusU+deltaDefocusU;
     double defocusV=old_defocusV+deltaDefocusV;
     double angle=old_defocusAngle+deltaDefocusAngle;
     if (defocusU!=currentDefocusU || defocusV!=currentDefocusV || angle!=currentAngle) {
         updateCTFImage(defocusU,defocusV,angle);
     }
     FilterCTF.ctf = ctf;
     FilterCTF.generateMask(P());
     if (phaseFlipped)
         FilterCTF.correctPhase();
     FilterCTF.applyMaskSpace(P());
 } 

◆ createWorkFiles()

virtual void ProgAngularSphAlignment::createWorkFiles ( )

inlineprotectedvirtual

Reimplemented in MpiProgAngularSphAlignment.

Definition at line 201 of file angular_sph_alignment.h.

                                    {
       return Rerunable::createWorkFiles(resume, getInputMd());
     }

◆ defineParams()

void ProgAngularSphAlignment::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippMetadataProgram.

Definition at line 96 of file angular_sph_alignment.cpp.

 {
     addUsageLine("Make a continuous angular assignment with deformations");
     defaultComments["-i"].clear();
     defaultComments["-i"].addComment("Metadata with initial alignment");
     defaultComments["-o"].clear();
     defaultComments["-o"].addComment("Metadata with the angular alignment and deformation parameters");
     XmippMetadataProgram::defineParams();
     addParamsLine("   --ref <volume>              : Reference volume");
     addParamsLine("  [--mask <m=\"\">]            : Reference volume mask");
     addParamsLine("  [--odir <outputDir=\".\">]   : Output directory");
     addParamsLine("  [--max_shift <s=-1>]         : Maximum shift allowed in pixels");
     addParamsLine("  [--max_angular_change <a=5>] : Maximum angular change allowed (in degrees)");
     addParamsLine("  [--max_resolution <f=4>]     : Maximum resolution (A)");
     addParamsLine("  [--sampling <Ts=1>]          : Sampling rate (A/pixel)");
     addParamsLine("  [--Rmax <R=-1>]              : Maximum radius (px). -1=Half of volume size");
     addParamsLine("  [--RDef <r=-1>]              : Maximum radius of the deformation (px). -1=Half of volume size");
     addParamsLine("  [--l1 <l1=3>]                : Degree Zernike Polynomials=1,2,3,...");
     addParamsLine("  [--l2 <l2=2>]                : Harmonical depth of the deformation=1,2,3,...");
     addParamsLine("  [--optimizeAlignment]        : Optimize alignment");
     addParamsLine("  [--optimizeDeformation]      : Optimize deformation");
     addParamsLine("  [--optimizeDefocus]          : Optimize defocus");
     addParamsLine("  [--phaseFlipped]             : Input images have been phase flipped");
     addParamsLine("  [--regularization <l=0.01>]  : Regularization weight");
     addParamsLine("  [--resume]                   : Resume processing");
     addExampleLine("A typical use is:",false);
     addExampleLine("xmipp_angular_sph_alignment -i anglesFromContinuousAssignment.xmd --ref reference.vol -o assigned_anglesAndDeformations.xmd --optimizeAlignment --optimizeDeformation --depth 1");
 }

◆ deformVol()

void ProgAngularSphAlignment::deformVol	(	MultidimArray< double > &	mVD,
		const MultidimArray< double > &	mV,
		double &	def,
		double	rot,
		double	tilt,
		double	psi
	)

Deform a volumen using Zernike-Spherical harmonic basis.

Definition at line 530 of file angular_sph_alignment.cpp.

 {
     size_t idxY0=(VEC_XSIZE(clnm)-8)/3;
     double Ncount=0.0;
     double modg=0.0;
     double diff2=0.0;
 
     def=0.0;
     size_t idxZ0=2*idxY0;
     sumVd=0.0;
     double RmaxF=RmaxDef;
     double RmaxF2=RmaxF*RmaxF;
     double iRmaxF=1.0/RmaxF;
     // Rotation Matrix
     Matrix2D<double> R;
     R.initIdentity(3);
     Euler_angles2matrix(rot, tilt, psi, R, false);
     R = R.inv();
     Matrix1D<double> pos;
     pos.initZeros(3);
 
     // TODO: Poner primero i y j en el loop, acumular suma y guardar al final
     for (int k=STARTINGZ(mV); k<=FINISHINGZ(mV); k++)
     {
         for (int i=STARTINGY(mV); i<=FINISHINGY(mV); i++)
         {
             for (int j=STARTINGX(mV); j<=FINISHINGX(mV); j++)
             {
                 ZZ(pos) = k; YY(pos) = i; XX(pos) = j;
                 pos = R * pos;
                 double gx=0.0;
                 double gy=0.0;
                 double gz=0.0;
                 // TODO: Sacar al bucle de z
                 double k2=ZZ(pos)*ZZ(pos);
                 double kr=ZZ(pos)*iRmaxF;
                 double k2i2=k2+YY(pos)*YY(pos);
                 double ir=YY(pos)*iRmaxF;
                 double r2=k2i2+XX(pos)*XX(pos);
                 double jr=XX(pos)*iRmaxF;
                 double rr=sqrt(r2)*iRmaxF;
                 if (r2<RmaxF2) {
                     for (size_t idx=0; idx<idxY0; idx++) {
                         if (VEC_ELEM(steps_cp,idx) == 1) {
                             double zsph=0.0;
                             int l1 = VEC_ELEM(vL1,idx);
                             int n = VEC_ELEM(vN,idx);
                             int l2 = VEC_ELEM(vL2,idx);
                             int m = VEC_ELEM(vM,idx);
                             zsph=ZernikeSphericalHarmonics(l1,n,l2,m,jr,ir,kr,rr);
                             if (rr>0 || l2==0) {
                                 gx += VEC_ELEM(clnm,idx)        *zsph;
                                 gy += VEC_ELEM(clnm,idx+idxY0)  *zsph;
                                 gz += VEC_ELEM(clnm,idx+idxZ0)  *zsph;
                             }
                         }
                     }
                     auto k_mask = (int)(ZZ(pos)+gz);
                     auto i_mask = (int)(YY(pos)+gy);
                     auto j_mask = (int)(XX(pos)+gx);
                     int voxelI_mask = 0;
                     if (!V_mask.outside(k_mask, i_mask, j_mask)) {
                         voxelI_mask = A3D_ELEM(V_mask, k_mask, i_mask, j_mask);
                     }
                     if (voxelI_mask == 1) {
                         double voxelI=mV.interpolatedElement3D(XX(pos)+gx,YY(pos)+gy,ZZ(pos)+gz);
                         A2D_ELEM(mP,i,j) += voxelI;
                         sumVd += voxelI;
                         modg += gx*gx+gy*gy+gz*gz;
                         Ncount++;
                     }
                 }
             }
         }
     }
 
     def = sqrt(modg/Ncount);
     totalDeformation = def;
 }

◆ fillVectorTerms()

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

Zernike and SPH coefficients allocation.

Definition at line 484 of file angular_sph_alignment.cpp.

 {
     int idx = 0;
     vL1.initZeros(vecSize);
     vN.initZeros(vecSize);
     vL2.initZeros(vecSize);
     vM.initZeros(vecSize);
     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++) {
                 VEC_ELEM(vL1,idx) = l;
                 VEC_ELEM(vN,idx) = h;
                 VEC_ELEM(vL2,idx) = h;
                 VEC_ELEM(vM,idx) = m-aux;
                 idx++;
             }
         }
     }
 }

◆ finishProcessing()

void ProgAngularSphAlignment::finishProcessing ( )

virtual

Write the final parameters.

Reimplemented from XmippMetadataProgram.

Reimplemented in MpiProgAngularSphAlignment.

Definition at line 190 of file angular_sph_alignment.cpp.

                                                {
     XmippMetadataProgram::finishProcessing();
     rename(Rerunable::getFileName().c_str(), fn_out.c_str());
 }

◆ minimizepos()

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

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

Definition at line 472 of file angular_sph_alignment.cpp.

 {
     int size = 0;
     numCoefficients(L1,l2,size);
     auto totalSize = (int)((steps.size()-8)/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 ProgAngularSphAlignment::numCoefficients	(	int	l1,
		int	l2,
		int &	nc
	)		const

Length of coefficients vector.

Definition at line 451 of file angular_sph_alignment.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*(count-numEven);
         }
     }
 }

◆ preProcess()

void ProgAngularSphAlignment::preProcess ( )

virtual

Produce side info. An exception is thrown if any of the files is not found

Reimplemented from XmippMetadataProgram.

Definition at line 125 of file angular_sph_alignment.cpp.

 {
     V.read(fnVolR);
     V().setXmippOrigin();
     Xdim=XSIZE(V());
     Vdeformed().initZeros(V());
 
     Ifilteredp().initZeros(Xdim,Xdim);
     Ifilteredp().setXmippOrigin();
 
     if (RmaxDef<0)
         RmaxDef = Xdim/2;
 
     // Read Reference mask if avalaible (otherwise sphere of radius RmaxDef is used)
     Mask mask;
     mask.type = BINARY_CIRCULAR_MASK;
     mask.mode = INNER_MASK;
     if (fnMaskR != "") {
         Image<double> aux;
         aux.read(fnMaskR);
         typeCast(aux(), V_mask);
         V_mask.setXmippOrigin();
     }
     else {
         mask.R1 = RmaxDef;
         mask.generate_mask(V());
         V_mask = mask.get_binary_mask();
         V_mask.setXmippOrigin();
     }
 
     // Total Volume Mass (Inside Mask)
     sumV = 0.0;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V_mask) {
         if (DIRECT_MULTIDIM_ELEM(V_mask,n) == 1) {
             sumV += DIRECT_MULTIDIM_ELEM(V(),n);
         }
     }
 
     // Construct mask
     if (Rmax<0)
         Rmax=Xdim/2;
     mask.R1 = Rmax;
     mask.generate_mask(Xdim,Xdim);
     mask2D=mask.get_binary_mask();
 
     // Low pass filter
     filter.FilterBand=LOWPASS;
     filter.w1=Ts/maxResol;
     filter.raised_w=0.02;
 
     // Transformation matrix
     A.initIdentity(3);
 
     // CTF Filter
     FilterCTF.FilterBand = CTF;
     FilterCTF.ctf.enable_CTFnoise = false;
     FilterCTF.ctf.produceSideInfo();
 
     vecSize = 0;
     numCoefficients(L1,L2,vecSize);
     fillVectorTerms(L1,L2);
 
     createWorkFiles();
 }

◆ processImage()

void ProgAngularSphAlignment::processImage	(	const FileName &	fnImg,
		const FileName &	fnImgOut,
		const MDRow &	rowIn,
		MDRow &	rowOut
	)

virtual

Predict angles and shift. At the input the pose parameters must have an initial guess of the parameters. At the output they have the estimated pose.

Implements XmippMetadataProgram.

Definition at line 293 of file angular_sph_alignment.cpp.

 {
     Matrix1D<double> steps;
     int totalSize = 3*vecSize+8;
     p.initZeros(totalSize);
     clnm.initZeros(totalSize);
 
     rowOut=rowIn;
 
     flagEnabled=1;
 
     rowIn.getValue(MDL_ANGLE_ROT,old_rot);
     rowIn.getValue(MDL_ANGLE_TILT,old_tilt);
     rowIn.getValue(MDL_ANGLE_PSI,old_psi);
     rowIn.getValue(MDL_SHIFT_X,old_shiftX);
     rowIn.getValue(MDL_SHIFT_Y,old_shiftY);
     if (rowIn.containsLabel(MDL_FLIP))
         rowIn.getValue(MDL_FLIP,old_flip);
     else
         old_flip = false;
     
     if (rowIn.containsLabel(MDL_CTF_DEFOCUSU) || rowIn.containsLabel(MDL_CTF_MODEL))
     {
         hasCTF=true;
         ctf.readFromMdRow(rowIn);
         ctf.Tm = Ts;
         ctf.produceSideInfo();
         old_defocusU=ctf.DeltafU;
         old_defocusV=ctf.DeltafV;
         old_defocusAngle=ctf.azimuthal_angle;
     }
     else
         hasCTF=false;
 
     if (verbose>=2)
         std::cout << "Processing " << fnImg << std::endl;
     I.read(fnImg);
     I().setXmippOrigin();
 
     Ifiltered()=I();
     filter.applyMaskSpace(Ifiltered());
 
     for (int h=1;h<=L2;h++)
     {
         if (verbose>=2)
         {
             std::cout<<std::endl;
             std::cout<<"------------------------------ Basis Degrees: ("<<L1<<","<<h<<") ----------------------------"<<std::endl;
         }
         steps.clear();
         steps.initZeros(totalSize);
 
         // Optimize
         double cost=-1;
         try
         {
             cost=1e38;
             int iter;
             if (optimizeAlignment)
                 steps(totalSize-8)=steps(totalSize-7)=steps(totalSize-6)=steps(totalSize-5)=steps(totalSize-4)=1.;
             if (optimizeDefocus)
                 steps(totalSize-3)=steps(totalSize-2)=steps(totalSize-1)=1.;
             if (optimizeDeformation)
             {
                 minimizepos(h,steps);
             }
             steps_cp = steps;
             powellOptimizer(p, 1, totalSize, &continuousSphCost, this, 0.01, cost, iter, steps, verbose>=2);
 
             if (verbose>=3)
             {
                 showOptimization = true;
                 continuousSphCost(p.adaptForNumericalRecipes(),this);
                 showOptimization = false;
             }
 
             if (cost>0)
             {
                 flagEnabled=-1;
                 p.initZeros();
             }
             cost=-cost;
             correlation=cost;
             if (verbose>=2)
             {
                 std::cout<<std::endl;
                 for (int j=1;j<4;j++)
                 {
                     switch (j)
                     {
                     case 1:
                         std::cout << "X Coefficients=(";
                         break;
                     case 2:
                         std::cout << "Y Coefficients=(";
                         break;
                     case 3:
                         std::cout << "Z Coefficients=(";
                         break;
                     }
                     for (int i=(j-1)*vecSize;i<j*vecSize;i++)
                     {
                         std::cout << p(i);
                         if (i<j*vecSize-1)
                             std::cout << ",";
                     }
                     std::cout << ")" << std::endl;
                 }
                 std::cout << "Radius=" << RmaxDef << std::endl;
                 std::cout << " Dshift=(" << p(totalSize-5) << "," << p(totalSize-4) << ") "
                           << "Drot=" << p(totalSize-3) << " Dtilt=" << p(totalSize-2) 
                           << " Dpsi=" << p(totalSize-1) << std::endl;
                 std::cout << " Total deformation=" << totalDeformation << std::endl;
                 std::cout<<std::endl;
             }
         }
         catch (XmippError &XE)
         {
             std::cerr << XE.what() << std::endl;
             std::cerr << "Warning: Cannot refine " << fnImg << std::endl;
             flagEnabled=-1;
         }
     }
 
         //AJ NEW
         writeImageParameters(fnImg);
         //END AJ
 
 }

◆ readParams()

void ProgAngularSphAlignment::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippMetadataProgram.

Definition at line 46 of file angular_sph_alignment.cpp.

 {
     XmippMetadataProgram::readParams();
     fnVolR = getParam("--ref");
     fnMaskR = getParam("--mask");
     fnOutDir = getParam("--odir");
     maxShift = getDoubleParam("--max_shift");
     maxAngularChange = getDoubleParam("--max_angular_change");
     maxResol = getDoubleParam("--max_resolution");
     Ts = getDoubleParam("--sampling");
     Rmax = getIntParam("--Rmax");
     RmaxDef = getIntParam("--RDef");
     optimizeAlignment = checkParam("--optimizeAlignment");
     optimizeDeformation = checkParam("--optimizeDeformation");
     optimizeDefocus = checkParam("--optimizeDefocus");
     phaseFlipped = checkParam("--phaseFlipped");
     L1 = getIntParam("--l1");
     L2 = getIntParam("--l2");
     lambda = getDoubleParam("--regularization");
     resume = checkParam("--resume");
     Rerunable::setFileName(fnOutDir + "/sphDone.xmd");
 }

◆ show()

void ProgAngularSphAlignment::show ( )

Show.

Definition at line 70 of file angular_sph_alignment.cpp.

 {
     if (!verbose)
         return;
     XmippMetadataProgram::show();
     std::cout
     << "Output directory:    " << fnOutDir           << std::endl
     << "Reference volume:    " << fnVolR             << std::endl
     << "Reference mask:      " << fnMaskR            << std::endl
     << "Max. Shift:          " << maxShift           << std::endl
     << "Max. Angular Change: " << maxAngularChange   << std::endl
     << "Max. Resolution:     " << maxResol           << std::endl
     << "Sampling:            " << Ts                 << std::endl
     << "Max. Radius:         " << Rmax               << std::endl
     << "Max. Radius Deform.  " << RmaxDef            << std::endl
     << "Zernike Degree:      " << L1                 << std::endl
     << "SH Degree:           " << L2                 << std::endl
     << "Optimize alignment:  " << optimizeAlignment  << std::endl
     << "Optimize deformation:" << optimizeDeformation<< std::endl
     << "Optimize defocus;    " << optimizeDefocus    << std::endl
     << "Phase flipped:       " << phaseFlipped       << std::endl
     << "Regularization:      " << lambda             << std::endl
     ;
 }

◆ tranformImageSph()

double ProgAngularSphAlignment::tranformImageSph	(	double *	pclnm,
		double	rot,
		double	tilt,
		double	psi,
		double	deltaDefocusU,
		double	deltaDefocusV,
		double	deltaDefocusAngle
	)

Definition at line 196 of file angular_sph_alignment.cpp.

 {
     const MultidimArray<double> &mV=V();
     FOR_ALL_ELEMENTS_IN_MATRIX1D(clnm)
         VEC_ELEM(clnm,i)=pclnm[i+1];
     double deformation=0.0;
     totalDeformation=0.0;
     P().initZeros((int)XSIZE(I()),(int)XSIZE(I()));
     P().setXmippOrigin();
     deformVol(P(), mV, deformation, rot, tilt, psi);
     if (hasCTF) {
         applyCTFImage(deltaDefocusU, deltaDefocusV, deltaDefocusAngle);
     }
     double cost=0;
     if (old_flip)
     {
         MAT_ELEM(A,0,0)*=-1;
         MAT_ELEM(A,0,1)*=-1;
         MAT_ELEM(A,0,2)*=-1;
     }
 
     applyGeometry(xmipp_transformation::LINEAR,Ifilteredp(),Ifiltered(),A,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP,0.);
     filter.applyMaskSpace(P());
     const MultidimArray<double> mP=P();
     const MultidimArray<int> &mMask2D=mask2D;
     MultidimArray<double> &mIfilteredp=Ifilteredp();
     double corr=correlationIndex(mIfilteredp,mP,&mMask2D);
     cost=-corr;
 
 #ifdef DEBUG
     std::cout << "A=" << A << std::endl;
     Image<double> save;
     save()=P();
     save.write("PPPtheo.xmp");
     save()=Ifilteredp();
     save.write("PPPfilteredp.xmp");
     save()=Ifiltered();
     save.write("PPPfiltered.xmp");
     // Vdeformed.write("PPPVdeformed.vol");
     std::cout << "Cost=" << cost << " deformation=" << deformation << std::endl;
     std::cout << "Press any key" << std::endl;
     char c; std::cin >> c;
 #endif
 
    if (showOptimization)
    {
         std::cout << "A=" << A << std::endl;
         Image<double> save(P);
         save.write("PPPtheo.xmp");
         save()=Ifilteredp();
         save.write("PPPfilteredp.xmp");
         save()=Ifiltered();
         save.write("PPPfiltered.xmp");
         std::cout << "Cost=" << cost << " corr=" << corr << std::endl;
         std::cout << "Deformation=" << totalDeformation << std::endl;
         std::cout << "Press any key" << std::endl;
         char c; std::cin >> c;
     }
 
     double massDiff=std::abs(sumV-sumVd)/sumV;
     double retval=cost+lambda*(deformation + massDiff);
     if (showOptimization)
         std::cout << cost << " " << deformation << " " << lambda*deformation << " " << sumV << " " << sumVd << " " << massDiff << " " << retval << std::endl;
     return retval;
 }

◆ updateCTFImage()

void ProgAngularSphAlignment::updateCTFImage	(	double	defocusU,
		double	defocusV,
		double	angle
	)

Definition at line 521 of file angular_sph_alignment.cpp.

 {
     ctf.K=1; // get pure CTF with no envelope
     currentDefocusU=ctf.DeltafU=defocusU;
     currentDefocusV=ctf.DeltafV=defocusV;
     currentAngle=ctf.azimuthal_angle=angle;
     ctf.produceSideInfo();
 }

◆ writeImageParameters()

void ProgAngularSphAlignment::writeImageParameters ( const FileName & fnImg )

virtual

Write the parameters found for one image

Reimplemented in MpiProgAngularSphAlignment.

Definition at line 424 of file angular_sph_alignment.cpp.

                                                                         {
     MetaDataVec md;
     int pos = 3*vecSize;
     size_t objId = md.addObject();
     md.setValue(MDL_IMAGE, fnImg, objId);
     if (flagEnabled==1) {
         md.setValue(MDL_ENABLED, 1, objId);
     }
     else {
         md.setValue(MDL_ENABLED, -1, objId);
     }
     md.setValue(MDL_ANGLE_ROT,   old_rot+p(pos+2), objId);
     md.setValue(MDL_ANGLE_TILT,  old_tilt+p(pos+3), objId);
     md.setValue(MDL_ANGLE_PSI,   old_psi+p(pos+4), objId);
     md.setValue(MDL_SHIFT_X,     old_shiftX+p(pos+0), objId);
     md.setValue(MDL_SHIFT_Y,     old_shiftY+p(pos+1), objId);
     md.setValue(MDL_FLIP,        old_flip, objId);
     md.setValue(MDL_SPH_DEFORMATION, totalDeformation, objId);
     std::vector<double> vectortemp;
     for (int j = 0; j < VEC_XSIZE(clnm); j++) {
         vectortemp.push_back(clnm(j));
     }
     md.setValue(MDL_SPH_COEFFICIENTS, vectortemp, objId);
     md.setValue(MDL_COST,        correlation, objId);
     md.append(Rerunable::getFileName());
 }

Member Data Documentation

◆ A

Matrix2D<double> ProgAngularSphAlignment::A

Definition at line 109 of file angular_sph_alignment.h.

◆ clnm

Matrix1D<double> ProgAngularSphAlignment::clnm

Definition at line 136 of file angular_sph_alignment.h.

◆ correlation

double ProgAngularSphAlignment::correlation

Definition at line 146 of file angular_sph_alignment.h.

◆ ctf

CTFDescription ProgAngularSphAlignment::ctf

Definition at line 130 of file angular_sph_alignment.h.

◆ currentAngle

double ProgAngularSphAlignment::currentAngle

Definition at line 128 of file angular_sph_alignment.h.

◆ currentDefocusU

double ProgAngularSphAlignment::currentDefocusU

Definition at line 126 of file angular_sph_alignment.h.

◆ currentDefocusV

double ProgAngularSphAlignment::currentDefocusV

Definition at line 127 of file angular_sph_alignment.h.

◆ filter

FourierFilter ProgAngularSphAlignment::filter

Definition at line 107 of file angular_sph_alignment.h.

◆ FilterCTF

FourierFilter ProgAngularSphAlignment::FilterCTF

Definition at line 132 of file angular_sph_alignment.h.

◆ flagEnabled

int ProgAngularSphAlignment::flagEnabled

Definition at line 87 of file angular_sph_alignment.h.

◆ fnMaskR

FileName ProgAngularSphAlignment::fnMaskR

Filename of the reference volume mask

Definition at line 48 of file angular_sph_alignment.h.

◆ fnOutDir

FileName ProgAngularSphAlignment::fnOutDir

Output directory.

Definition at line 50 of file angular_sph_alignment.h.

◆ fnVolR

FileName ProgAngularSphAlignment::fnVolR

Filename of the reference volume

Definition at line 46 of file angular_sph_alignment.h.

◆ hasCTF

bool ProgAngularSphAlignment::hasCTF

Definition at line 120 of file angular_sph_alignment.h.

◆ I

Image<double> ProgAngularSphAlignment::I

Definition at line 100 of file angular_sph_alignment.h.

◆ Ifiltered

Image<double> ProgAngularSphAlignment::Ifiltered

Definition at line 102 of file angular_sph_alignment.h.

◆ Ifilteredp

Image<double> ProgAngularSphAlignment::Ifilteredp

Definition at line 103 of file angular_sph_alignment.h.

◆ ignoreCTF

bool ProgAngularSphAlignment::ignoreCTF

Definition at line 76 of file angular_sph_alignment.h.

◆ Ip

Image<double> ProgAngularSphAlignment::Ip

Definition at line 101 of file angular_sph_alignment.h.

◆ L1

int ProgAngularSphAlignment::L1

Degrees of Zernike polynomials and spherical harmonics

Definition at line 52 of file angular_sph_alignment.h.

◆ L2

int ProgAngularSphAlignment::L2

Definition at line 53 of file angular_sph_alignment.h.

◆ lambda

double ProgAngularSphAlignment::lambda

Definition at line 82 of file angular_sph_alignment.h.

◆ mask2D

MultidimArray<int> ProgAngularSphAlignment::mask2D

Definition at line 93 of file angular_sph_alignment.h.

◆ maxAngularChange

double ProgAngularSphAlignment::maxAngularChange

Maximum angular change allowed

Definition at line 62 of file angular_sph_alignment.h.

◆ maxResol

double ProgAngularSphAlignment::maxResol

Maximum frequency (A)

Definition at line 64 of file angular_sph_alignment.h.

◆ maxShift

double ProgAngularSphAlignment::maxShift

Maximum shift allowed

Definition at line 60 of file angular_sph_alignment.h.

◆ old_defocusAngle

double ProgAngularSphAlignment::old_defocusAngle

Definition at line 124 of file angular_sph_alignment.h.

◆ old_defocusU

double ProgAngularSphAlignment::old_defocusU

Definition at line 122 of file angular_sph_alignment.h.

◆ old_defocusV

double ProgAngularSphAlignment::old_defocusV

Definition at line 123 of file angular_sph_alignment.h.

◆ old_flip

bool ProgAngularSphAlignment::old_flip

Definition at line 118 of file angular_sph_alignment.h.

◆ old_psi

double ProgAngularSphAlignment::old_psi

Definition at line 113 of file angular_sph_alignment.h.

◆ old_rot

double ProgAngularSphAlignment::old_rot

Definition at line 111 of file angular_sph_alignment.h.

◆ old_shiftX

double ProgAngularSphAlignment::old_shiftX

Definition at line 115 of file angular_sph_alignment.h.

◆ old_shiftY

double ProgAngularSphAlignment::old_shiftY

Definition at line 116 of file angular_sph_alignment.h.

◆ old_tilt

double ProgAngularSphAlignment::old_tilt

Definition at line 112 of file angular_sph_alignment.h.

◆ optimizeAlignment

bool ProgAngularSphAlignment::optimizeAlignment

Definition at line 70 of file angular_sph_alignment.h.

◆ optimizeDefocus

bool ProgAngularSphAlignment::optimizeDefocus

Definition at line 74 of file angular_sph_alignment.h.

◆ optimizeDeformation

bool ProgAngularSphAlignment::optimizeDeformation

Definition at line 72 of file angular_sph_alignment.h.

◆ optimizeRadius

bool ProgAngularSphAlignment::optimizeRadius

Definition at line 78 of file angular_sph_alignment.h.

◆ p

Matrix1D<double> ProgAngularSphAlignment::p

Definition at line 86 of file angular_sph_alignment.h.

◆ P

Image<double> ProgAngularSphAlignment::P

Definition at line 105 of file angular_sph_alignment.h.

◆ phaseFlipped

bool ProgAngularSphAlignment::phaseFlipped

Definition at line 80 of file angular_sph_alignment.h.

◆ resume

bool ProgAngularSphAlignment::resume

Resume computations

Definition at line 91 of file angular_sph_alignment.h.

◆ Rmax

int ProgAngularSphAlignment::Rmax

Maximum radius

Definition at line 68 of file angular_sph_alignment.h.

◆ RmaxDef

int ProgAngularSphAlignment::RmaxDef

Definition at line 84 of file angular_sph_alignment.h.

◆ showOptimization

bool ProgAngularSphAlignment::showOptimization

Definition at line 144 of file angular_sph_alignment.h.

◆ steps_cp

Matrix1D<double> ProgAngularSphAlignment::steps_cp

Definition at line 138 of file angular_sph_alignment.h.

◆ sumV

double ProgAngularSphAlignment::sumV

Definition at line 141 of file angular_sph_alignment.h.

◆ sumVd

double ProgAngularSphAlignment::sumVd

Definition at line 142 of file angular_sph_alignment.h.

◆ totalDeformation

double ProgAngularSphAlignment::totalDeformation

Definition at line 140 of file angular_sph_alignment.h.

◆ Ts

double ProgAngularSphAlignment::Ts

Sampling rate

Definition at line 66 of file angular_sph_alignment.h.

◆ V

Image<double> ProgAngularSphAlignment::V

Definition at line 98 of file angular_sph_alignment.h.

◆ V_mask

MultidimArray<int> ProgAngularSphAlignment::V_mask

Definition at line 94 of file angular_sph_alignment.h.

◆ Vdeformed

Image<double> ProgAngularSphAlignment::Vdeformed

Definition at line 99 of file angular_sph_alignment.h.

◆ vecSize

int ProgAngularSphAlignment::vecSize

Definition at line 134 of file angular_sph_alignment.h.

◆ vL1

Matrix1D<int> ProgAngularSphAlignment::vL1

Zernike and SPH coefficients vectors

Definition at line 55 of file angular_sph_alignment.h.

◆ vL2

Matrix1D<int> ProgAngularSphAlignment::vL2

Definition at line 57 of file angular_sph_alignment.h.

◆ vM

Matrix1D<int> ProgAngularSphAlignment::vM

Definition at line 58 of file angular_sph_alignment.h.

◆ vN

Matrix1D<int> ProgAngularSphAlignment::vN

Definition at line 56 of file angular_sph_alignment.h.

◆ Xdim

size_t ProgAngularSphAlignment::Xdim

Definition at line 96 of file angular_sph_alignment.h.

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

xmipp/libraries/reconstruction/angular_sph_alignment.h
xmipp/libraries/reconstruction/angular_sph_alignment.cpp

Public Member Functions

Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgAngularSphAlignment()

◆ ~ProgAngularSphAlignment()

Member Function Documentation

◆ applyCTFImage()

◆ createWorkFiles()

◆ defineParams()

◆ deformVol()

◆ fillVectorTerms()

◆ finishProcessing()

◆ minimizepos()

◆ numCoefficients()

◆ preProcess()

◆ processImage()

◆ readParams()

◆ show()

◆ tranformImageSph()

◆ updateCTFImage()

◆ writeImageParameters()

Member Data Documentation

◆ A

◆ clnm

◆ correlation

◆ ctf

◆ currentAngle

◆ currentDefocusU

◆ currentDefocusV

◆ filter

◆ FilterCTF

◆ flagEnabled

◆ fnMaskR

◆ fnOutDir

◆ fnVolR

◆ hasCTF

◆ I

◆ Ifiltered

◆ Ifilteredp

◆ ignoreCTF

◆ Ip

◆ L1

◆ L2

◆ lambda

◆ mask2D

◆ maxAngularChange

◆ maxResol

◆ maxShift

◆ old_defocusAngle

◆ old_defocusU

◆ old_defocusV

◆ old_flip

◆ old_psi

◆ old_rot

◆ old_shiftX

◆ old_shiftY

◆ old_tilt

◆ optimizeAlignment

◆ optimizeDefocus

◆ optimizeDeformation

◆ optimizeRadius

◆ p

◆ P

◆ phaseFlipped

◆ resume

◆ Rmax

◆ RmaxDef

◆ showOptimization

◆ steps_cp

◆ sumV

◆ sumVd

◆ totalDeformation

◆ Ts

◆ V

◆ V_mask

◆ Vdeformed

◆ vecSize