#include <forward_zernike_subtomos.h>

Inheritance diagram for ProgForwardZernikeSubtomos:

Collaboration diagram for ProgForwardZernikeSubtomos:

Public Types
enum	Direction { Direction::ROTATE, Direction::UNROTATE }

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

	~ProgForwardZernikeSubtomos ()
	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 &vecSize)
	Length of coefficients vector. More...

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

void	fillVectorTerms (int l1, int l2, Matrix1D< int > &vL1, Matrix1D< int > &vN, Matrix1D< int > &vL2, Matrix1D< int > &vM)
	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	updateCTFImage (double defocusU, double defocusV, double angle)

double	transformImageSph (double *pclnm)

template<Direction DIRECTION>
void	rotateCoefficients ()

virtual void	finishProcessing ()

virtual void	writeImageParameters (MDRow &row)

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

Matrix1D< double >	weightsInterpolation3D (double x, double y, double z)

void	splattingAtPos (std::array< double, 3 > r, double weight, MultidimArray< double > &mP, const MultidimArray< double > &mV)

void	rotatePositions (double rot, double tilt, double psi)

void	preComputeDF ()

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

int	num_images

int	algn_params

int	ctf_params

Matrix1D< double >	p

int	flagEnabled

bool	useCTF

double	t1

double	t2

bool	resume

MultidimArray< int >	mask3D

MultidimArray< int >	V_mask

size_t	Xdim

FileName	fnImage

Image< double >	V

Image< double >	Vdeformed

Image< double >	I

Image< double >	Ifiltered

Image< double >	Ifilteredp

Image< double >	P

FourierFilter	filter

FourierFilter	filterMW

Matrix2D< double >	A

double	old_rot

double	old_tilt

double	old_psi

double	deltaRot

double	deltaTilt

double	deltaPsi

double	old_shiftX

double	old_shiftY

double	old_shiftZ

double	deltaX

double	deltaY

double	deltaZ

bool	old_flip

bool	hasCTF

double	old_defocusU

double	old_defocusV

double	old_defocusAngle

double	deltaDefocusU

double	deltaDefocusV

double	deltaDefocusAngle

double	currentDefocusU

double	currentDefocusV

double	currentAngle

FourierFilter	FilterCTF

int	vecSize

Matrix1D< double >	clnm

Matrix1D< double >	prev_clnm

Matrix1D< double >	steps_cp

double	totalDeformation

double	prior_deformation

int	sumV

bool	showOptimization

double	correlation

int	loop_step

struct blobtype	blob

double	blob_r

double	sigma4

Matrix1D< double >	gaussianProjectionTable

Matrix1D< double >	gaussianProjectionTable2

MultidimArray< double >	vpos

MultidimArray< double >	df

std::vector< size_t >	idx_z_clnm

std::vector< double >	z_clnm_diff

Matrix2D< double >	R

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
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 43 of file forward_zernike_subtomos.h.

Member Enumeration Documentation

◆ Direction

enum ProgForwardZernikeSubtomos::Direction

strong

Enumerator
ROTATE
UNROTATE

Definition at line 162 of file forward_zernike_subtomos.h.

162 { ROTATE, UNROTATE };

Constructor & Destructor Documentation

◆ ProgForwardZernikeSubtomos()

ProgForwardZernikeSubtomos::ProgForwardZernikeSubtomos ( )

Empty constructor.

Definition at line 35 of file forward_zernike_subtomos.cpp.

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

◆ ~ProgForwardZernikeSubtomos()

ProgForwardZernikeSubtomos::~ProgForwardZernikeSubtomos ( )

default

Destructor.

Member Function Documentation

◆ checkPoint()

void ProgForwardZernikeSubtomos::checkPoint ( )

virtual

For very long programs, it may be needed to write checkpoints.

Reimplemented from XmippMetadataProgram.

Reimplemented in MpiProgForwardZernikeSubtomos.

Definition at line 622 of file forward_zernike_subtomos.cpp.

                                             {
     MDRowVec rowAppend;
     MetaDataVec checkPoint;
     getOutputMd().getRow(rowAppend, getOutputMd().lastRowId());
     checkPoint.addRow(rowAppend);
     checkPoint.append(Rerunable::getFileName());
 }

◆ createWorkFiles()

void ProgForwardZernikeSubtomos::createWorkFiles ( )

inlineprotected

Definition at line 231 of file forward_zernike_subtomos.h.

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

Rerunable::createWorkFiles

virtual void createWorkFiles(bool resume, MetaData *md)

Definition: rerunable_program.h:49

ProgForwardZernikeSubtomos::resume

bool resume

Definition: forward_zernike_subtomos.h:96

XmippMetadataProgram::getInputMd

MetaData * getInputMd()

Definition: xmipp_metadata_program.h:58

◆ defineParams()

void ProgForwardZernikeSubtomos::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippMetadataProgram.

Definition at line 105 of file forward_zernike_subtomos.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");
     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("  [--step <step=1>]            : Voxel index step");
     addParamsLine("  [--useCTF]                   : Correct CTF");
     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("  [--blobr <b=4>]              : Blob radius for forward mapping splatting");
     addParamsLine("  [--t1 <t1=-60>]              : First tilt angle range");
     addParamsLine("  [--t2 <t2=60>]               : Second tilt angle range");
     addParamsLine("  [--resume]                   : Resume processing");
     addExampleLine("A typical use is:",false);
     addExampleLine("xmipp_forward_zernike_subtomos -i anglesFromContinuousAssignment.xmd --ref reference.vol -o assigned_anglesAndDeformations.xmd --optimizeAlignment --optimizeDeformation --depth 1");
 }

◆ deformVol()

void ProgForwardZernikeSubtomos::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 729 of file forward_zernike_subtomos.cpp.

 {
     size_t idxY0=(VEC_XSIZE(clnm)-9)/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_inv;
     // Matrix2D<double> R;
     // R.initIdentity(3);
     // Euler_angles2matrix(rot, tilt, psi, R, false);
     // Matrix2D<double> R_inv = R.inv();
     // R_inv = R.inv();
 
     // auto stepsMask = std::vector<size_t>();
     // if (optimizeDeformation)
     // {
     //  for (size_t idx = 0; idx < idxY0; idx++)
     //  {
     //      if (1 == VEC_ELEM(steps_cp, idx))
     //      {
     //          stepsMask.emplace_back(idx);
     //      }
     //  }
     // }
     // else {
     //  for (size_t idx = 0; idx < idxY0; idx++)
     //  {
     //      stepsMask.emplace_back(idx);
     //  }
     // }
 
     auto sz = idx_z_clnm.size();
     Matrix1D<int> l1, l2, n, m, idx_v;
 
     if (!idx_z_clnm.empty())
     {
         l1.initZeros(sz);
         l2.initZeros(sz);
         n.initZeros(sz);
         m.initZeros(sz);
         idx_v.initZeros(sz);
         for (auto j=0; j<sz; j++)
         {
             auto idx = idx_z_clnm[j];
             if (idx >= idxY0 && idx < idxZ0)
                 idx -= idxY0;
             else if (idx >= idxZ0)
                 idx -= idxZ0;
 
             VEC_ELEM(idx_v,j) = idx;
             VEC_ELEM(l1,j) = VEC_ELEM(vL1, idx);
             VEC_ELEM(n,j) = VEC_ELEM(vN, idx);
             VEC_ELEM(l2,j) = VEC_ELEM(vL2, idx);
             VEC_ELEM(m,j) = VEC_ELEM(vM, idx);
         }
     }
 
     // // TODO: Poner primero i y j en el loop, acumular suma y guardar al final
     // const auto lastZ = FINISHINGZ(mV);
     // const auto lastY = FINISHINGY(mV);
     // const auto lastX = FINISHINGX(mV);
     // for (int k=STARTINGZ(mV); k<=lastZ; k+=loop_step)
     // {
     //  for (int i=STARTINGY(mV); i<=lastY; i+=loop_step)
     //  {
     //      for (int j=STARTINGX(mV); j<=lastX; j+=loop_step)
     //      {
     //          if (A3D_ELEM(V_mask,k,i,j) == 1) {
     //              // ZZ(p) = k; YY(p) = i; XX(p) = j;
     //              // pos = R_inv * pos;
     //              // pos = R * pos;
     //              double gx=0.0, gy=0.0, gz=0.0;
     //              double k2=k*k;
     //              double kr=k*iRmaxF;
     //              double k2i2=k2+i*i;
     //              double ir=i*iRmaxF;
     //              double r2=k2i2+j*j;
     //              double jr=j*iRmaxF;
     //              double rr=sqrt(r2)*iRmaxF;
     //              for (auto idx : stepsMask) {
     //                  auto l1 = VEC_ELEM(vL1,idx);
     //                  auto n = VEC_ELEM(vN,idx);
     //                  auto l2 = VEC_ELEM(vL2,idx);
     //                  auto m = VEC_ELEM(vM,idx);
     //                  auto zsph=ZernikeSphericalHarmonics(l1,n,l2,m,jr,ir,kr,rr);
     //                  auto c = std::array<double, 3>{};
     //                  // XX(c_rot) = VEC_ELEM(clnm,idx); YY(c_rot) = VEC_ELEM(clnm,idx+idxY0); ZZ(c_rot) = VEC_ELEM(clnm,idx+idxZ0);
     //                  // if (num_images == 1 && optimizeDeformation)
     //                  // {
     //                  //  //? Hacer algun check para ahorrarnos cuentas si no usamos priors (en ese 
     //                  //  //? caso podemos usar las lineas comentadas)
     //                  //  double c_x = VEC_ELEM(clnm,idx);
     //                  //  double c_y = VEC_ELEM(clnm,idx+idxY0);
     //                  //  // c[0] = R_inv.mdata[0] * c_x + R_inv.mdata[1] * c_y;
     //                  //  // c[1] = R_inv.mdata[3] * c_x + R_inv.mdata[4] * c_y;
     //                  //  // c[2] = R_inv.mdata[6] * c_x + R_inv.mdata[7] * c_y;
     //                  //  double c_z = VEC_ELEM(clnm, idx + idxZ0);
     //                  //  c[0] = R_inv.mdata[0] * c_x + R_inv.mdata[1] * c_y + R_inv.mdata[2] * c_z;
     //                  //  c[1] = R_inv.mdata[3] * c_x + R_inv.mdata[4] * c_y + R_inv.mdata[5] * c_z;
     //                  //  c[2] = R_inv.mdata[6] * c_x + R_inv.mdata[7] * c_y + R_inv.mdata[8] * c_z;
     //                  // }
     //                  // else {
     //                  c[0] = VEC_ELEM(clnm,idx);
     //                  c[1] = VEC_ELEM(clnm,idx+idxY0);
     //                  c[2] = VEC_ELEM(clnm,idx+idxZ0);
     //                  // }
     //                  if (rr>0 || l2==0) {
     //                      gx += c[0]  *(zsph);
     //                      gy += c[1]  *(zsph);
     //                      gz += c[2]  *(zsph);
     //                  }
     //              }
 
     //              auto pos = std::array<double, 3>{};
     //              double r_x = j + gx;
     //              double r_y = i + gy;
     //              double r_z = k + gz;
     //              pos[0] = R.mdata[0] * r_x + R.mdata[1] * r_y + R.mdata[2] * r_z;
     //              pos[1] = R.mdata[3] * r_x + R.mdata[4] * r_y + R.mdata[5] * r_z;
     //              pos[2] = R.mdata[6] * r_x + R.mdata[7] * r_y + R.mdata[8] * r_z;
                     
     //              double voxel_mV = A3D_ELEM(mV,k,i,j);
     //              splattingAtPos(pos, voxel_mV, mP, mV);
     //              if (!mV.outside(pos[2], pos[1], pos[0]))
     //                  sumVd += voxel_mV;
     //              modg += gx*gx+gy*gy+gz*gz;
     //              Ncount++;
     //          }
     //      }
     //  }
     // }
 
     const auto &mVpos = vpos;
     const auto lastY = FINISHINGY(mVpos);
     for (int i=STARTINGY(mVpos); i<=lastY; i++)
     {
         double &gx = A2D_ELEM(df, i, 0);
         double &gy = A2D_ELEM(df, i, 1);
         double &gz = A2D_ELEM(df, i, 2);
         double r_x = A2D_ELEM(mVpos, i, 0);
         double r_y = A2D_ELEM(mVpos, i, 1);
         double r_z = A2D_ELEM(mVpos, i, 2);
         double xr = A2D_ELEM(mVpos, i, 3);
         double yr = A2D_ELEM(mVpos, i, 4);
         double zr = A2D_ELEM(mVpos, i, 5);
         double rr = A2D_ELEM(mVpos, i, 6);
 
         if (!idx_z_clnm.empty())
         {
             for (auto j = 0; j < sz; j++)
             {   
                 auto idx = VEC_ELEM(idx_v, j);
                 auto aux_l2 = VEC_ELEM(l2, j);
                 auto zsph = ZernikeSphericalHarmonics(VEC_ELEM(l1, j), VEC_ELEM(n, j),
                                                       aux_l2, VEC_ELEM(m, j), xr, yr, zr, rr);
 
                 auto diff_c_x = VEC_ELEM(clnm, idx) - VEC_ELEM(prev_clnm, idx);
                 auto diff_c_y = VEC_ELEM(clnm, idx + idxY0) - VEC_ELEM(prev_clnm, idx + idxY0);
                 auto diff_c_z = VEC_ELEM(clnm, idx + idxZ0) - VEC_ELEM(prev_clnm, idx + idxZ0);
                 // auto i_diff_c_x = R_inv.mdata[0] * diff_c_x + R_inv.mdata[1] * diff_c_y;
                 // auto i_diff_c_y = R_inv.mdata[3] * diff_c_x + R_inv.mdata[4] * diff_c_y;
                 // auto i_diff_c_z = R_inv.mdata[6] * diff_c_x + R_inv.mdata[7] * diff_c_y;
                 if (rr > 0 || aux_l2 == 0)
                 {
                     gx += diff_c_x * (zsph);
                     gy += diff_c_y * (zsph);
                     gz += diff_c_z * (zsph);
                 }
             }
         }
 
         auto r_gx = R.mdata[0] * gx + R.mdata[1] * gy + R.mdata[2] * gz;
         auto r_gy = R.mdata[3] * gx + R.mdata[4] * gy + R.mdata[5] * gz;
         auto r_gz = R.mdata[6] * gx + R.mdata[7] * gy + R.mdata[8] * gz;
 
         auto pos = std::array<double, 3>{};
         pos[0] = r_x + r_gx;
         pos[1] = r_y + r_gy;
         pos[2] = r_z + r_gz;
 
         double voxel_mV = A2D_ELEM(mVpos, i, 7);
         splattingAtPos(pos, voxel_mV, mP, mV);
 
         modg += gx * gx + gy * gy + gz * gz;
         Ncount++;
     }
 
     // def=sqrt(modg/Ncount);
     def = sqrt(modg/Ncount);
     totalDeformation = def;
 }

◆ fillVectorTerms()

void ProgForwardZernikeSubtomos::fillVectorTerms	(	int	l1,
		int	l2,
		Matrix1D< int > &	vL1,
		Matrix1D< int > &	vN,
		Matrix1D< int > &	vL2,
		Matrix1D< int > &	vM
	)

Zernike and SPH coefficients allocation.

Definition at line 673 of file forward_zernike_subtomos.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;
         int aux = 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 ProgForwardZernikeSubtomos::finishProcessing ( )

virtual

Write the final parameters.

Reimplemented from XmippMetadataProgram.

Reimplemented in MpiProgForwardZernikeSubtomos.

Definition at line 249 of file forward_zernike_subtomos.cpp.

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

◆ minimizepos()

void ProgForwardZernikeSubtomos::minimizepos	(	int	L1,
		int	l2,
		Matrix1D< double > &	steps
	)

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

Definition at line 646 of file forward_zernike_subtomos.cpp.

 {
     int size = 0;
     int prevSize = 0;
     numCoefficients(L1,l2,size);
     numCoefficients(L1,l2-1,prevSize);
     int totalSize = (steps.size()-9)/3;
     if (l2 > 1)
     {
         for (int idx = prevSize; idx < size; idx++)
         {
             VEC_ELEM(steps, idx) = 1.;
             VEC_ELEM(steps, idx + totalSize) = 1.;
             VEC_ELEM(steps, idx + 2 * totalSize) = 1.;
         }
     }
     else
     {
         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 ProgForwardZernikeSubtomos::numCoefficients	(	int	l1,
		int	l2,
		int &	vecSize
	)

Length of coefficients vector.

Definition at line 630 of file forward_zernike_subtomos.cpp.

 {
     for (int h=0; h<=l2; h++)
     {
         int numSPH = 2*h+1;
         int count=l1-h+1;
         int numEven=(count>>1)+(count&1 && !(h&1));
         if (h%2 == 0) {
             vecSize += numSPH*numEven;
         }
         else {
             vecSize += numSPH*(l1-h+1-numEven);
         }
     }
 }

◆ preComputeDF()

void ProgForwardZernikeSubtomos::preComputeDF ( )

Definition at line 1072 of file forward_zernike_subtomos.cpp.

 {   
     size_t idxY0=(VEC_XSIZE(clnm)-9)/3;
     size_t idxZ0=2*idxY0;
     // Matrix2D<double> R_inv = R.inv();
     const auto &mVpos = vpos;
     const auto lastY = FINISHINGY(mVpos);
     for (int i=STARTINGY(mVpos); i<=lastY; i++)
     {
         double &gx = A2D_ELEM(df, i, 0);
         double &gy = A2D_ELEM(df, i, 1);
         double &gz = A2D_ELEM(df, i, 2);
         double r_x = A2D_ELEM(mVpos, i, 0);
         double r_y = A2D_ELEM(mVpos, i, 1);
         double r_z = A2D_ELEM(mVpos, i, 2);
         double xr = A2D_ELEM(mVpos, i, 3);
         double yr = A2D_ELEM(mVpos, i, 4);
         double zr = A2D_ELEM(mVpos, i, 5);
         double rr = A2D_ELEM(mVpos, i, 6);
 
         if (!idx_z_clnm.empty())
         {
             for (int idx = 0; idx < idxY0; idx++)
             {
                 auto aux_l2 = VEC_ELEM(vL2, idx);
                 auto zsph = ZernikeSphericalHarmonics(VEC_ELEM(vL1, idx), VEC_ELEM(vN, idx),
                                                       aux_l2, VEC_ELEM(vM, idx), xr, yr, zr, rr);
                 auto c_x = VEC_ELEM(clnm, idx);
                 auto c_y = VEC_ELEM(clnm, idx + idxY0);
                 auto c_z = VEC_ELEM(clnm, idx + idxZ0);
                 // auto i_c_x = R_inv.mdata[0] * c_x + R_inv.mdata[1] * c_y + R_inv.mdata[2] * c_z;
                 // auto i_c_y = R_inv.mdata[3] * c_x + R_inv.mdata[4] * c_y + R_inv.mdata[5] * c_z;
                 // auto i_c_z = R_inv.mdata[6] * c_x + R_inv.mdata[7] * c_y + R_inv.mdata[8] * c_z;
                 if (rr > 0 || aux_l2 == 0)
                 {
                     gx += c_x * (zsph);
                     gy += c_y * (zsph);
                     gz += c_z * (zsph);
                 }
             }
         }
     }   
 }

◆ preProcess()

void ProgForwardZernikeSubtomos::preProcess ( )

virtual

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

Reimplemented from XmippMetadataProgram.

Definition at line 139 of file forward_zernike_subtomos.cpp.

 {
     V.read(fnVolR);
     V().setXmippOrigin();
     Xdim=XSIZE(V());
     Vdeformed().initZeros(V());
     Vdeformed().setXmippOrigin();
     // sumV=V().sum();
 
     Ifilteredp().initZeros(Xdim,Xdim,Xdim);
     Ifilteredp().setXmippOrigin();
     P().initZeros(Xdim,Xdim,Xdim);
 
     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();
         double Rmax2 = RmaxDef*RmaxDef;
         for (int k=STARTINGZ(V_mask); k<=FINISHINGZ(V_mask); k++) {
             for (int i=STARTINGY(V_mask); i<=FINISHINGY(V_mask); i++) {
                 for (int j=STARTINGX(V_mask); j<=FINISHINGX(V_mask); j++) {
                     double r2 = k*k + i*i + j*j;
                     if (r2>=Rmax2)
                         A3D_ELEM(V_mask,k,i,j) = 0;
                 }
             }
         }
     }
     else {
         mask.R1 = RmaxDef;
         mask.generate_mask(V());
         V_mask = mask.get_binary_mask();
         V_mask.setXmippOrigin();
     }
 
     // Total Volume Mass (Inside Mask)
     sumV = 0;
     for (int k = STARTINGZ(V()); k <= FINISHINGZ(V()); k += loop_step)
     {
         for (int i = STARTINGY(V()); i <= FINISHINGY(V()); i += loop_step)
         {
             for (int j = STARTINGX(V()); j <= FINISHINGX(V()); j += loop_step)
             {
                 if (A3D_ELEM(V_mask, k, i, j) == 1)
                 {
                     sumV++;
                 }
             }
         }
     }
 
     // Construct mask
     if (Rmax<0)
         Rmax=Xdim/2;
     mask.R1 = Rmax;
     mask.generate_mask(Xdim,Xdim,Xdim);
     mask3D=mask.get_binary_mask();
 
     // Low pass filter
     filter.FilterBand=LOWPASS;
     filter.do_generate_3dmask = true;
     filter.w1=Ts/maxResol;
     filter.raised_w=0.02;
 
     // MW filter
     // filterMW.FilterBand=WEDGE;
     // filterMW.FilterShape=WEDGE;
     // filterMW.do_generate_3dmask = true;
     // filterMW.t1 = t1;
     // filterMW.t2 = t2;
     // filterMW.rot=filterMW.tilt=filterMW.psi=0.0;
     filterMW.FilterBand=LOWPASS;
     filterMW.FilterShape=WEDGE_RC;
     filterMW.do_generate_3dmask = true;
     filterMW.t1 = t1;
     filterMW.t2 = t2;
     filterMW.rot=filterMW.tilt=filterMW.psi=0.0;
     filterMW.w1=Ts/maxResol;
     filterMW.raised_w=0.02;
 
     // Transformation matrix
     A.initIdentity(4);
 
     // CTF Filter
     FilterCTF.FilterBand = CTF;
     FilterCTF.do_generate_3dmask = true;
     FilterCTF.ctf.enable_CTFnoise = false;
     FilterCTF.ctf.produceSideInfo();
 
     vecSize = 0;
     numCoefficients(L1,L2,vecSize);
     fillVectorTerms(L1,L2,vL1,vN,vL2,vM);
 
     createWorkFiles();
 
     // Blob
     blob.radius = blob_r;   // Blob radius in voxels
     blob.order  = 2;        // Order of the Bessel function
     blob.alpha  = 7.05;     // Smoothness parameter
 
 }

◆ processImage()

void ProgForwardZernikeSubtomos::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 416 of file forward_zernike_subtomos.cpp.

 {
     Matrix1D<double> steps;
     // totalSize = 3*num_Z_coeff + 3 shifts + 3 angles + 3 CTF
     int totalSize = 3*vecSize + 9;
     p.initZeros(totalSize);
     clnm.initZeros(totalSize);
     prev_clnm.initZeros(totalSize);
 
     // Init positions and deformation field
     vpos.initZeros(sumV, 8);
     df.initZeros(sumV, 3);
 
     flagEnabled=1;
 
     rowIn.getValueOrDefault(MDL_IMAGE,       fnImage, "");
     rowIn.getValueOrDefault(MDL_ANGLE_ROT,   old_rot, 0.0);
     rowIn.getValueOrDefault(MDL_ANGLE_TILT,  old_tilt, 0.0);
     rowIn.getValueOrDefault(MDL_ANGLE_PSI,   old_psi, 0.0);
     rowIn.getValueOrDefault(MDL_SHIFT_X,     old_shiftX, 0.0);
     rowIn.getValueOrDefault(MDL_SHIFT_Y,     old_shiftY, 0.0);
     rowIn.getValueOrDefault(MDL_SHIFT_Z,     old_shiftZ, 0.0);
     I.read(fnImage); 
     I().setXmippOrigin();
     Ifiltered() = I(); 
     filter.applyMaskSpace(Ifiltered());
     rotatePositions(old_rot, old_tilt, old_psi);
     
     if (verbose >= 2)
         std::cout << "Processing Image (" << fnImage << ")" << std::endl;
 
     if (rowIn.containsLabel(MDL_FLIP))
         rowIn.getValue(MDL_FLIP,old_flip);
     else
         old_flip = false;
 
     prior_deformation = 0.0;
     if (rowIn.containsLabel(MDL_SPH_COEFFICIENTS))
     {
         std::vector<double> vectortemp;
         rowIn.getValue(MDL_SPH_COEFFICIENTS, vectortemp);
         rowIn.getValueOrDefault(MDL_SPH_DEFORMATION, prior_deformation, 0.0);
         for (int i=0; i<3*vecSize; i++)
         {
             clnm[i] = vectortemp[i];
         }
         // if (optimizeDeformation)
         //  rotateCoefficients<Direction::ROTATE>();
         p = clnm;
         prev_clnm = clnm;
         preComputeDF();
     }   
     
     // FIXME: Add defocus per image and make CTF correction available
     if ((rowIn.containsLabel(MDL_CTF_DEFOCUSU) || rowIn.containsLabel(MDL_CTF_MODEL)) && useCTF)
     {
         hasCTF=true;
         FilterCTF.ctf.readFromMdRow(rowIn);
         FilterCTF.ctf.Tm = Ts;
         FilterCTF.ctf.produceSideInfo();
         old_defocusU=FilterCTF.ctf.DeltafU;
         old_defocusV=FilterCTF.ctf.DeltafV;
         old_defocusAngle=FilterCTF.ctf.azimuthal_angle;
     }
     else
         hasCTF=false;
 
     // If deformation is not optimized, do a single iteration
     //? Si usamos priors es mejor ir poco a poco, ir poco a poco pero usar todos los coeffs cada vez (mas lento)
     //? O dar solo una iteracion con todos los coeffs?
     int h = 1;
     // if (!optimizeDeformation)
     // if (rowIn.containsLabel(MDL_SPH_COEFFICIENTS))
     //  h = L2;
 
     for (;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;
         correlation = 0.0;
         try
         {
             cost=1e38;
             int iter;
             if (optimizeAlignment)
             {
                 int init = steps.size() - 9;
                 int end = steps.size() - 3;
                 for (int i = init; i < end; i++)
                     steps(i)=1.;
             }
             if (optimizeDefocus) 
             {
                 int init = steps.size() - 3;
                 int end = steps.size();
                 for (int i = init; i < end; i++)
                     steps(i)=1.;
             }
             if (optimizeDeformation)
             {
                 minimizepos(L1,h,steps);
             }
             steps_cp = steps;
             powellOptimizer(p, 1, totalSize, &continuousZernikeSubtomoCost, this, 0.1, cost, iter, steps, verbose>=2);
 
             if (verbose>=3)
             {
                 showOptimization = true;
                 continuousZernikeSubtomoCost(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;
         }
     }
 
     // if (optimizeDeformation)
     // {
     //  rotateCoefficients<Direction::UNROTATE>();
     // }
 
     //AJ NEW
     writeImageParameters(rowOut);
     //END AJ
 
 }

◆ readParams()

void ProgForwardZernikeSubtomos::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippMetadataProgram.

Definition at line 46 of file forward_zernike_subtomos.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");
     useCTF = checkParam("--useCTF");
     L1 = getIntParam("--l1");
     L2 = getIntParam("--l2");
     loop_step = getIntParam("--step");
     lambda = getDoubleParam("--regularization");
     resume = checkParam("--resume");
     Rerunable::setFileName(fnOutDir + "/sphDone.xmd");
     blob_r = getDoubleParam("--blobr");
     t1 = getDoubleParam("--t1");
     t2 = getDoubleParam("--t2");
     keep_input_columns = true;
 }

◆ rotateCoefficients()

template<ProgForwardZernikeSubtomos::Direction DIRECTION>

void ProgForwardZernikeSubtomos::rotateCoefficients ( )

Definition at line 706 of file forward_zernike_subtomos.cpp.

                                                     {
     int pos = 3*vecSize;
     size_t idxY0=(VEC_XSIZE(clnm)-9)/3;
     size_t idxZ0=2*idxY0;
 
     double rot = old_rot+p(pos+2);
     double tilt = old_tilt+p(pos+3);
     double psi = old_psi+p(pos+4);
 
     Matrix2D<double> R;
     R.initIdentity(3);
     Euler_angles2matrix(rot, tilt, psi, R, false);
     if (DIRECTION == Direction::UNROTATE)
         R = R.inv();
     Matrix1D<double> c;
     c.initZeros(3);
     for (size_t idx=0; idx<idxY0; idx++) {
         XX(c) = VEC_ELEM(clnm,idx); YY(c) = VEC_ELEM(clnm,idx+idxY0); ZZ(c) = VEC_ELEM(clnm,idx+idxZ0);
         c = R * c;
         VEC_ELEM(clnm,idx) = XX(c); VEC_ELEM(clnm,idx+idxY0) = YY(c); VEC_ELEM(clnm,idx+idxZ0) = ZZ(c);
     }
 }

◆ rotatePositions()

void ProgForwardZernikeSubtomos::rotatePositions	(	double	rot,
		double	tilt,
		double	psi
	)

Definition at line 1028 of file forward_zernike_subtomos.cpp.

 {
     // Matrix2D<double> R;
     R.initIdentity(3);
     Euler_angles2matrix(rot, tilt, psi, R, false);
 
     const MultidimArray<double> &mV=V();
 
     const auto lastZ = FINISHINGZ(mV);
     const auto lastY = FINISHINGY(mV);
     const auto lastX = FINISHINGX(mV);
     size_t count = 0;
     double iRmaxF=1.0/RmaxDef;
     for (int k=STARTINGZ(mV); k<=lastZ; k+=loop_step)
     {
         for (int i=STARTINGY(mV); i<=lastY; i+=loop_step)
         {
             for (int j=STARTINGX(mV); j<=lastX; j+=loop_step)
             {
                 if (A3D_ELEM(V_mask,k,i,j) == 1) {
                     double x = j;
                     double y = i;
                     double z = k;
                     double r_x = R.mdata[0] * x + R.mdata[1] * y + R.mdata[2] * z;
                     double r_y = R.mdata[3] * x + R.mdata[4] * y + R.mdata[5] * z;
                     double r_z = R.mdata[6] * x + R.mdata[7] * y + R.mdata[8] * z;
 
                     A2D_ELEM(vpos, count, 0) = r_x;
                     A2D_ELEM(vpos, count, 1) = r_y;
                     A2D_ELEM(vpos, count, 2) = r_z;
                     A2D_ELEM(vpos, count, 3) = j * iRmaxF;
                     A2D_ELEM(vpos, count, 4) = i * iRmaxF;
                     A2D_ELEM(vpos, count, 5) = z * iRmaxF;
                     A2D_ELEM(vpos, count, 6) = sqrt(x*x + y*y + z*z) * iRmaxF;
                     A2D_ELEM(vpos, count, 7) = A3D_ELEM(mV, k, i, j);
 
                     count++;
                 }
             }
         }
     }
 } 

◆ show()

void ProgForwardZernikeSubtomos::show ( )

Show.

Definition at line 76 of file forward_zernike_subtomos.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
     << "Step:                " << loop_step          << std::endl
     << "Correct CTF:         " << useCTF             << 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
     << "Blob radius:         " << blob_r             << std::endl
     ;
 }

◆ splattingAtPos()

void ProgForwardZernikeSubtomos::splattingAtPos	(	std::array< double, 3 >	r,
		double	weight,
		MultidimArray< double > &	mP,
		const MultidimArray< double > &	mV
	)

Definition at line 961 of file forward_zernike_subtomos.cpp.

                                                                                                                                                 {
     // Find the part of the volume that must be updated
     // double x_pos = r[0];
     // double y_pos = r[1];
     // double z_pos = r[2];
     // // int k0 = XMIPP_MAX(FLOOR(z_pos - blob_r), STARTINGZ(mV));
     // // int kF = XMIPP_MIN(CEIL(z_pos + blob_r), FINISHINGZ(mV));
     // // int i0 = XMIPP_MAX(FLOOR(y_pos - blob_r), STARTINGY(mV));
     // // int iF = XMIPP_MIN(CEIL(y_pos + blob_r), FINISHINGY(mV));
     // // int j0 = XMIPP_MAX(FLOOR(x_pos - blob_r), STARTINGX(mV));
     // // int jF = XMIPP_MIN(CEIL(x_pos + blob_r), FINISHINGX(mV));
     // int k0 = XMIPP_MAX(FLOOR(z_pos - sigma4), STARTINGZ(mV));
     // int kF = XMIPP_MIN(CEIL(z_pos + sigma4), FINISHINGZ(mV));
     // int i0 = XMIPP_MAX(FLOOR(y_pos - sigma4), STARTINGY(mV));
     // int iF = XMIPP_MIN(CEIL(y_pos + sigma4), FINISHINGY(mV));
     // int j0 = XMIPP_MAX(FLOOR(x_pos - sigma4), STARTINGX(mV));
     // int jF = XMIPP_MIN(CEIL(x_pos + sigma4), FINISHINGX(mV));
     // int size = gaussianProjectionTable.size();
     // for (int k = k0; k <= kF; k++)
     // {
     //  double k2 = (z_pos - k) * (z_pos - k);
     //  for (int i = i0; i <= iF; i++)
     //  {
     //      double y2 = (y_pos - i) * (y_pos - i);
     //      for (int j = j0; j <= jF; j++)
     //      {
     //          // double mod = sqrt((x_pos - j) * (x_pos - j) + y2 + k2);
     //          // // A3D_ELEM(Vdeformed(),k, i, j) += weight * blob_val(rdiff.module(), blob);
     //          // A3D_ELEM(mP, k, i, j) += weight * kaiser_value(mod, blob.radius, blob.alpha, blob.order);
     //          double mod = sqrt((x_pos - j) * (x_pos - j) + y2 + k2);
     //          double didx = mod * 1000;
     //          int idx = ROUND(didx);
     //          if (idx < size)
     //          {
     //              double gw = gaussianProjectionTable.vdata[idx];
     //              A3D_ELEM(mP, k, i, j) += weight * gw;
     //          }
     //      }
     //  }
     // }
 
     double x_pos = r[0];
     double y_pos = r[1];
     double z_pos = r[2];
     int i0 = XMIPP_MAX(CEIL(y_pos - loop_step), STARTINGY(mV));
     int iF = XMIPP_MIN(FLOOR(y_pos + loop_step), FINISHINGY(mV));
     int j0 = XMIPP_MAX(CEIL(x_pos - loop_step), STARTINGX(mV));
     int jF = XMIPP_MIN(FLOOR(x_pos + loop_step), FINISHINGX(mV));
     int k0 = XMIPP_MAX(CEIL(z_pos - loop_step), STARTINGZ(mV));
     int kF = XMIPP_MIN(FLOOR(z_pos + loop_step), FINISHINGZ(mV));
 
     double m = 1. / loop_step;
     for (int k = k0; k <= kF; k++)
     {
         double z_val = 1. - m * ABS(k - z_pos);
         for (int i = i0; i <= iF; i++)
         {
             double y_val = 1. - m * ABS(i - y_pos);
             for (int j = j0; j <= jF; j++)
             {
                 double x_val = 1. - m * ABS(j - x_pos);
                 A3D_ELEM(mP, k, i, j) += weight * x_val * y_val * z_val;
             }
         }
     }
 }

◆ transformImageSph()

double ProgForwardZernikeSubtomos::transformImageSph ( double * pclnm )

Definition at line 255 of file forward_zernike_subtomos.cpp.

 {
     const MultidimArray<double> &mV=V();
     idx_z_clnm.clear();
     z_clnm_diff.clear();
     FOR_ALL_ELEMENTS_IN_MATRIX1D(clnm)
     {
         VEC_ELEM(clnm,i)=pclnm[i+1];
         if (VEC_ELEM(clnm,i) != VEC_ELEM(prev_clnm,i))
         {
             idx_z_clnm.push_back(i);
             z_clnm_diff.push_back(VEC_ELEM(clnm, i) - VEC_ELEM(prev_clnm, i));
             // std::cout << i << std::endl;
         }
     }
     double deformation=0.0;
     totalDeformation=0.0;
 
     P().initZeros(Xdim,Xdim,Xdim);
     P().setXmippOrigin();
     double currentRot=old_rot + deltaRot;
     double currentTilt=old_tilt + deltaTilt;
     double currentPsi=old_psi + deltaPsi;
     deformVol(P(), mV, deformation, currentRot, currentTilt, currentPsi);
 
     double cost=0.0;
     const MultidimArray<int> &mMask3D=mask3D;
 
     if (old_flip)
     {
         MAT_ELEM(A, 0, 0) *= -1;
         MAT_ELEM(A, 0, 1) *= -1;
         MAT_ELEM(A, 0, 2) *= -1;
     }
 
     // Image<double> save;
     // save()=Ifiltered();
     // save.write("PPPbefore.mrc");
 
     applyGeometry(xmipp_transformation::LINEAR,Ifilteredp(),Ifiltered(),A,
                     xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP,0.);
 
     // Image<double> save;
     // save()=P();
     // save.write("PPPtheo.mrc");
 
     // save()=Ifilteredp();
     // save.write("PPPexp.mrc");
 
     // std::cout << "Press any key" << std::endl;
     // char c; std::cin >> c;
 
     filterMW.generateMask(P());
     filterMW.applyMaskSpace(P());
 
     if (hasCTF)
     {
         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);
         }
         // FilterCTF1.ctf = ctf;
         FilterCTF.generateMask(P());
         if (phaseFlipped)
             FilterCTF.correctPhase();
         FilterCTF.applyMaskSpace(P());
     }
     // filter.applyMaskSpace(P());
 
 
     const MultidimArray<double> &mP=P();
     MultidimArray<double> &mIfilteredp=Ifilteredp();
     double corr1=correlationIndex(mIfilteredp,mP,&mMask3D);
 
     cost=-corr1;
 
 #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;
         save()=P();
         save.write("PPPtheo.mrc");
         save()=Ifilteredp();
         save.write("PPPfilteredp.mrc");
         save()=Ifiltered();
         save.write("PPPfiltered.mrc");
         Vdeformed.write("PPPVdeformed.vol");
         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;
     prev_clnm = clnm;
     double retval=cost+lambda*abs(deformation - prior_deformation);
     if (showOptimization)
         std::cout << cost << " " << deformation << " " << lambda*deformation << " " << sumV << " " << retval << std::endl;
     return retval;
 }

◆ updateCTFImage()

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

Definition at line 696 of file forward_zernike_subtomos.cpp.

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

◆ weightsInterpolation3D()

Matrix1D< double > ProgForwardZernikeSubtomos::weightsInterpolation3D	(	double	x,
		double	y,
		double	z
	)

Definition at line 930 of file forward_zernike_subtomos.cpp.

                                                                                                 {
     Matrix1D<double> w;
     w.initZeros(8);
 
     int x0 = FLOOR(x);
     double fx0 = x - x0;
     int x1 = x0 + 1;
     double fx1 = x1 - x;
 
     int y0 = FLOOR(y);
     double fy0 = y - y0;
     int y1 = y0 + 1;
     double fy1 = y1 - y;
 
     int z0 = FLOOR(z);
     double fz0 = z - z0;
     int z1 = z0 + 1;
     double fz1 = z1 - z;
 
     VEC_ELEM(w,0) = fx1 * fy1 * fz1;  // w000 (x0,y0,z0)
     VEC_ELEM(w,1) = fx1 * fy1 * fz0;  // w001 (x0,y0,z1)
     VEC_ELEM(w,2) = fx1 * fy0 * fz1;  // w010 (x0,y1,z0)
     VEC_ELEM(w,3) = fx1 * fy0 * fz0;  // w011 (x0,y1,z1)
     VEC_ELEM(w,4) = fx0 * fy1 * fz1;  // w100 (x1,y0,z0)
     VEC_ELEM(w,5) = fx0 * fy1 * fz0;  // w101 (x1,y0,z1)
     VEC_ELEM(w,6) = fx0 * fy0 * fz1;  // w110 (x1,y1,z0)
     VEC_ELEM(w,7) = fx0 * fy0 * fz0;  // w111 (x1,y1,z1)
 
     return w;
 }

◆ writeImageParameters()

void ProgForwardZernikeSubtomos::writeImageParameters ( MDRow & row )

virtual

Write the parameters found for one image

Definition at line 596 of file forward_zernike_subtomos.cpp.

                                                                 {
     int pos = 3*vecSize;
     if (flagEnabled==1) {
         row.setValue(MDL_ENABLED, 1);
     }
     else {
         row.setValue(MDL_ENABLED, -1);
     }
     
     row.setValue(MDL_ANGLE_ROT,   old_rot+p(pos+3));
     row.setValue(MDL_ANGLE_TILT,  old_tilt+p(pos+4));
     row.setValue(MDL_ANGLE_PSI,   old_psi+p(pos+5));
     row.setValue(MDL_SHIFT_X,     old_shiftX+p(pos));
     row.setValue(MDL_SHIFT_Y,     old_shiftY+p(pos+1));
     row.setValue(MDL_SHIFT_Z,     old_shiftZ+p(pos+2));
 
     row.setValue(MDL_SPH_DEFORMATION, totalDeformation);
     std::vector<double> vectortemp;
     size_t end_clnm = VEC_XSIZE(clnm)-9;
     for (int j = 0; j < end_clnm; j++) {
         vectortemp.push_back(clnm(j));
     }
     row.setValue(MDL_SPH_COEFFICIENTS, vectortemp);
     // row.setValue(MDL_COST, correlation);
 }

Member Data Documentation

◆ A

Matrix2D<double> ProgForwardZernikeSubtomos::A

Definition at line 114 of file forward_zernike_subtomos.h.

◆ algn_params

int ProgForwardZernikeSubtomos::algn_params

Definition at line 85 of file forward_zernike_subtomos.h.

◆ blob

struct blobtype ProgForwardZernikeSubtomos::blob

Definition at line 145 of file forward_zernike_subtomos.h.

◆ blob_r

double ProgForwardZernikeSubtomos::blob_r

Definition at line 146 of file forward_zernike_subtomos.h.

◆ clnm

Matrix1D<double> ProgForwardZernikeSubtomos::clnm

Definition at line 132 of file forward_zernike_subtomos.h.

◆ correlation

double ProgForwardZernikeSubtomos::correlation

Definition at line 141 of file forward_zernike_subtomos.h.

◆ ctf_params

int ProgForwardZernikeSubtomos::ctf_params

Definition at line 87 of file forward_zernike_subtomos.h.

◆ currentAngle

double ProgForwardZernikeSubtomos::currentAngle

Definition at line 126 of file forward_zernike_subtomos.h.

◆ currentDefocusU

double ProgForwardZernikeSubtomos::currentDefocusU

Definition at line 126 of file forward_zernike_subtomos.h.

◆ currentDefocusV

double ProgForwardZernikeSubtomos::currentDefocusV

Definition at line 126 of file forward_zernike_subtomos.h.

◆ deltaDefocusAngle

double ProgForwardZernikeSubtomos::deltaDefocusAngle

Definition at line 124 of file forward_zernike_subtomos.h.

◆ deltaDefocusU

double ProgForwardZernikeSubtomos::deltaDefocusU

Definition at line 124 of file forward_zernike_subtomos.h.

◆ deltaDefocusV

double ProgForwardZernikeSubtomos::deltaDefocusV

Definition at line 124 of file forward_zernike_subtomos.h.

◆ deltaPsi

double ProgForwardZernikeSubtomos::deltaPsi

Definition at line 116 of file forward_zernike_subtomos.h.

◆ deltaRot

double ProgForwardZernikeSubtomos::deltaRot

Definition at line 116 of file forward_zernike_subtomos.h.

◆ deltaTilt

double ProgForwardZernikeSubtomos::deltaTilt

Definition at line 116 of file forward_zernike_subtomos.h.

◆ deltaX

double ProgForwardZernikeSubtomos::deltaX

Definition at line 118 of file forward_zernike_subtomos.h.

◆ deltaY

double ProgForwardZernikeSubtomos::deltaY

Definition at line 118 of file forward_zernike_subtomos.h.

◆ deltaZ

double ProgForwardZernikeSubtomos::deltaZ

Definition at line 118 of file forward_zernike_subtomos.h.

◆ df

MultidimArray<double> ProgForwardZernikeSubtomos::df

Definition at line 156 of file forward_zernike_subtomos.h.

◆ filter

FourierFilter ProgForwardZernikeSubtomos::filter

Definition at line 111 of file forward_zernike_subtomos.h.

◆ FilterCTF

FourierFilter ProgForwardZernikeSubtomos::FilterCTF

Definition at line 128 of file forward_zernike_subtomos.h.

◆ filterMW

FourierFilter ProgForwardZernikeSubtomos::filterMW

Definition at line 112 of file forward_zernike_subtomos.h.

◆ flagEnabled

int ProgForwardZernikeSubtomos::flagEnabled

Definition at line 89 of file forward_zernike_subtomos.h.

◆ fnImage

FileName ProgForwardZernikeSubtomos::fnImage

Definition at line 102 of file forward_zernike_subtomos.h.

◆ fnMaskR

FileName ProgForwardZernikeSubtomos::fnMaskR

Filename of the reference volume mask

Definition at line 49 of file forward_zernike_subtomos.h.

◆ fnOutDir

FileName ProgForwardZernikeSubtomos::fnOutDir

Output directory.

Definition at line 51 of file forward_zernike_subtomos.h.

◆ fnVolR

FileName ProgForwardZernikeSubtomos::fnVolR

Filename of the reference volume

Definition at line 47 of file forward_zernike_subtomos.h.

◆ gaussianProjectionTable

Matrix1D<double> ProgForwardZernikeSubtomos::gaussianProjectionTable

Definition at line 150 of file forward_zernike_subtomos.h.

◆ gaussianProjectionTable2

Matrix1D<double> ProgForwardZernikeSubtomos::gaussianProjectionTable2

Definition at line 153 of file forward_zernike_subtomos.h.

◆ hasCTF

bool ProgForwardZernikeSubtomos::hasCTF

Definition at line 122 of file forward_zernike_subtomos.h.

◆ I

Image<double> ProgForwardZernikeSubtomos::I

Definition at line 105 of file forward_zernike_subtomos.h.

◆ idx_z_clnm

std::vector<size_t> ProgForwardZernikeSubtomos::idx_z_clnm

Definition at line 157 of file forward_zernike_subtomos.h.

◆ Ifiltered

Image<double> ProgForwardZernikeSubtomos::Ifiltered

Definition at line 106 of file forward_zernike_subtomos.h.

◆ Ifilteredp

Image<double> ProgForwardZernikeSubtomos::Ifilteredp

Definition at line 107 of file forward_zernike_subtomos.h.

◆ ignoreCTF

bool ProgForwardZernikeSubtomos::ignoreCTF

Definition at line 73 of file forward_zernike_subtomos.h.

◆ L1

int ProgForwardZernikeSubtomos::L1

Degrees of Zernike polynomials and spherical harmonics

Definition at line 53 of file forward_zernike_subtomos.h.

◆ L2

int ProgForwardZernikeSubtomos::L2

Definition at line 53 of file forward_zernike_subtomos.h.

◆ lambda

double ProgForwardZernikeSubtomos::lambda

Definition at line 79 of file forward_zernike_subtomos.h.

◆ loop_step

int ProgForwardZernikeSubtomos::loop_step

Definition at line 143 of file forward_zernike_subtomos.h.

◆ mask3D

MultidimArray<int> ProgForwardZernikeSubtomos::mask3D

Definition at line 98 of file forward_zernike_subtomos.h.

◆ maxAngularChange

double ProgForwardZernikeSubtomos::maxAngularChange

Maximum angular change allowed

Definition at line 59 of file forward_zernike_subtomos.h.

◆ maxResol

double ProgForwardZernikeSubtomos::maxResol

Maximum frequency (A)

Definition at line 61 of file forward_zernike_subtomos.h.

◆ maxShift

double ProgForwardZernikeSubtomos::maxShift

Maximum shift allowed

Definition at line 57 of file forward_zernike_subtomos.h.

◆ num_images

int ProgForwardZernikeSubtomos::num_images

Definition at line 83 of file forward_zernike_subtomos.h.

◆ old_defocusAngle

double ProgForwardZernikeSubtomos::old_defocusAngle

Definition at line 124 of file forward_zernike_subtomos.h.

◆ old_defocusU

double ProgForwardZernikeSubtomos::old_defocusU

Definition at line 124 of file forward_zernike_subtomos.h.

◆ old_defocusV

double ProgForwardZernikeSubtomos::old_defocusV

Definition at line 124 of file forward_zernike_subtomos.h.

◆ old_flip

bool ProgForwardZernikeSubtomos::old_flip

Definition at line 120 of file forward_zernike_subtomos.h.

◆ old_psi

double ProgForwardZernikeSubtomos::old_psi

Definition at line 116 of file forward_zernike_subtomos.h.

◆ old_rot

double ProgForwardZernikeSubtomos::old_rot

Definition at line 116 of file forward_zernike_subtomos.h.

◆ old_shiftX

double ProgForwardZernikeSubtomos::old_shiftX

Definition at line 118 of file forward_zernike_subtomos.h.

◆ old_shiftY

double ProgForwardZernikeSubtomos::old_shiftY

Definition at line 118 of file forward_zernike_subtomos.h.

◆ old_shiftZ

double ProgForwardZernikeSubtomos::old_shiftZ

Definition at line 118 of file forward_zernike_subtomos.h.

◆ old_tilt

double ProgForwardZernikeSubtomos::old_tilt

Definition at line 116 of file forward_zernike_subtomos.h.

◆ optimizeAlignment

bool ProgForwardZernikeSubtomos::optimizeAlignment

Definition at line 67 of file forward_zernike_subtomos.h.

◆ optimizeDefocus

bool ProgForwardZernikeSubtomos::optimizeDefocus

Definition at line 71 of file forward_zernike_subtomos.h.

◆ optimizeDeformation

bool ProgForwardZernikeSubtomos::optimizeDeformation

Definition at line 69 of file forward_zernike_subtomos.h.

◆ optimizeRadius

bool ProgForwardZernikeSubtomos::optimizeRadius

Definition at line 75 of file forward_zernike_subtomos.h.

◆ p

Matrix1D<double> ProgForwardZernikeSubtomos::p

Definition at line 88 of file forward_zernike_subtomos.h.

◆ P

Image<double> ProgForwardZernikeSubtomos::P

Definition at line 109 of file forward_zernike_subtomos.h.

◆ phaseFlipped

bool ProgForwardZernikeSubtomos::phaseFlipped

Definition at line 77 of file forward_zernike_subtomos.h.

◆ prev_clnm

Matrix1D<double> ProgForwardZernikeSubtomos::prev_clnm

Definition at line 132 of file forward_zernike_subtomos.h.

◆ prior_deformation

double ProgForwardZernikeSubtomos::prior_deformation

Definition at line 136 of file forward_zernike_subtomos.h.

◆ R

Matrix2D<double> ProgForwardZernikeSubtomos::R

Definition at line 159 of file forward_zernike_subtomos.h.

◆ resume

bool ProgForwardZernikeSubtomos::resume

Resume computations

Definition at line 96 of file forward_zernike_subtomos.h.

◆ Rmax

int ProgForwardZernikeSubtomos::Rmax

Maximum radius

Definition at line 65 of file forward_zernike_subtomos.h.

◆ RmaxDef

int ProgForwardZernikeSubtomos::RmaxDef

Definition at line 81 of file forward_zernike_subtomos.h.

◆ showOptimization

bool ProgForwardZernikeSubtomos::showOptimization

Definition at line 139 of file forward_zernike_subtomos.h.

◆ sigma4

double ProgForwardZernikeSubtomos::sigma4

Definition at line 148 of file forward_zernike_subtomos.h.

◆ steps_cp

Matrix1D<double> ProgForwardZernikeSubtomos::steps_cp

Definition at line 134 of file forward_zernike_subtomos.h.

◆ sumV

int ProgForwardZernikeSubtomos::sumV

Definition at line 137 of file forward_zernike_subtomos.h.

◆ t1

double ProgForwardZernikeSubtomos::t1

Definition at line 92 of file forward_zernike_subtomos.h.

◆ t2

double ProgForwardZernikeSubtomos::t2

Definition at line 92 of file forward_zernike_subtomos.h.

◆ totalDeformation

double ProgForwardZernikeSubtomos::totalDeformation

Definition at line 136 of file forward_zernike_subtomos.h.

◆ Ts

double ProgForwardZernikeSubtomos::Ts

Sampling rate

Definition at line 63 of file forward_zernike_subtomos.h.

◆ useCTF

bool ProgForwardZernikeSubtomos::useCTF

Definition at line 90 of file forward_zernike_subtomos.h.

◆ V

Image<double> ProgForwardZernikeSubtomos::V

Definition at line 104 of file forward_zernike_subtomos.h.

◆ V_mask

MultidimArray<int> ProgForwardZernikeSubtomos::V_mask

Definition at line 98 of file forward_zernike_subtomos.h.

◆ Vdeformed

Image<double> ProgForwardZernikeSubtomos::Vdeformed

Definition at line 104 of file forward_zernike_subtomos.h.

◆ vecSize

int ProgForwardZernikeSubtomos::vecSize

Definition at line 130 of file forward_zernike_subtomos.h.

◆ vL1

Matrix1D<int> ProgForwardZernikeSubtomos::vL1

Zernike and SPH coefficients vectors

Definition at line 55 of file forward_zernike_subtomos.h.

◆ vL2

Matrix1D<int> ProgForwardZernikeSubtomos::vL2

Definition at line 55 of file forward_zernike_subtomos.h.

◆ vM

Matrix1D<int> ProgForwardZernikeSubtomos::vM

Definition at line 55 of file forward_zernike_subtomos.h.

◆ vN

Matrix1D<int> ProgForwardZernikeSubtomos::vN

Definition at line 55 of file forward_zernike_subtomos.h.

◆ vpos

MultidimArray<double> ProgForwardZernikeSubtomos::vpos

Definition at line 156 of file forward_zernike_subtomos.h.

◆ Xdim

size_t ProgForwardZernikeSubtomos::Xdim

Definition at line 100 of file forward_zernike_subtomos.h.

◆ z_clnm_diff

std::vector<double> ProgForwardZernikeSubtomos::z_clnm_diff

Definition at line 158 of file forward_zernike_subtomos.h.

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

xmipp/libraries/reconstruction/forward_zernike_subtomos.h
xmipp/libraries/reconstruction/forward_zernike_subtomos.cpp

Public Types

Public Member Functions

Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Member Enumeration Documentation

◆ Direction

Constructor & Destructor Documentation

◆ ProgForwardZernikeSubtomos()

◆ ~ProgForwardZernikeSubtomos()

Member Function Documentation

◆ checkPoint()

◆ createWorkFiles()

◆ defineParams()

◆ deformVol()

◆ fillVectorTerms()

◆ finishProcessing()

◆ minimizepos()

◆ numCoefficients()

◆ preComputeDF()

◆ preProcess()

◆ processImage()

◆ readParams()

◆ rotateCoefficients()

◆ rotatePositions()

◆ show()

◆ splattingAtPos()

◆ transformImageSph()

◆ updateCTFImage()

◆ weightsInterpolation3D()

◆ writeImageParameters()

Member Data Documentation

◆ A

◆ algn_params

◆ blob

◆ blob_r

◆ clnm

◆ correlation

◆ ctf_params

◆ currentAngle

◆ currentDefocusU

◆ currentDefocusV

◆ deltaDefocusAngle

◆ deltaDefocusU

◆ deltaDefocusV

◆ deltaPsi

◆ deltaRot

◆ deltaTilt

◆ deltaX

◆ deltaY

◆ deltaZ

◆ df

◆ filter

◆ FilterCTF

◆ filterMW

◆ flagEnabled

◆ fnImage

◆ fnMaskR

◆ fnOutDir

◆ fnVolR

◆ gaussianProjectionTable

◆ gaussianProjectionTable2

◆ hasCTF

◆ I

◆ idx_z_clnm

◆ Ifiltered

◆ Ifilteredp

◆ ignoreCTF

◆ L1

◆ L2

◆ lambda

◆ loop_step

◆ mask3D

◆ maxAngularChange

◆ maxResol

◆ maxShift

◆ num_images

◆ old_defocusAngle

◆ old_defocusU