Xmipp  v3.23.11-Nereus
Public Types | Public Member Functions | Public Attributes | Protected Member Functions | List of all members
ProgForwardZernikeImages Class Reference
Libraries » Reconstruction library » angular_continuous_assign2 (Continuous angular assignment)

#include <forward_zernike_images.h>

Inheritance diagram for ProgForwardZernikeImages:
Inheritance graph
[legend]
Collaboration diagram for ProgForwardZernikeImages:
Collaboration graph
[legend]

Public Types

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

Public Member Functions

 ProgForwardZernikeImages ()
 Empty constructor. More...
 
 ~ProgForwardZernikeImages ()
 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)
 
double bspline1 (double x)
 
double bspline3 (double x)
 
void removePixels ()
 
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
 
int image_mode
 
bool useCTF
 
bool resume
 
MultidimArray< int > mask2D
 
MultidimArray< int > V_mask
 
size_t Xdim
 
std::vector< FileName > fnImage
 
Image< double > V
 
Image< double > Vdeformed
 
std::vector< Image< double > > I
 
std::vector< Image< double > > Ifiltered
 
std::vector< Image< double > > Ifilteredp
 
std::vector< Image< double > > P
 
FourierFilter filter
 
FourierFilter filterp
 
Matrix2D< double > A1
 
Matrix2D< double > A2
 
Matrix2D< double > A3
 
std::vector< double > old_rot
 
std::vector< double > old_tilt
 
std::vector< double > old_psi
 
std::vector< double > deltaRot
 
std::vector< double > deltaTilt
 
std::vector< double > deltaPsi
 
std::vector< double > old_shiftX
 
std::vector< double > old_shiftY
 
std::vector< double > deltaX
 
std::vector< double > deltaY
 
bool old_flip
 
bool hasCTF
 
std::vector< double > old_defocusU
 
std::vector< double > old_defocusV
 
std::vector< double > old_defocusAngle
 
std::vector< double > deltaDefocusU
 
std::vector< double > deltaDefocusV
 
std::vector< double > deltaDefocusAngle
 
std::vector< double > currentDefocusU
 
std::vector< double > currentDefocusV
 
std::vector< double > currentAngle
 
FourierFilter FilterCTF1
 
FourierFilter FilterCTF2
 
FourierFilter FilterCTF3
 
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
 
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_images.h.

Member Enumeration Documentation

◆ Direction

enum ProgForwardZernikeImages::Direction
strong
Enumerator
ROTATE 
UNROTATE 

Definition at line 155 of file forward_zernike_images.h.

155 { ROTATE, UNROTATE };

Constructor & Destructor Documentation

◆ ProgForwardZernikeImages()

ProgForwardZernikeImages::ProgForwardZernikeImages ( )

Empty constructor.

Definition at line 35 of file forward_zernike_images.cpp.

35  : Rerunable("")
36 {
37  resume = false;
38  produces_a_metadata = true;
39  each_image_produces_an_output = false;
40  showOptimization = false;
41 }
XmippMetadataProgram::each_image_produces_an_output
bool each_image_produces_an_output
Indicate that an output is produced for each image in the input.
Definition: xmipp_metadata_program.h:90
Rerunable::Rerunable
Rerunable(const FileName &fn)
Definition: rerunable_program.h:43
XmippMetadataProgram::produces_a_metadata
bool produces_a_metadata
Indicate that the unique final output file is a Metadata.
Definition: xmipp_metadata_program.h:88

◆ ~ProgForwardZernikeImages()

ProgForwardZernikeImages::~ProgForwardZernikeImages ( )
default

Destructor.

Member Function Documentation

◆ bspline1()

double ProgForwardZernikeImages::bspline1 ( double  x)

Definition at line 1194 of file forward_zernike_images.cpp.

1195 {
1196  double m = 1. / blob_r;
1197  if (0. < x && x < blob_r)
1198  return m * (blob_r - x);
1199  else if (-blob_r < x && x <= 0.)
1200  return m * (blob_r + x);
1201  else
1202  return 0.;
1203 }
x
doublereal * x
Definition: numerical_recipes.cpp:2230
m
int m
Definition: image_find_center.cpp:57

◆ bspline3()

double ProgForwardZernikeImages::bspline3 ( double  x)

Definition at line 1205 of file forward_zernike_images.cpp.

1206 {
1207  double c = 1. / 6.;
1208  double sx = x / blob_r;
1209  double sx2 = sx*sx;
1210  double sx3 = sx2*sx;
1211  double double_int = 2. * blob_r;
1212  if (-double_int < x && x < -blob_r)
1213  return c * (sx3 + 6.*sx2 + 12.*sx + 8.);
1214  else if (-blob_r <= x && x< 0.)
1215  return c * (-3.*sx3 - 6.*sx2 + 4.);
1216  else if (0. <= x && x < blob_r)
1217  return c * (3.*sx3 - 6.*sx2 + 4.);
1218  else if (blob_r <= x && x < double_int)
1219  return c * (-sx3 + 6.*sx2 - 12.*sx + 8.);
1220  else
1221  return 0.;
1222 }
c
doublereal * c
Definition: numerical_recipes.cpp:2230
x
doublereal * x
Definition: numerical_recipes.cpp:2230

◆ checkPoint()

void ProgForwardZernikeImages::checkPoint ( )
virtual

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

Reimplemented from XmippMetadataProgram.

Reimplemented in MpiProgForwardZernikeImages.

Definition at line 934 of file forward_zernike_images.cpp.

934  {
935  MDRowVec rowAppend;
936  MetaDataVec checkPoint;
937  getOutputMd().getRow(rowAppend, getOutputMd().lastRowId());
938  checkPoint.addRow(rowAppend);
939  checkPoint.append(Rerunable::getFileName());
940 }
MetaDataVec::getRow
std::unique_ptr< MDRow > getRow(size_t id) override
Definition: metadata_vec.cpp:310
MetaDataVec::addRow
size_t addRow(const MDRow &row) override
Definition: metadata_vec.cpp:225
Rerunable::getFileName
const FileName & getFileName() const
Definition: rerunable_program.h:80
MetaData::append
void append(const FileName &outFile) const
Definition: metadata_base.cpp:647
ProgForwardZernikeImages::checkPoint
virtual void checkPoint()
For very long programs, it may be needed to write checkpoints.
Definition: forward_zernike_images.cpp:934

◆ createWorkFiles()

void ProgForwardZernikeImages::createWorkFiles ( )
inlineprotected

Definition at line 230 of file forward_zernike_images.h.

230 { return Rerunable::createWorkFiles(resume, getInputMd()); }
Rerunable::createWorkFiles
virtual void createWorkFiles(bool resume, MetaData *md)
Definition: rerunable_program.h:49

◆ defineParams()

void ProgForwardZernikeImages::defineParams ( )
virtual

Define parameters.

Reimplemented from XmippMetadataProgram.

Definition at line 105 of file forward_zernike_images.cpp.

106 {
107  addUsageLine("Make a continuous angular assignment with deformations");
108  defaultComments["-i"].clear();
109  defaultComments["-i"].addComment("Metadata with initial alignment");
110  defaultComments["-o"].clear();
111  defaultComments["-o"].addComment("Metadata with the angular alignment and deformation parameters");
112  XmippMetadataProgram::defineParams();
113  addParamsLine(" --ref <volume> : Reference volume");
114  addParamsLine(" [--mask <m=\"\">] : Reference volume");
115  addParamsLine(" [--odir <outputDir=\".\">] : Output directory");
116  addParamsLine(" [--max_shift <s=-1>] : Maximum shift allowed in pixels");
117  addParamsLine(" [--max_angular_change <a=5>] : Maximum angular change allowed (in degrees)");
118  addParamsLine(" [--max_resolution <f=4>] : Maximum resolution (A)");
119  addParamsLine(" [--sampling <Ts=1>] : Sampling rate (A/pixel)");
120  addParamsLine(" [--Rmax <R=-1>] : Maximum radius (px). -1=Half of volume size");
121  addParamsLine(" [--RDef <r=-1>] : Maximum radius of the deformation (px). -1=Half of volume size");
122  addParamsLine(" [--l1 <l1=3>] : Degree Zernike Polynomials=1,2,3,...");
123  addParamsLine(" [--l2 <l2=2>] : Harmonical depth of the deformation=1,2,3,...");
124  addParamsLine(" [--step <step=1>] : Voxel index step");
125  addParamsLine(" [--useCTF] : Correct CTF");
126  addParamsLine(" [--optimizeAlignment] : Optimize alignment");
127  addParamsLine(" [--optimizeDeformation] : Optimize deformation");
128  addParamsLine(" [--optimizeDefocus] : Optimize defocus");
129  addParamsLine(" [--phaseFlipped] : Input images have been phase flipped");
130  addParamsLine(" [--regularization <l=0.01>] : Regularization weight");
131  addParamsLine(" [--blobr <b=4>] : Blob radius for forward mapping splatting");
132  addParamsLine(" [--image_mode <im=-1>] : Image mode (single, pairs, triplets). By default, it will be automatically identified.");
133  addParamsLine(" [--resume] : Resume processing");
134  addExampleLine("A typical use is:",false);
135  addExampleLine("xmipp_forward_zernike_images_priors -i anglesFromContinuousAssignment.xmd --ref reference.vol -o assigned_anglesAndDeformations.xmd --optimizeAlignment --optimizeDeformation --depth 1");
136 }
XmippProgram::addExampleLine
void addExampleLine(const char *example, bool verbatim=true)
Definition: xmipp_program.cpp:286
XmippProgram::addUsageLine
void addUsageLine(const char *line, bool verbatim=false)
Definition: xmipp_program.cpp:282
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305
XmippProgram::defaultComments
std::map< String, CommentList > defaultComments
Definition: xmipp_program.h:83

◆ deformVol()

void ProgForwardZernikeImages::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 1047 of file forward_zernike_images.cpp.

1049 {
1050  size_t idxY0=(VEC_XSIZE(clnm)-algn_params-ctf_params)/3;
1051  double Ncount=0.0;
1052  double modg=0.0;
1053  double diff2=0.0;
1054 
1055  def=0.0;
1056  size_t idxZ0=2*idxY0;
1057  double RmaxF=RmaxDef;
1058  double RmaxF2=RmaxF*RmaxF;
1059  double iRmaxF=1.0/RmaxF;
1060  Matrix2D<double> R_inv = R.inv();
1061 
1062  auto sz = idx_z_clnm.size();
1063  Matrix1D<int> l1, l2, n, m, idx_v;
1064 
1065  if (!idx_z_clnm.empty())
1066  {
1067  l1.initZeros(sz);
1068  l2.initZeros(sz);
1069  n.initZeros(sz);
1070  m.initZeros(sz);
1071  idx_v.initZeros(sz);
1072  for (auto j=0; j<sz; j++)
1073  {
1074  auto idx = idx_z_clnm[j];
1075  if (idx >= idxY0)
1076  idx -= idxY0;
1077 
1078  VEC_ELEM(idx_v,j) = idx;
1079  VEC_ELEM(l1,j) = VEC_ELEM(vL1, idx);
1080  VEC_ELEM(n,j) = VEC_ELEM(vN, idx);
1081  VEC_ELEM(l2,j) = VEC_ELEM(vL2, idx);
1082  VEC_ELEM(m,j) = VEC_ELEM(vM, idx);
1083  }
1084  }
1085 
1086  const auto &mVpos = vpos;
1087  const auto lastY = FINISHINGY(mVpos);
1088  for (int i=STARTINGY(mVpos); i<=lastY; i++)
1089  {
1090  double &gx = A2D_ELEM(df, i, 0);
1091  double &gy = A2D_ELEM(df, i, 1);
1092  double &gz = A2D_ELEM(df, i, 2);
1093  double r_x = A2D_ELEM(mVpos, i, 0);
1094  double r_y = A2D_ELEM(mVpos, i, 1);
1095  double r_z = A2D_ELEM(mVpos, i, 2);
1096  double xr = A2D_ELEM(mVpos, i, 3);
1097  double yr = A2D_ELEM(mVpos, i, 4);
1098  double zr = A2D_ELEM(mVpos, i, 5);
1099  double rr = A2D_ELEM(mVpos, i, 6);
1100 
1101  if (!idx_z_clnm.empty())
1102  {
1103  for (auto j = 0; j < sz; j++)
1104  {
1105  auto idx = VEC_ELEM(idx_v, j);
1106  auto aux_l2 = VEC_ELEM(l2, j);
1107  auto zsph = ZernikeSphericalHarmonics(VEC_ELEM(l1, j), VEC_ELEM(n, j),
1108  aux_l2, VEC_ELEM(m, j), xr, yr, zr, rr);
1109 
1110  auto diff_c_x = VEC_ELEM(clnm, idx) - VEC_ELEM(prev_clnm, idx);
1111  auto diff_c_y = VEC_ELEM(clnm, idx + idxY0) - VEC_ELEM(prev_clnm, idx + idxY0);
1112  auto i_diff_c_x = R_inv.mdata[0] * diff_c_x + R_inv.mdata[1] * diff_c_y;
1113  auto i_diff_c_y = R_inv.mdata[3] * diff_c_x + R_inv.mdata[4] * diff_c_y;
1114  auto i_diff_c_z = R_inv.mdata[6] * diff_c_x + R_inv.mdata[7] * diff_c_y;
1115  if (rr > 0 || aux_l2 == 0)
1116  {
1117  gx += i_diff_c_x * (zsph);
1118  gy += i_diff_c_y * (zsph);
1119  gz += i_diff_c_z * (zsph);
1120  }
1121  }
1122  }
1123 
1124  auto r_gx = R.mdata[0] * gx + R.mdata[1] * gy + R.mdata[2] * gz;
1125  auto r_gy = R.mdata[3] * gx + R.mdata[4] * gy + R.mdata[5] * gz;
1126 
1127  auto pos = std::array<double, 3>{};
1128  pos[0] = r_x + r_gx;
1129  pos[1] = r_y + r_gy;
1130  pos[2] = r_z;
1131 
1132  double voxel_mV = A2D_ELEM(mVpos, i, 7);
1133  splattingAtPos(pos, voxel_mV, mP, mV);
1134 
1135  modg += gx * gx + gy * gy + gz * gz;
1136  Ncount++;
1137  }
1138 
1139  def = sqrt(modg/Ncount);
1140  totalDeformation = def;
1141 }
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
Matrix2D::mdata
T * mdata
Definition: matrix2d.h:395
i
#define i
Definition: numerical_recipes.cpp:2493
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
ProgForwardZernikeImages::vpos
MultidimArray< double > vpos
Definition: forward_zernike_images.h:149
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
j
#define j
Definition: numerical_recipes.cpp:2493
m
int m
Definition: image_find_center.cpp:57
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
ZernikeSphericalHarmonics
double ZernikeSphericalHarmonics(int n, int m, double xr, double yr, double zr, double r)
Definition: numerical_tools.cpp:1054
ProgForwardZernikeImages::splattingAtPos
void splattingAtPos(std::array< double, 3 > r, double weight, MultidimArray< double > &mP, const MultidimArray< double > &mV)
Definition: forward_zernike_images.cpp:1174
ProgForwardZernikeImages::idx_z_clnm
std::vector< size_t > idx_z_clnm
Definition: forward_zernike_images.h:150
n
int * n
Definition: numerical_recipes.cpp:2229
ProgForwardZernikeImages::df
MultidimArray< double > df
Definition: forward_zernike_images.h:149

◆ fillVectorTerms()

void ProgForwardZernikeImages::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 991 of file forward_zernike_images.cpp.

993 {
994  int idx = 0;
995  vL1.initZeros(vecSize);
996  vN.initZeros(vecSize);
997  vL2.initZeros(vecSize);
998  vM.initZeros(vecSize);
999  for (int h=0; h<=l2; h++) {
1000  int totalSPH = 2*h+1;
1001  int aux = std::floor(totalSPH/2);
1002  for (int l=h; l<=l1; l+=2) {
1003  for (int m=0; m<totalSPH; m++) {
1004  VEC_ELEM(vL1,idx) = l;
1005  VEC_ELEM(vN,idx) = h;
1006  VEC_ELEM(vL2,idx) = h;
1007  VEC_ELEM(vM,idx) = m-aux;
1008  idx++;
1009  }
1010  }
1011  }
1012 }
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
m
int m
Definition: image_find_center.cpp:57

◆ finishProcessing()

void ProgForwardZernikeImages::finishProcessing ( )
virtual

Write the final parameters.

Reimplemented from XmippMetadataProgram.

Reimplemented in MpiProgForwardZernikeImages.

Definition at line 285 of file forward_zernike_images.cpp.

285  {
286  XmippMetadataProgram::finishProcessing();
287  rename(Rerunable::getFileName().c_str(), (fnOutDir + fn_out).c_str());
288 }
ProgForwardZernikeImages::fnOutDir
FileName fnOutDir
Output directory.
Definition: forward_zernike_images.h:51
Rerunable::getFileName
const FileName & getFileName() const
Definition: rerunable_program.h:80

◆ minimizepos()

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

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

Definition at line 958 of file forward_zernike_images.cpp.

959 {
960  int size = 0;
961  int prevSize = 0;
962  numCoefficients(L1,l2,size);
963  numCoefficients(L1,l2-1,prevSize);
964  int totalSize = (steps.size()-algn_params-ctf_params)/3;
965  if (l2 > 1)
966  {
967  for (int idx = prevSize; idx < size; idx++)
968  {
969  VEC_ELEM(steps, idx) = 1.;
970  VEC_ELEM(steps, idx + totalSize) = 1.;
971  if (num_images > 1)
972  {
973  VEC_ELEM(steps, idx + 2 * totalSize) = 1.;
974  }
975  }
976  }
977  else
978  {
979  for (int idx = 0; idx < size; idx++)
980  {
981  VEC_ELEM(steps, idx) = 1.;
982  VEC_ELEM(steps, idx + totalSize) = 1.;
983  if (num_images > 1)
984  {
985  VEC_ELEM(steps, idx + 2 * totalSize) = 1.;
986  }
987  }
988  }
989 }
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
ProgForwardZernikeImages::numCoefficients
void numCoefficients(int l1, int l2, int &vecSize)
Length of coefficients vector.
Definition: forward_zernike_images.cpp:942

◆ numCoefficients()

void ProgForwardZernikeImages::numCoefficients ( int  l1,
int  l2,
int &  vecSize 
)

Length of coefficients vector.

Definition at line 942 of file forward_zernike_images.cpp.

943 {
944  for (int h=0; h<=l2; h++)
945  {
946  int numSPH = 2*h+1;
947  int count=l1-h+1;
948  int numEven=(count>>1)+(count&1 && !(h&1));
949  if (h%2 == 0) {
950  vecSize += numSPH*numEven;
951  }
952  else {
953  vecSize += numSPH*(l1-h+1-numEven);
954  }
955  }
956 }

◆ preComputeDF()

void ProgForwardZernikeImages::preComputeDF ( )

Definition at line 1300 of file forward_zernike_images.cpp.

1301 {
1302  size_t idxY0=(VEC_XSIZE(clnm)-algn_params-ctf_params)/3;
1303  size_t idxZ0=2*idxY0;
1304  Matrix2D<double> R_inv = R.inv();
1305  const auto &mVpos = vpos;
1306  const auto lastY = FINISHINGY(mVpos);
1307  for (int i=STARTINGY(mVpos); i<=lastY; i++)
1308  {
1309  double &gx = A2D_ELEM(df, i, 0);
1310  double &gy = A2D_ELEM(df, i, 1);
1311  double &gz = A2D_ELEM(df, i, 2);
1312  double r_x = A2D_ELEM(mVpos, i, 0);
1313  double r_y = A2D_ELEM(mVpos, i, 1);
1314  double r_z = A2D_ELEM(mVpos, i, 2);
1315  double xr = A2D_ELEM(mVpos, i, 3);
1316  double yr = A2D_ELEM(mVpos, i, 4);
1317  double zr = A2D_ELEM(mVpos, i, 5);
1318  double rr = A2D_ELEM(mVpos, i, 6);
1319 
1320  if (!idx_z_clnm.empty())
1321  {
1322  for (int idx = 0; idx < idxY0; idx++)
1323  {
1324  auto aux_l2 = VEC_ELEM(vL2, idx);
1325  auto zsph = ZernikeSphericalHarmonics(VEC_ELEM(vL1, idx), VEC_ELEM(vN, idx),
1326  aux_l2, VEC_ELEM(vM, idx), xr, yr, zr, rr);
1327  auto c_x = VEC_ELEM(clnm, idx);
1328  auto c_y = VEC_ELEM(clnm, idx + idxY0);
1329  auto c_z = VEC_ELEM(clnm, idx + idxZ0);
1330  auto i_c_x = R_inv.mdata[0] * c_x + R_inv.mdata[1] * c_y + R_inv.mdata[2] * c_z;
1331  auto i_c_y = R_inv.mdata[3] * c_x + R_inv.mdata[4] * c_y + R_inv.mdata[5] * c_z;
1332  auto i_c_z = R_inv.mdata[6] * c_x + R_inv.mdata[7] * c_y + R_inv.mdata[8] * c_z;
1333  if (rr > 0 || aux_l2 == 0)
1334  {
1335  gx += i_c_x * (zsph);
1336  gy += i_c_y * (zsph);
1337  gz += i_c_z * (zsph);
1338  }
1339  }
1340  }
1341  }
1342 }
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
Matrix2D::mdata
T * mdata
Definition: matrix2d.h:395
i
#define i
Definition: numerical_recipes.cpp:2493
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
ProgForwardZernikeImages::vpos
MultidimArray< double > vpos
Definition: forward_zernike_images.h:149
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
ZernikeSphericalHarmonics
double ZernikeSphericalHarmonics(int n, int m, double xr, double yr, double zr, double r)
Definition: numerical_tools.cpp:1054
ProgForwardZernikeImages::idx_z_clnm
std::vector< size_t > idx_z_clnm
Definition: forward_zernike_images.h:150
ProgForwardZernikeImages::df
MultidimArray< double > df
Definition: forward_zernike_images.h:149

◆ preProcess()

void ProgForwardZernikeImages::preProcess ( )
virtual

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

Reimplemented from XmippMetadataProgram.

Definition at line 138 of file forward_zernike_images.cpp.

139 {
140  V.read(fnVolR);
141  V().setXmippOrigin();
142  Xdim=XSIZE(V());
143  Vdeformed().initZeros(V());
144  Vdeformed().setXmippOrigin();
145 
146  // Check execution mode (single, pair, or triplet)
147  if (image_mode > 0)
148  {
149  num_images = image_mode;
150  }
151  else{
152  num_images = 1; // Single image
153  if (getInputMd()->containsLabel(MDL_IMAGE1) && !getInputMd()->containsLabel(MDL_IMAGE2))
154  {
155  num_images = 2; // Pair
156  }
157  else if (getInputMd()->containsLabel(MDL_IMAGE1) && getInputMd()->containsLabel(MDL_IMAGE2))
158  {
159  num_images = 3; // Triplet
160  }
161  }
162 
163  // Preallocate vectors (Size depends on image number)
164  fnImage.resize(num_images, "");
165  I.resize(num_images); Ifiltered.resize(num_images); Ifilteredp.resize(num_images);
166  P.resize(num_images);
167  old_rot.resize(num_images, 0.); old_tilt.resize(num_images, 0.); old_psi.resize(num_images, 0.);
168  deltaRot.resize(num_images, 0.); deltaTilt.resize(num_images, 0.); deltaPsi.resize(num_images, 0.);
169  old_shiftX.resize(num_images, 0.); old_shiftY.resize(num_images, 0.);
170  deltaX.resize(num_images, 0.); deltaY.resize(num_images, 0.);
171  old_defocusU.resize(num_images, 0.); old_defocusV.resize(num_images, 0.); old_defocusAngle.resize(num_images, 0.);
172  deltaDefocusU.resize(num_images, 0.); deltaDefocusV.resize(num_images, 0.); deltaDefocusAngle.resize(num_images, 0.);
173  currentDefocusU.resize(num_images, 0.); currentDefocusV.resize(num_images, 0.); currentAngle.resize(num_images, 0.);
174 
175  switch (num_images)
176  {
177  case 2:
178  Ifilteredp[0]().initZeros(Xdim,Xdim); Ifilteredp[1]().initZeros(Xdim,Xdim);
179  Ifilteredp[0]().setXmippOrigin(); Ifilteredp[1]().setXmippOrigin();
180  P[0]().initZeros(Xdim,Xdim); P[1]().initZeros(Xdim,Xdim);
181  break;
182  case 3:
183  Ifilteredp[0]().initZeros(Xdim,Xdim); Ifilteredp[1]().initZeros(Xdim,Xdim); Ifilteredp[2]().initZeros(Xdim,Xdim);
184  Ifilteredp[0]().setXmippOrigin(); Ifilteredp[1]().setXmippOrigin(); Ifilteredp[2]().setXmippOrigin();
185  P[0]().initZeros(Xdim,Xdim); P[1]().initZeros(Xdim,Xdim); P[2]().initZeros(Xdim,Xdim);
186  break;
187 
188  default:
189  Ifilteredp[0]().initZeros(Xdim,Xdim);
190  Ifilteredp[0]().setXmippOrigin();
191  P[0]().initZeros(Xdim,Xdim);
192  break;
193  }
194 
195  if (RmaxDef<0)
196  RmaxDef = Xdim/2;
197 
198  // Read Reference mask if avalaible (otherwise sphere of radius RmaxDef is used)
199  Mask mask;
200  mask.type = BINARY_CIRCULAR_MASK;
201  mask.mode = INNER_MASK;
202  if (fnMaskR != "") {
203  Image<double> aux;
204  aux.read(fnMaskR);
205  typeCast(aux(), V_mask);
206  V_mask.setXmippOrigin();
207  double Rmax2 = RmaxDef*RmaxDef;
208  for (int k=STARTINGZ(V_mask); k<=FINISHINGZ(V_mask); k++) {
209  for (int i=STARTINGY(V_mask); i<=FINISHINGY(V_mask); i++) {
210  for (int j=STARTINGX(V_mask); j<=FINISHINGX(V_mask); j++) {
211  double r2 = k*k + i*i + j*j;
212  if (r2>=Rmax2)
213  A3D_ELEM(V_mask,k,i,j) = 0;
214  }
215  }
216  }
217  }
218  else {
219  mask.R1 = RmaxDef;
220  mask.generate_mask(V());
221  V_mask = mask.get_binary_mask();
222  V_mask.setXmippOrigin();
223  }
224 
225  // Total Volume Mass (Inside Mask)
226  sumV = 0;
227  for (int k = STARTINGZ(V()); k <= FINISHINGZ(V()); k += loop_step)
228  {
229  for (int i = STARTINGY(V()); i <= FINISHINGY(V()); i += loop_step)
230  {
231  for (int j = STARTINGX(V()); j <= FINISHINGX(V()); j += loop_step)
232  {
233  if (A3D_ELEM(V_mask, k, i, j) == 1)
234  {
235  sumV++;
236  }
237  }
238  }
239  }
240 
241  // Construct mask
242  if (Rmax<0)
243  Rmax=Xdim/2;
244  mask.R1 = Rmax;
245  mask.generate_mask(Xdim,Xdim);
246  mask2D=mask.get_binary_mask();
247 
248  // Low pass filter
249  filter.FilterBand=LOWPASS;
250  filter.w1=Ts/maxResol;
251  filter.raised_w=0.02;
252  filterp.FilterBand=LOWPASS;
253  filterp.FilterShape=REALGAUSSIAN;
254  filterp.w1=1;
255 
256  // Transformation matrix
257  A1.initIdentity(3);
258  A2.initIdentity(3);
259  A3.initIdentity(3);
260 
261  // CTF Filter
262  FilterCTF1.FilterBand = CTF;
263  FilterCTF1.ctf.enable_CTFnoise = false;
264  FilterCTF1.ctf.produceSideInfo();
265  FilterCTF2.FilterBand = CTF;
266  FilterCTF2.ctf.enable_CTFnoise = false;
267  FilterCTF2.ctf.produceSideInfo();
268  FilterCTF3.FilterBand = CTF;
269  FilterCTF3.ctf.enable_CTFnoise = false;
270  FilterCTF3.ctf.produceSideInfo();
271 
272  vecSize = 0;
273  numCoefficients(L1,L2,vecSize);
274  fillVectorTerms(L1,L2,vL1,vN,vL2,vM);
275 
276  createWorkFiles();
277 
278  // Blob
279  blob.radius = blob_r; // Blob radius in voxels
280  blob.order = 2; // Order of the Bessel function
281  blob.alpha = 7.05; // Smoothness parameter
282 
283 }
ProgForwardZernikeImages::fillVectorTerms
void fillVectorTerms(int l1, int l2, Matrix1D< int > &vL1, Matrix1D< int > &vN, Matrix1D< int > &vL2, Matrix1D< int > &vM)
Zernike and SPH coefficients allocation.
Definition: forward_zernike_images.cpp:991
ProgForwardZernikeImages::Ifiltered
std::vector< Image< double > > Ifiltered
Definition: forward_zernike_images.h:105
ProgForwardZernikeImages::fnImage
std::vector< FileName > fnImage
Definition: forward_zernike_images.h:101
ProgForwardZernikeImages::old_shiftY
std::vector< double > old_shiftY
Definition: forward_zernike_images.h:116
blobtype::alpha
double alpha
Smoothness parameter.
Definition: blobs.h:121
ProgForwardZernikeImages::old_psi
std::vector< double > old_psi
Definition: forward_zernike_images.h:114
FourierFilter::ctf
CTFDescription ctf
Definition: fourier_filter.h:134
ProgForwardZernikeImages::old_tilt
std::vector< double > old_tilt
Definition: forward_zernike_images.h:114
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
ProgForwardZernikeImages::old_defocusV
std::vector< double > old_defocusV
Definition: forward_zernike_images.h:122
Mask
Definition: mask.h:360
ProgForwardZernikeImages::deltaPsi
std::vector< double > deltaPsi
Definition: forward_zernike_images.h:114
ProgForwardZernikeImages::I
std::vector< Image< double > > I
Definition: forward_zernike_images.h:104
ProgForwardZernikeImages::deltaTilt
std::vector< double > deltaTilt
Definition: forward_zernike_images.h:114
ProgForwardZernikeImages::old_defocusU
std::vector< double > old_defocusU
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::Ifilteredp
std::vector< Image< double > > Ifilteredp
Definition: forward_zernike_images.h:106
ProgForwardZernikeImages::deltaDefocusAngle
std::vector< double > deltaDefocusAngle
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::deltaDefocusU
std::vector< double > deltaDefocusU
Definition: forward_zernike_images.h:122
FINISHINGZ
#define FINISHINGZ(v)
Definition: multidim_array_base.h:71
ProgForwardZernikeImages::deltaX
std::vector< double > deltaX
Definition: forward_zernike_images.h:116
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
i
#define i
Definition: numerical_recipes.cpp:2493
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
ProgForwardZernikeImages::deltaY
std::vector< double > deltaY
Definition: forward_zernike_images.h:116
ProgForwardZernikeImages::P
std::vector< Image< double > > P
Definition: forward_zernike_images.h:108
ProgForwardZernikeImages::currentAngle
std::vector< double > currentAngle
Definition: forward_zernike_images.h:124
CTF
#define CTF
Definition: fourier_filter.h:81
ProgForwardZernikeImages::old_rot
std::vector< double > old_rot
Definition: forward_zernike_images.h:114
Mask::R1
double R1
Definition: mask.h:413
ProgForwardZernikeImages::old_defocusAngle
std::vector< double > old_defocusAngle
Definition: forward_zernike_images.h:122
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
MDL_IMAGE2
Image associated to this object (std::string)
Definition: metadata_label.h:256
Mask::type
int type
Definition: mask.h:402
MDL_IMAGE1
Image associated to this object (std::string)
Definition: metadata_label.h:255
j
#define j
Definition: numerical_recipes.cpp:2493
ProgForwardZernikeImages::numCoefficients
void numCoefficients(int l1, int l2, int &vecSize)
Length of coefficients vector.
Definition: forward_zernike_images.cpp:942
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
ProgForwardZernikeImages::old_shiftX
std::vector< double > old_shiftX
Definition: forward_zernike_images.h:116
Mask::generate_mask
void generate_mask(bool apply_geo=false)
Definition: mask.cpp:1577
typeCast
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
BINARY_CIRCULAR_MASK
#define BINARY_CIRCULAR_MASK
Definition: mask.h:365
CTFDescription::produceSideInfo
void produceSideInfo()
Produce Side information.
Definition: ctf.cpp:1392
blobtype::order
int order
Derivation order and Bessel function order.
Definition: blobs.h:118
ImageBase::read
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
Definition: xmipp_image_base.cpp:119
r2
float r2
Definition: image_find_center.cpp:54
REALGAUSSIAN
#define REALGAUSSIAN
Definition: fourier_filter.h:87
CTFDescription1D::enable_CTFnoise
bool enable_CTFnoise
Enable CTFnoise part.
Definition: ctf.h:273
Mask::get_binary_mask
const MultidimArray< int > & get_binary_mask() const
Definition: mask.h:707
ProgForwardZernikeImages::currentDefocusU
std::vector< double > currentDefocusU
Definition: forward_zernike_images.h:124
ProgForwardZernikeImages::deltaDefocusV
std::vector< double > deltaDefocusV
Definition: forward_zernike_images.h:122
STARTINGZ
#define STARTINGZ(v)
Definition: multidim_array_base.h:67
blobtype::radius
double radius
Spatial radius in Universal System units.
Definition: blobs.h:115
LOWPASS
#define LOWPASS
Definition: fourier_filter.h:77
ProgForwardZernikeImages::currentDefocusV
std::vector< double > currentDefocusV
Definition: forward_zernike_images.h:124
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673
Mask::mode
int mode
Definition: mask.h:407
ProgForwardZernikeImages::deltaRot
std::vector< double > deltaRot
Definition: forward_zernike_images.h:114
INNER_MASK
constexpr int INNER_MASK
Definition: mask.h:47

◆ processImage()

void ProgForwardZernikeImages::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 641 of file forward_zernike_images.cpp.

642 {
643  Matrix1D<double> steps;
644  // totalSize = 3*num_Z_coeff + num_images * 2 shifts + num_images * 3 angles + num_images * 3 CTF
645  algn_params = num_images * 5;
646  ctf_params = num_images * 3;
647  int totalSize = 3*vecSize + algn_params + ctf_params;
648  p.initZeros(totalSize);
649  clnm.initZeros(totalSize);
650  prev_clnm.initZeros(totalSize);
651 
652  // Init positions and deformation field
653  vpos.initZeros(sumV, 8);
654  df.initZeros(sumV, 3);
655 
656  flagEnabled=1;
657 
658  rowIn.getValueOrDefault(MDL_IMAGE, fnImage[0], "");
659  rowIn.getValueOrDefault(MDL_ANGLE_ROT, old_rot[0], 0.0);
660  rowIn.getValueOrDefault(MDL_ANGLE_TILT, old_tilt[0], 0.0);
661  rowIn.getValueOrDefault(MDL_ANGLE_PSI, old_psi[0], 0.0);
662  rowIn.getValueOrDefault(MDL_SHIFT_X, old_shiftX[0], 0.0);
663  rowIn.getValueOrDefault(MDL_SHIFT_Y, old_shiftY[0], 0.0);
664  I[0].read(fnImage[0]);
665  I[0]().setXmippOrigin();
666  Ifiltered[0]() = I[0]();
667  filter.applyMaskSpace(Ifiltered[0]());
668  rotatePositions(old_rot[0], old_tilt[0], old_psi[0]);
669 
670  switch (num_images)
671  {
672  case 2:
673  rowIn.getValueOrDefault(MDL_IMAGE1, fnImage[1], "");
674  rowIn.getValueOrDefault(MDL_ANGLE_ROT2, old_rot[1], 0.0);
675  rowIn.getValueOrDefault(MDL_ANGLE_TILT2, old_tilt[1], 0.0);
676  rowIn.getValueOrDefault(MDL_ANGLE_PSI2, old_psi[1], 0.0);
677  rowIn.getValueOrDefault(MDL_SHIFT_X2, old_shiftX[1], 0.0);
678  rowIn.getValueOrDefault(MDL_SHIFT_Y2, old_shiftY[1], 0.0);
679  I[1].read(fnImage[1]);
680  I[1]().setXmippOrigin();
681  Ifiltered[1]() = I[1]();
682  filter.applyMaskSpace(Ifiltered[1]());
683 
684  if (verbose >= 2)
685  std::cout << "Processing Pair (" << fnImage[0] << "," << fnImage[1] << ")" << std::endl;
686  break;
687 
688  case 3:
689  rowIn.getValueOrDefault(MDL_IMAGE1, fnImage[1], "");
690  rowIn.getValueOrDefault(MDL_ANGLE_ROT2, old_rot[1], 0.0);
691  rowIn.getValueOrDefault(MDL_ANGLE_TILT2, old_tilt[1], 0.0);
692  rowIn.getValueOrDefault(MDL_ANGLE_PSI2, old_psi[1], 0.0);
693  rowIn.getValueOrDefault(MDL_SHIFT_X2, old_shiftX[1], 0.0);
694  rowIn.getValueOrDefault(MDL_SHIFT_Y2, old_shiftY[1], 0.0);
695  I[1].read(fnImage[1]);
696  I[1]().setXmippOrigin();
697  Ifiltered[1]() = I[1]();
698  filter.applyMaskSpace(Ifiltered[1]());
699 
700  rowIn.getValueOrDefault(MDL_IMAGE2, fnImage[2], "");
701  rowIn.getValueOrDefault(MDL_ANGLE_ROT3, old_rot[2], 0.0);
702  rowIn.getValueOrDefault(MDL_ANGLE_TILT3, old_tilt[2], 0.0);
703  rowIn.getValueOrDefault(MDL_ANGLE_PSI3, old_psi[2], 0.0);
704  rowIn.getValueOrDefault(MDL_SHIFT_X3, old_shiftX[2], 0.0);
705  rowIn.getValueOrDefault(MDL_SHIFT_Y3, old_shiftY[2], 0.0);
706  I[2].read(fnImage[2]);
707  I[2]().setXmippOrigin();
708  Ifiltered[2]() = I[2]();
709  filter.applyMaskSpace(Ifiltered[2]());
710 
711  if (verbose >= 2)
712  std::cout << "Processing Triplet (" << fnImage[0] << "," << fnImage[1] << "," << fnImage[2] << ")" << std::endl;
713  break;
714 
715  default:
716  if (verbose >= 2)
717  std::cout << "Processing Image (" << fnImage[0] << ")" << std::endl;
718  break;
719  }
720 
721  if (rowIn.containsLabel(MDL_FLIP))
722  rowIn.getValue(MDL_FLIP,old_flip);
723  else
724  old_flip = false;
725 
726  prior_deformation = 0.0;
727  if (rowIn.containsLabel(MDL_SPH_COEFFICIENTS))
728  {
729  std::vector<double> vectortemp;
730  rowIn.getValue(MDL_SPH_COEFFICIENTS, vectortemp);
731  rowIn.getValueOrDefault(MDL_SPH_DEFORMATION, prior_deformation, 0.0);
732  for (int i=0; i<3*vecSize; i++)
733  {
734  clnm[i] = vectortemp[i];
735  }
736  if (optimizeDeformation)
737  rotateCoefficients<Direction::ROTATE>();
738  p = clnm;
739  prev_clnm = clnm;
740  preComputeDF();
741  }
742 
743  // FIXME: Add defocus per image and make CTF correction available
744  if ((rowIn.containsLabel(MDL_CTF_DEFOCUSU) || rowIn.containsLabel(MDL_CTF_MODEL)) && useCTF)
745  {
746  hasCTF=true;
747  FilterCTF1.ctf.readFromMdRow(rowIn);
748  FilterCTF1.ctf.Tm = Ts;
749  FilterCTF1.ctf.produceSideInfo();
750  old_defocusU[0]=FilterCTF1.ctf.DeltafU;
751  old_defocusV[0]=FilterCTF1.ctf.DeltafV;
752  old_defocusAngle[0]=FilterCTF1.ctf.azimuthal_angle;
753  }
754  else
755  hasCTF=false;
756 
757  // If deformation is not optimized, do a single iteration
758  //? Si usamos priors es mejor ir poco a poco, ir poco a poco pero usar todos los coeffs cada vez (mas lento)
759  //? O dar solo una iteracion con todos los coeffs?
760  int h = 1;
761  // if (!optimizeDeformation)
762  // if (rowIn.containsLabel(MDL_SPH_COEFFICIENTS))
763  // h = L2;
764 
765  for (;h<=L2;h++)
766  {
767  if (verbose>=2)
768  {
769  std::cout<<std::endl;
770  std::cout<<"------------------------------ Basis Degrees: ("<<L1<<","<<h<<") ----------------------------"<<std::endl;
771  }
772  steps.clear();
773  steps.initZeros(totalSize);
774 
775  // Optimize
776  double cost=-1;
777  try
778  {
779  cost=1e38;
780  int iter;
781  if (optimizeAlignment)
782  {
783  int init = steps.size() - algn_params - ctf_params;
784  int end = steps.size() - ctf_params;
785  for (int i = init; i < end; i++)
786  steps(i)=1.;
787  }
788  if (optimizeDefocus)
789  {
790  int init = steps.size() - ctf_params;
791  int end = steps.size();
792  for (int i = init; i < end; i++)
793  steps(i)=1.;
794  }
795  if (optimizeDeformation)
796  {
797  minimizepos(L1,h,steps);
798  }
799  steps_cp = steps;
800  powellOptimizer(p, 1, totalSize, &continuousZernikeCost, this, 0.1, cost, iter, steps, verbose>=2);
801 
802  if (verbose>=3)
803  {
804  showOptimization = true;
805  continuousZernikeCost(p.adaptForNumericalRecipes(),this);
806  showOptimization = false;
807  }
808 
809  if (cost>0)
810  {
811  flagEnabled=-1;
812  p.initZeros();
813  }
814  cost=-cost;
815  correlation=cost;
816  if (verbose>=2)
817  {
818  std::cout<<std::endl;
819  for (int j=1;j<4;j++)
820  {
821  switch (j)
822  {
823  case 1:
824  std::cout << "X Coefficients=(";
825  break;
826  case 2:
827  std::cout << "Y Coefficients=(";
828  break;
829  case 3:
830  std::cout << "Z Coefficients=(";
831  break;
832  default:
833  // The code will never reach the default case
834  break;
835  }
836  for (int i=(j-1)*vecSize;i<j*vecSize;i++)
837  {
838  std::cout << p(i);
839  if (i<j*vecSize-1)
840  std::cout << ",";
841  }
842  std::cout << ")" << std::endl;
843  }
844  std::cout << "Radius=" << RmaxDef << std::endl;
845  std::cout << " Dshift=(" << p(totalSize-5) << "," << p(totalSize-4) << ") "
846  << "Drot=" << p(totalSize-3) << " Dtilt=" << p(totalSize-2)
847  << " Dpsi=" << p(totalSize-1) << std::endl;
848  std::cout << " Total deformation=" << totalDeformation << std::endl;
849  std::cout<<std::endl;
850  }
851  }
852  catch (XmippError &XE)
853  {
854  std::cerr << XE.what() << std::endl;
855  std::cerr << "Warning: Cannot refine " << fnImg << std::endl;
856  flagEnabled=-1;
857  }
858  }
859 
860  clnm = p;
861  if (num_images == 1 && optimizeDeformation)
862  {
863  rotateCoefficients<Direction::UNROTATE>();
864  }
865 
866  //AJ NEW
867  writeImageParameters(rowOut);
868  //END AJ
869 
870 }
ProgForwardZernikeImages::Ifiltered
std::vector< Image< double > > Ifiltered
Definition: forward_zernike_images.h:105
ProgForwardZernikeImages::minimizepos
void minimizepos(int L1, int l2, Matrix1D< double > &steps)
Determine the positions to be minimize of a vector containing spherical harmonic coefficients.
Definition: forward_zernike_images.cpp:958
ProgForwardZernikeImages::fnImage
std::vector< FileName > fnImage
Definition: forward_zernike_images.h:101
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
ProgForwardZernikeImages::old_shiftY
std::vector< double > old_shiftY
Definition: forward_zernike_images.h:116
MDL_CTF_DEFOCUSU
Defocus U (Angstroms)
Definition: metadata_label.h:135
MDL_ANGLE_TILT2
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:60
MDL_SHIFT_X3
Shift for the image in the X axis (double)
Definition: metadata_label.h:433
MDL_SHIFT_Y2
Shift for the image in the Y axis (double)
Definition: metadata_label.h:436
ProgForwardZernikeImages::old_psi
std::vector< double > old_psi
Definition: forward_zernike_images.h:114
FourierFilter::ctf
CTFDescription ctf
Definition: fourier_filter.h:134
ProgForwardZernikeImages::old_tilt
std::vector< double > old_tilt
Definition: forward_zernike_images.h:114
Matrix1D::clear
void clear()
Definition: matrix1d.cpp:67
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
ProgForwardZernikeImages::rotatePositions
void rotatePositions(double rot, double tilt, double psi)
Definition: forward_zernike_images.cpp:1257
ProgForwardZernikeImages::old_defocusV
std::vector< double > old_defocusV
Definition: forward_zernike_images.h:122
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
ProgForwardZernikeImages::I
std::vector< Image< double > > I
Definition: forward_zernike_images.h:104
CTFDescription::DeltafV
double DeltafV
Defocus in V (in Angstroms). Negative values are underfocused.
Definition: ctf.h:830
MDL_SHIFT_X
Shift for the image in the X axis (double)
Definition: metadata_label.h:431
ProgForwardZernikeImages::old_defocusU
std::vector< double > old_defocusU
Definition: forward_zernike_images.h:122
MDL_ANGLE_TILT3
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:61
MDL_ANGLE_PSI
Special label to be used when gathering MDs in MpiMetadataPrograms.
Definition: metadata_label.h:51
MDL_CTF_MODEL
Name for the CTF Model (std::string)
Definition: metadata_label.h:129
powellOptimizer
void powellOptimizer(Matrix1D< double > &p, int i0, int n, double(*f)(double *x, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show)
Definition: numerical_tools.cpp:44
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
ProgForwardZernikeImages::writeImageParameters
virtual void writeImageParameters(MDRow &row)
Definition: forward_zernike_images.cpp:873
i
#define i
Definition: numerical_recipes.cpp:2493
MDL_SHIFT_Y3
Shift for the image in the Y axis (double)
Definition: metadata_label.h:437
CTFDescription::DeltafU
double DeltafU
Global gain. By default, 1.
Definition: ctf.h:828
continuousZernikeCost
double continuousZernikeCost(double *x, void *_prm)
Definition: forward_zernike_images.cpp:520
CTFDescription::azimuthal_angle
double azimuthal_angle
Azimuthal angle (between X and U) in degrees.
Definition: ctf.h:832
MDL_ANGLE_ROT2
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:56
MDRow::getValue
T & getValue(MDLabel label)
Definition: metadata_row_base.h:97
MDL_SPH_DEFORMATION
Deformation in voxels.
Definition: metadata_label.h:449
ProgForwardZernikeImages::old_rot
std::vector< double > old_rot
Definition: forward_zernike_images.h:114
MDL_FLIP
Flip the image? (bool)
Definition: metadata_label.h:236
ProgForwardZernikeImages::old_defocusAngle
std::vector< double > old_defocusAngle
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::vpos
MultidimArray< double > vpos
Definition: forward_zernike_images.h:149
MDL_IMAGE2
Image associated to this object (std::string)
Definition: metadata_label.h:256
CTFDescription1D::Tm
double Tm
Sampling rate (A/pixel)
Definition: ctf.h:240
MDL_IMAGE1
Image associated to this object (std::string)
Definition: metadata_label.h:255
CTFDescription::readFromMdRow
void readFromMdRow(const MDRow &row, bool disable_if_not_K=true)
Definition: ctf.cpp:1172
XmippProgram::verbose
int verbose
Verbosity level.
Definition: xmipp_program.h:161
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
MDL_ANGLE_PSI2
Psi angle of an image (double,degrees)
Definition: metadata_label.h:52
j
#define j
Definition: numerical_recipes.cpp:2493
steps
double steps
Definition: numerical_recipes.cpp:5166
MDRow::getValueOrDefault
const T & getValueOrDefault(MDLabel label, const T &def) const
Definition: metadata_row_base.h:129
MDL_SHIFT_X2
Shift for the image in the X axis (double)
Definition: metadata_label.h:432
ProgForwardZernikeImages::old_shiftX
std::vector< double > old_shiftX
Definition: forward_zernike_images.h:116
MDRow::containsLabel
virtual bool containsLabel(MDLabel label) const =0
CTFDescription::produceSideInfo
void produceSideInfo()
Produce Side information.
Definition: ctf.cpp:1392
MDL_SPH_COEFFICIENTS
Deformation coefficients.
Definition: metadata_label.h:448
MDL_ANGLE_ROT3
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:57
MDL_SHIFT_Y
Shift for the image in the Y axis (double)
Definition: metadata_label.h:435
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
Matrix1D::adaptForNumericalRecipes
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
MDL_ANGLE_PSI3
Psi angle of an image (double,degrees)
Definition: metadata_label.h:53
MDL_IMAGE
Name of an image (std::string)
Definition: metadata_label.h:246
ProgForwardZernikeImages::df
MultidimArray< double > df
Definition: forward_zernike_images.h:149
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710

◆ readParams()

void ProgForwardZernikeImages::readParams ( )
virtual

Read argument from command line.

Reimplemented from XmippMetadataProgram.

Definition at line 46 of file forward_zernike_images.cpp.

47 {
48  XmippMetadataProgram::readParams();
49  fnVolR = getParam("--ref");
50  fnMaskR = getParam("--mask");
51  fnOutDir = getParam("--odir");
52  maxShift = getDoubleParam("--max_shift");
53  maxAngularChange = getDoubleParam("--max_angular_change");
54  maxResol = getDoubleParam("--max_resolution");
55  Ts = getDoubleParam("--sampling");
56  Rmax = getIntParam("--Rmax");
57  RmaxDef = getIntParam("--RDef");
58  optimizeAlignment = checkParam("--optimizeAlignment");
59  optimizeDeformation = checkParam("--optimizeDeformation");
60  optimizeDefocus = checkParam("--optimizeDefocus");
61  phaseFlipped = checkParam("--phaseFlipped");
62  useCTF = checkParam("--useCTF");
63  L1 = getIntParam("--l1");
64  L2 = getIntParam("--l2");
65  loop_step = getIntParam("--step");
66  lambda = getDoubleParam("--regularization");
67  image_mode = getIntParam("--image_mode");
68  resume = checkParam("--resume");
69  Rerunable::setFileName(fnOutDir + "/sphDone.xmd");
70  blob_r = getDoubleParam("--blobr");
71  keep_input_columns = true;
72 }
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
ProgForwardZernikeImages::fnOutDir
FileName fnOutDir
Output directory.
Definition: forward_zernike_images.h:51
Rerunable::setFileName
void setFileName(const FileName &fn)
Definition: rerunable_program.h:82
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
XmippMetadataProgram::keep_input_columns
bool keep_input_columns
Keep input metadata columns.
Definition: xmipp_metadata_program.h:105
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331

◆ removePixels()

void ProgForwardZernikeImages::removePixels ( )

Definition at line 1224 of file forward_zernike_images.cpp.

1225 {
1226  auto &mI = I[0]();
1227  MultidimArray<double> aux;
1228  aux.initZeros(mI);
1229  aux.setXmippOrigin();
1230  const MultidimArray<double> &mV=V();
1231 
1232  for (int i = STARTINGY(mI); i <= FINISHINGY(mI); i += loop_step)
1233  {
1234  for (int j = STARTINGX(mI); j <= FINISHINGX(mI); j += loop_step)
1235  {
1236 
1237  int i0 = XMIPP_MAX(FLOOR(i - loop_step), STARTINGY(mI));
1238  int iF = XMIPP_MIN(CEIL(i + loop_step), FINISHINGY(mI));
1239  int j0 = XMIPP_MAX(FLOOR(j - loop_step), STARTINGX(mI));
1240  int jF = XMIPP_MIN(CEIL(j + loop_step), FINISHINGX(mI));
1241  double weight = A2D_ELEM(mI, i, j);
1242 
1243  for (int img_i = i0; img_i <= iF; img_i++)
1244  {
1245  double y_val = bspline1(img_i - i);
1246  for (int img_j = j0; img_j <= jF; img_j++)
1247  {
1248  A2D_ELEM(aux, img_i, img_j) += weight * bspline1(img_j - j) * y_val;
1249  }
1250  }
1251  }
1252  }
1253 
1254  mI = aux;
1255 }
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
XMIPP_MAX
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
ProgForwardZernikeImages::I
std::vector< Image< double > > I
Definition: forward_zernike_images.h:104
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
i
#define i
Definition: numerical_recipes.cpp:2493
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
CEIL
#define CEIL(x)
Definition: xmipp_macros.h:225
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
j
#define j
Definition: numerical_recipes.cpp:2493
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723

◆ rotateCoefficients()

template<ProgForwardZernikeImages::Direction DIRECTION>
void ProgForwardZernikeImages::rotateCoefficients ( )

Definition at line 1024 of file forward_zernike_images.cpp.

1024  {
1025  int pos = 3*vecSize;
1026  size_t idxY0=(VEC_XSIZE(clnm)-algn_params-ctf_params)/3;
1027  size_t idxZ0=2*idxY0;
1028 
1029  double rot = old_rot[0]+p(pos+2);
1030  double tilt = old_tilt[0]+p(pos+3);
1031  double psi = old_psi[0]+p(pos+4);
1032 
1033  Matrix2D<double> R;
1034  R.initIdentity(3);
1035  Euler_angles2matrix(rot, tilt, psi, R, false);
1036  if (DIRECTION == Direction::UNROTATE)
1037  R = R.inv();
1038  Matrix1D<double> c;
1039  c.initZeros(3);
1040  for (size_t idx=0; idx<idxY0; idx++) {
1041  XX(c) = VEC_ELEM(clnm,idx); YY(c) = VEC_ELEM(clnm,idx+idxY0); ZZ(c) = VEC_ELEM(clnm,idx+idxZ0);
1042  c = R * c;
1043  VEC_ELEM(clnm,idx) = XX(c); VEC_ELEM(clnm,idx+idxY0) = YY(c); VEC_ELEM(clnm,idx+idxZ0) = ZZ(c);
1044  }
1045 }
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
ProgForwardZernikeImages::old_psi
std::vector< double > old_psi
Definition: forward_zernike_images.h:114
ProgForwardZernikeImages::old_tilt
std::vector< double > old_tilt
Definition: forward_zernike_images.h:114
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
Euler_angles2matrix
void Euler_angles2matrix(T alpha, T beta, T gamma, Matrix2D< T > &A, bool homogeneous)
Definition: geometry.cpp:624
c
doublereal * c
Definition: numerical_recipes.cpp:2230
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
XX
#define XX(v)
Definition: matrix1d.h:85
ProgForwardZernikeImages::old_rot
std::vector< double > old_rot
Definition: forward_zernike_images.h:114
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
YY
#define YY(v)
Definition: matrix1d.h:93
alglib::psi
double psi(const double x)
Definition: specialfunctions.cpp:2492
ZZ
#define ZZ(v)
Definition: matrix1d.h:101
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673

◆ rotatePositions()

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

Definition at line 1257 of file forward_zernike_images.cpp.

1258 {
1259  R.initIdentity(3);
1260  Euler_angles2matrix(rot, tilt, psi, R, false);
1261 
1262  const MultidimArray<double> &mV=V();
1263 
1264  const auto lastZ = FINISHINGZ(mV);
1265  const auto lastY = FINISHINGY(mV);
1266  const auto lastX = FINISHINGX(mV);
1267  size_t count = 0;
1268  double iRmaxF=1.0/RmaxDef;
1269  for (int k=STARTINGZ(mV); k<=lastZ; k+=loop_step)
1270  {
1271  for (int i=STARTINGY(mV); i<=lastY; i+=loop_step)
1272  {
1273  for (int j=STARTINGX(mV); j<=lastX; j+=loop_step)
1274  {
1275  if (A3D_ELEM(V_mask,k,i,j) == 1) {
1276  double x = j;
1277  double y = i;
1278  double z = k;
1279  double r_x = R.mdata[0] * x + R.mdata[1] * y + R.mdata[2] * z;
1280  double r_y = R.mdata[3] * x + R.mdata[4] * y + R.mdata[5] * z;
1281  double r_z = R.mdata[6] * x + R.mdata[7] * y + R.mdata[8] * z;
1282 
1283  A2D_ELEM(vpos, count, 0) = r_x;
1284  A2D_ELEM(vpos, count, 1) = r_y;
1285  A2D_ELEM(vpos, count, 2) = r_z;
1286  A2D_ELEM(vpos, count, 3) = j * iRmaxF;
1287  A2D_ELEM(vpos, count, 4) = i * iRmaxF;
1288  A2D_ELEM(vpos, count, 5) = z * iRmaxF;
1289  A2D_ELEM(vpos, count, 6) = sqrt(x*x + y*y + z*z) * iRmaxF;
1290  A2D_ELEM(vpos, count, 7) = A3D_ELEM(mV, k, i, j);
1291 
1292  count++;
1293  }
1294  }
1295  }
1296  }
1297 }
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
Euler_angles2matrix
void Euler_angles2matrix(T alpha, T beta, T gamma, Matrix2D< T > &A, bool homogeneous)
Definition: geometry.cpp:624
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
y
static double * y
Definition: numerical_recipes.cpp:8487
Matrix2D::mdata
T * mdata
Definition: matrix2d.h:395
FINISHINGZ
#define FINISHINGZ(v)
Definition: multidim_array_base.h:71
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
ProgForwardZernikeImages::vpos
MultidimArray< double > vpos
Definition: forward_zernike_images.h:149
z
double z
Definition: numerical_recipes.cpp:8490
j
#define j
Definition: numerical_recipes.cpp:2493
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
alglib::psi
double psi(const double x)
Definition: specialfunctions.cpp:2492
STARTINGZ
#define STARTINGZ(v)
Definition: multidim_array_base.h:67
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673

◆ show()

void ProgForwardZernikeImages::show ( )

Show.

Definition at line 75 of file forward_zernike_images.cpp.

76 {
77  if (!verbose)
78  return;
79  XmippMetadataProgram::show();
80  std::cout
81  << "Output directory: " << fnOutDir << std::endl
82  << "Reference volume: " << fnVolR << std::endl
83  << "Reference mask: " << fnMaskR << std::endl
84  << "Max. Shift: " << maxShift << std::endl
85  << "Max. Angular Change: " << maxAngularChange << std::endl
86  << "Max. Resolution: " << maxResol << std::endl
87  << "Sampling: " << Ts << std::endl
88  << "Max. Radius: " << Rmax << std::endl
89  << "Max. Radius Deform. " << RmaxDef << std::endl
90  << "Zernike Degree: " << L1 << std::endl
91  << "SH Degree: " << L2 << std::endl
92  << "Step: " << loop_step << std::endl
93  << "Correct CTF: " << useCTF << std::endl
94  << "Optimize alignment: " << optimizeAlignment << std::endl
95  << "Optimize deformation:" << optimizeDeformation<< std::endl
96  << "Optimize defocus: " << optimizeDefocus << std::endl
97  << "Phase flipped: " << phaseFlipped << std::endl
98  << "Regularization: " << lambda << std::endl
99  << "Blob radius: " << blob_r << std::endl
100  << "Image mode: " << image_mode << std::endl
101  ;
102 }
ProgForwardZernikeImages::fnOutDir
FileName fnOutDir
Output directory.
Definition: forward_zernike_images.h:51
XmippProgram::verbose
int verbose
Verbosity level.
Definition: xmipp_program.h:161
XmippMetadataProgram::show
void show() const override
Definition: xmipp_metadata_program.cpp:232

◆ splattingAtPos()

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

Definition at line 1174 of file forward_zernike_images.cpp.

1174  {
1175  double x_pos = r[0];
1176  double y_pos = r[1];
1177  int i0 = XMIPP_MAX(CEIL(y_pos - loop_step), STARTINGY(mV));
1178  int iF = XMIPP_MIN(FLOOR(y_pos + loop_step), FINISHINGY(mV));
1179  int j0 = XMIPP_MAX(CEIL(x_pos - loop_step), STARTINGX(mV));
1180  int jF = XMIPP_MIN(FLOOR(x_pos + loop_step), FINISHINGX(mV));
1181 
1182  double m = 1. / loop_step;
1183  for (int i = i0; i <= iF; i++)
1184  {
1185  double y_val = 1. - m * ABS(i - y_pos);
1186  for (int j = j0; j <= jF; j++)
1187  {
1188  double x_val = 1. - m * ABS(j - x_pos);
1189  A2D_ELEM(mP, i, j) += weight * x_val * y_val;
1190  }
1191  }
1192 }
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
XMIPP_MAX
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
i
#define i
Definition: numerical_recipes.cpp:2493
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
CEIL
#define CEIL(x)
Definition: xmipp_macros.h:225
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
j
#define j
Definition: numerical_recipes.cpp:2493
m
int m
Definition: image_find_center.cpp:57
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63

◆ transformImageSph()

double ProgForwardZernikeImages::transformImageSph ( double *  pclnm)

Definition at line 291 of file forward_zernike_images.cpp.

292 {
293  const MultidimArray<double> &mV=V();
294  idx_z_clnm.clear();
295  z_clnm_diff.clear();
296  FOR_ALL_ELEMENTS_IN_MATRIX1D(clnm)
297  {
298  VEC_ELEM(clnm,i)=pclnm[i+1];
299  if (VEC_ELEM(clnm,i) != VEC_ELEM(prev_clnm,i))
300  {
301  idx_z_clnm.push_back(i);
302  z_clnm_diff.push_back(VEC_ELEM(clnm, i) - VEC_ELEM(prev_clnm, i));
303  }
304  }
305  // std::cout << std::endl;
306  double deformation=0.0;
307  totalDeformation=0.0;
308 
309  P[0]().initZeros(Xdim,Xdim);
310  P[0]().setXmippOrigin();
311  double currentRot=old_rot[0] + deltaRot[0];
312  double currentTilt=old_tilt[0] + deltaTilt[0];
313  double currentPsi=old_psi[0] + deltaPsi[0];
314  deformVol(P[0](), mV, deformation, currentRot, currentTilt, currentPsi);
315 
316  double cost=0.0;
317  const MultidimArray<int> &mMask2D=mask2D;
318  double corr2 = 0.0;
319  double corr3 = 0.0;
320 
321  switch (num_images)
322  {
323  case 2:
324  {
325  P[1]().initZeros(Xdim,Xdim);
326  P[1]().setXmippOrigin();
327  currentRot = old_rot[1] + deltaRot[1];
328  currentTilt = old_tilt[1] + deltaTilt[1];
329  currentPsi = old_psi[1] + deltaPsi[1];
330  deformVol(P[1](), mV, deformation, currentRot, currentTilt, currentPsi);
331 
332  if (old_flip)
333  {
334  MAT_ELEM(A1, 0, 0) *= -1;
335  MAT_ELEM(A1, 0, 1) *= -1;
336  MAT_ELEM(A1, 0, 2) *= -1;
337  MAT_ELEM(A2, 0, 0) *= -1;
338  MAT_ELEM(A2, 0, 1) *= -1;
339  MAT_ELEM(A2, 0, 2) *= -1;
340  }
341  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[0](), Ifiltered[0](), A1,
342  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
343  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[1](), Ifiltered[1](), A2,
344  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
345  // filter.applyMaskSpace(P[1]());
346  const MultidimArray<double> mP2 = P[1]();
347  MultidimArray<double> &mI2filteredp = Ifilteredp[1]();
348  corr2 = correlationIndex(mI2filteredp, mP2, &mMask2D);
349  }
350  break;
351 
352  case 3:
353  {
354  P[1]().initZeros(Xdim,Xdim);
355  P[1]().setXmippOrigin();
356  currentRot = old_rot[1] + deltaRot[1];
357  currentTilt = old_tilt[1] + deltaTilt[1];
358  currentPsi = old_psi[1] + deltaPsi[1];
359  deformVol(P[1](), mV, deformation, currentRot, currentTilt, currentPsi);
360 
361  P[2]().initZeros(Xdim,Xdim);
362  P[2]().setXmippOrigin();
363  currentRot = old_rot[2] + deltaRot[2];
364  currentTilt = old_tilt[2] + deltaTilt[2];
365  currentPsi = old_psi[2] + deltaPsi[2];
366  deformVol(P[2](), mV, deformation, currentRot, currentTilt, currentPsi);
367 
368  if (old_flip)
369  {
370  MAT_ELEM(A1, 0, 0) *= -1;
371  MAT_ELEM(A1, 0, 1) *= -1;
372  MAT_ELEM(A1, 0, 2) *= -1;
373  MAT_ELEM(A2, 0, 0) *= -1;
374  MAT_ELEM(A2, 0, 1) *= -1;
375  MAT_ELEM(A2, 0, 2) *= -1;
376  MAT_ELEM(A3, 0, 0) *= -1;
377  MAT_ELEM(A3, 0, 1) *= -1;
378  MAT_ELEM(A3, 0, 2) *= -1;
379  }
380  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[0](), Ifiltered[0](), A1,
381  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
382  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[1](), Ifiltered[1](), A2,
383  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
384  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[2](), Ifiltered[2](), A3,
385  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
386  // filter.applyMaskSpace(P[1]());
387  // filter.applyMaskSpace(P[2]());
388  const MultidimArray<double> mP2 = P[1]();
389  const MultidimArray<double> mP3 = P[2]();
390  MultidimArray<double> &mI2filteredp = Ifilteredp[1]();
391  MultidimArray<double> &mI3filteredp = Ifilteredp[2]();
392  corr2 = correlationIndex(mI2filteredp, mP2, &mMask2D);
393  corr3 = correlationIndex(mI3filteredp, mP3, &mMask2D);
394  }
395  break;
396 
397  default:
398  if (old_flip)
399  {
400  MAT_ELEM(A1, 0, 0) *= -1;
401  MAT_ELEM(A1, 0, 1) *= -1;
402  MAT_ELEM(A1, 0, 2) *= -1;
403  }
404  applyGeometry(xmipp_transformation::LINEAR,Ifilteredp[0](),Ifiltered[0](),A1,
405  xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP,0.);
406  break;
407  }
408 
409  filter.generateMask(P[0]());
410  filter.applyMaskSpace(P[0]());
411  if (hasCTF)
412  {
413  double defocusU=old_defocusU[0]+deltaDefocusU[0];
414  double defocusV=old_defocusV[0]+deltaDefocusV[0];
415  double angle=old_defocusAngle[0]+deltaDefocusAngle[0];
416  if (defocusU!=currentDefocusU[0] || defocusV!=currentDefocusV[0] || angle!=currentAngle[0]) {
417  updateCTFImage(defocusU,defocusV,angle);
418  }
419  // FilterCTF1.ctf = ctf;
420  FilterCTF1.generateMask(P[0]());
421  if (phaseFlipped)
422  FilterCTF1.correctPhase();
423  FilterCTF1.applyMaskSpace(P[0]());
424  }
425 
426 
427  const MultidimArray<double> mP1=P[0]();
428  MultidimArray<double> &mI1filteredp=Ifilteredp[0]();
429  double corr1=correlationIndex(mI1filteredp,mP1,&mMask2D);
430 
431  switch (num_images)
432  {
433  case 2:
434  cost=-(corr1+corr2+corr3) / 2.0;
435  break;
436  case 3:
437  cost=-(corr1+corr2+corr3) / 3.0;
438  break;
439  default:
440  cost=-(corr1+corr2+corr3);
441  break;
442  }
443 
444 #ifdef DEBUG
445  std::cout << "A=" << A << std::endl;
446  Image<double> save;
447  save()=P();
448  save.write("PPPtheo.xmp");
449  save()=Ifilteredp();
450  save.write("PPPfilteredp.xmp");
451  save()=Ifiltered();
452  save.write("PPPfiltered.xmp");
453  // Vdeformed.write("PPPVdeformed.vol");
454  std::cout << "Cost=" << cost << " deformation=" << deformation << std::endl;
455  std::cout << "Press any key" << std::endl;
456  char c; std::cin >> c;
457 #endif
458 
459  if (showOptimization)
460  {
461  std::cout << "A1=" << A1 << std::endl;
462  Image<double> save;
463  save()=P[0]();
464  save.write("PPPtheo1.xmp");
465  save()=Ifilteredp[0]();
466  save.write("PPPfilteredp1.xmp");
467  save()=Ifiltered[0]();
468  save.write("PPPfiltered1.xmp");
469 
470  switch (num_images)
471  {
472  case 2:
473  {
474  std::cout << "A2=" << A2 << std::endl;
475  save()=P[1]();
476  save.write("PPPtheo2.xmp");
477  save()=Ifilteredp[1]();
478  save.write("PPPfilteredp2.xmp");
479  save()=Ifiltered[1]();
480  save.write("PPPfiltered2.xmp");
481  }
482  break;
483 
484  case 3:
485  {
486  std::cout << "A2=" << A2 << std::endl;
487  save()=P[1]();
488  save.write("PPPtheo2.xmp");
489  save()=Ifilteredp[1]();
490  save.write("PPPfilteredp2.xmp");
491  save()=Ifiltered[1]();
492  save.write("PPPfiltered2.xmp");
493 
494  std::cout << "A3=" << A3 << std::endl;
495  save()=P[2]();
496  save.write("PPPtheo3.xmp");
497  save()=Ifilteredp[2]();
498  save.write("PPPfilteredp3.xmp");
499  save()=Ifiltered[2]();
500  save.write("PPPfiltered3.xmp");
501  }
502  break;
503 
504  default:
505  break;
506  }
507  Vdeformed.write("PPPVdeformed.vol");
508  std::cout << "Deformation=" << totalDeformation << std::endl;
509  std::cout << "Press any key" << std::endl;
510  char c; std::cin >> c;
511  }
512 
513  prev_clnm = clnm;
514  double retval=cost+lambda*abs(deformation - prior_deformation);
515  if (showOptimization)
516  std::cout << cost << " " << deformation << " " << lambda*deformation << " " << sumV << " " << retval << std::endl;
517  return retval;
518 }
ProgForwardZernikeImages::Ifiltered
std::vector< Image< double > > Ifiltered
Definition: forward_zernike_images.h:105
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
ProgForwardZernikeImages::old_psi
std::vector< double > old_psi
Definition: forward_zernike_images.h:114
ProgForwardZernikeImages::old_tilt
std::vector< double > old_tilt
Definition: forward_zernike_images.h:114
c
doublereal * c
Definition: numerical_recipes.cpp:2230
ProgForwardZernikeImages::old_defocusV
std::vector< double > old_defocusV
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::deltaPsi
std::vector< double > deltaPsi
Definition: forward_zernike_images.h:114
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
ProgForwardZernikeImages::deltaTilt
std::vector< double > deltaTilt
Definition: forward_zernike_images.h:114
applyGeometry
void applyGeometry(int SplineDegree, MultidimArray< std::complex< double > > &V2, const MultidimArray< std::complex< double > > &V1, const Matrix2D< double > &A, bool inv, bool wrap, std::complex< double > outside, MultidimArray< double > *BcoeffsPtr)
Definition: transformations.cpp:668
ProgForwardZernikeImages::old_defocusU
std::vector< double > old_defocusU
Definition: forward_zernike_images.h:122
ImageBase::write
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
Definition: xmipp_image_base.cpp:292
correlationIndex
double correlationIndex(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)
Definition: multidim_array.h:3974
ProgForwardZernikeImages::Ifilteredp
std::vector< Image< double > > Ifilteredp
Definition: forward_zernike_images.h:106
ProgForwardZernikeImages::deltaDefocusAngle
std::vector< double > deltaDefocusAngle
Definition: forward_zernike_images.h:122
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
ProgForwardZernikeImages::deltaDefocusU
std::vector< double > deltaDefocusU
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::updateCTFImage
void updateCTFImage(double defocusU, double defocusV, double angle)
Definition: forward_zernike_images.cpp:1014
i
#define i
Definition: numerical_recipes.cpp:2493
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
FOR_ALL_ELEMENTS_IN_MATRIX1D
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
ProgForwardZernikeImages::deformVol
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.
Definition: forward_zernike_images.cpp:1047
ProgForwardZernikeImages::P
std::vector< Image< double > > P
Definition: forward_zernike_images.h:108
ProgForwardZernikeImages::currentAngle
std::vector< double > currentAngle
Definition: forward_zernike_images.h:124
ProgForwardZernikeImages::old_rot
std::vector< double > old_rot
Definition: forward_zernike_images.h:114
ProgForwardZernikeImages::old_defocusAngle
std::vector< double > old_defocusAngle
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::z_clnm_diff
std::vector< double > z_clnm_diff
Definition: forward_zernike_images.h:151
ProgForwardZernikeImages::idx_z_clnm
std::vector< size_t > idx_z_clnm
Definition: forward_zernike_images.h:150
ProgForwardZernikeImages::currentDefocusU
std::vector< double > currentDefocusU
Definition: forward_zernike_images.h:124
FourierFilter::generateMask
void generateMask(MultidimArray< double > &v)
Definition: fourier_filter.cpp:605
ProgForwardZernikeImages::deltaDefocusV
std::vector< double > deltaDefocusV
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::currentDefocusV
std::vector< double > currentDefocusV
Definition: forward_zernike_images.h:124
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710
ProgForwardZernikeImages::deltaRot
std::vector< double > deltaRot
Definition: forward_zernike_images.h:114

◆ updateCTFImage()

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

Definition at line 1014 of file forward_zernike_images.cpp.

1015 {
1016  FilterCTF1.ctf.K=1; // get pure CTF with no envelope
1017  currentDefocusU[0]=FilterCTF1.ctf.DeltafU=defocusU;
1018  currentDefocusV[0]=FilterCTF1.ctf.DeltafV=defocusV;
1019  currentAngle[0]=FilterCTF1.ctf.azimuthal_angle=angle;
1020  FilterCTF1.ctf.produceSideInfo();
1021 }
FourierFilter::ctf
CTFDescription ctf
Definition: fourier_filter.h:134
CTFDescription::DeltafV
double DeltafV
Defocus in V (in Angstroms). Negative values are underfocused.
Definition: ctf.h:830
CTFDescription::DeltafU
double DeltafU
Global gain. By default, 1.
Definition: ctf.h:828
CTFDescription::azimuthal_angle
double azimuthal_angle
Azimuthal angle (between X and U) in degrees.
Definition: ctf.h:832
ProgForwardZernikeImages::currentAngle
std::vector< double > currentAngle
Definition: forward_zernike_images.h:124
CTFDescription1D::K
double K
Global gain. By default, 1.
Definition: ctf.h:238
CTFDescription::produceSideInfo
void produceSideInfo()
Produce Side information.
Definition: ctf.cpp:1392
ProgForwardZernikeImages::currentDefocusU
std::vector< double > currentDefocusU
Definition: forward_zernike_images.h:124
ProgForwardZernikeImages::currentDefocusV
std::vector< double > currentDefocusV
Definition: forward_zernike_images.h:124

◆ weightsInterpolation3D()

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

Definition at line 1143 of file forward_zernike_images.cpp.

1143  {
1144  Matrix1D<double> w;
1145  w.initZeros(8);
1146 
1147  int x0 = FLOOR(x);
1148  double fx0 = x - x0;
1149  int x1 = x0 + 1;
1150  double fx1 = x1 - x;
1151 
1152  int y0 = FLOOR(y);
1153  double fy0 = y - y0;
1154  int y1 = y0 + 1;
1155  double fy1 = y1 - y;
1156 
1157  int z0 = FLOOR(z);
1158  double fz0 = z - z0;
1159  int z1 = z0 + 1;
1160  double fz1 = z1 - z;
1161 
1162  VEC_ELEM(w,0) = fx1 * fy1 * fz1; // w000 (x0,y0,z0)
1163  VEC_ELEM(w,1) = fx1 * fy1 * fz0; // w001 (x0,y0,z1)
1164  VEC_ELEM(w,2) = fx1 * fy0 * fz1; // w010 (x0,y1,z0)
1165  VEC_ELEM(w,3) = fx1 * fy0 * fz0; // w011 (x0,y1,z1)
1166  VEC_ELEM(w,4) = fx0 * fy1 * fz1; // w100 (x1,y0,z0)
1167  VEC_ELEM(w,5) = fx0 * fy1 * fz0; // w101 (x1,y0,z1)
1168  VEC_ELEM(w,6) = fx0 * fy0 * fz1; // w110 (x1,y1,z0)
1169  VEC_ELEM(w,7) = fx0 * fy0 * fz0; // w111 (x1,y1,z1)
1170 
1171  return w;
1172 }
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
y
static double * y
Definition: numerical_recipes.cpp:8487
w
doublereal * w
Definition: numerical_recipes.cpp:2477
z0
#define z0
x
doublereal * x
Definition: numerical_recipes.cpp:2230
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
y0
#define y0
x0
#define x0
z
double z
Definition: numerical_recipes.cpp:8490
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592

◆ writeImageParameters()

void ProgForwardZernikeImages::writeImageParameters ( MDRow &  row)
virtual

Write the parameters found for one image

Definition at line 873 of file forward_zernike_images.cpp.

873  {
874  int pos = 3*vecSize;
875  if (flagEnabled==1) {
876  row.setValue(MDL_ENABLED, 1);
877  }
878  else {
879  row.setValue(MDL_ENABLED, -1);
880  }
881 
882  switch (num_images)
883  {
884  case 2:
885  row.setValue(MDL_ANGLE_ROT, old_rot[0]+p(pos+4));
886  row.setValue(MDL_ANGLE_ROT2, old_rot[1]+p(pos+5));
887  row.setValue(MDL_ANGLE_TILT, old_tilt[0]+p(pos+6));
888  row.setValue(MDL_ANGLE_TILT2, old_tilt[1]+p(pos+7));
889  row.setValue(MDL_ANGLE_PSI, old_psi[0]+p(pos+8));
890  row.setValue(MDL_ANGLE_PSI2, old_psi[1]+p(pos+9));
891  row.setValue(MDL_SHIFT_X, old_shiftX[0]+p(pos));
892  row.setValue(MDL_SHIFT_X2, old_shiftX[1]+p(pos+1));
893  row.setValue(MDL_SHIFT_Y, old_shiftY[0]+p(pos+2));
894  row.setValue(MDL_SHIFT_Y2, old_shiftY[1]+p(pos+3));
895  break;
896 
897  case 3:
898  row.setValue(MDL_ANGLE_ROT, old_rot[0]+p(pos+6));
899  row.setValue(MDL_ANGLE_ROT2, old_rot[1]+p(pos+7));
900  row.setValue(MDL_ANGLE_ROT3, old_rot[2]+p(pos+8));
901  row.setValue(MDL_ANGLE_TILT, old_tilt[0]+p(pos+9));
902  row.setValue(MDL_ANGLE_TILT2, old_tilt[1]+p(pos+10));
903  row.setValue(MDL_ANGLE_TILT3, old_tilt[2]+p(pos+11));
904  row.setValue(MDL_ANGLE_PSI, old_psi[0]+p(pos+12));
905  row.setValue(MDL_ANGLE_PSI2, old_psi[1]+p(pos+13));
906  row.setValue(MDL_ANGLE_PSI3, old_psi[2]+p(pos+14));
907  row.setValue(MDL_SHIFT_X, old_shiftX[0]+p(pos));
908  row.setValue(MDL_SHIFT_X2, old_shiftX[1]+p(pos+1));
909  row.setValue(MDL_SHIFT_X3, old_shiftX[2]+p(pos+2));
910  row.setValue(MDL_SHIFT_Y, old_shiftY[0]+p(pos+3));
911  row.setValue(MDL_SHIFT_Y2, old_shiftY[1]+p(pos+4));
912  row.setValue(MDL_SHIFT_Y3, old_shiftY[2]+p(pos+5));
913  break;
914 
915  default:
916  row.setValue(MDL_ANGLE_ROT, old_rot[0]+p(pos+2));
917  row.setValue(MDL_ANGLE_TILT, old_tilt[0]+p(pos+3));
918  row.setValue(MDL_ANGLE_PSI, old_psi[0]+p(pos+4));
919  row.setValue(MDL_SHIFT_X, old_shiftX[0]+p(pos));
920  row.setValue(MDL_SHIFT_Y, old_shiftY[0]+p(pos+1));
921  break;
922  }
923 
924  row.setValue(MDL_SPH_DEFORMATION, totalDeformation);
925  std::vector<double> vectortemp;
926  size_t end_clnm = VEC_XSIZE(clnm)-algn_params-ctf_params;
927  for (int j = 0; j < end_clnm; j++) {
928  vectortemp.push_back(clnm(j));
929  }
930  row.setValue(MDL_SPH_COEFFICIENTS, vectortemp);
931  row.setValue(MDL_COST, correlation);
932 }
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
ProgForwardZernikeImages::old_shiftY
std::vector< double > old_shiftY
Definition: forward_zernike_images.h:116
MDL_ANGLE_TILT2
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:60
MDL_SHIFT_X3
Shift for the image in the X axis (double)
Definition: metadata_label.h:433
MDL_SHIFT_Y2
Shift for the image in the Y axis (double)
Definition: metadata_label.h:436
ProgForwardZernikeImages::old_psi
std::vector< double > old_psi
Definition: forward_zernike_images.h:114
ProgForwardZernikeImages::old_tilt
std::vector< double > old_tilt
Definition: forward_zernike_images.h:114
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
MDL_SHIFT_X
Shift for the image in the X axis (double)
Definition: metadata_label.h:431
MDL_ANGLE_TILT3
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:61
MDL_ANGLE_PSI
Special label to be used when gathering MDs in MpiMetadataPrograms.
Definition: metadata_label.h:51
MDL_ENABLED
Is this image enabled? (int [-1 or 1])
Definition: metadata_label.h:231
MDL_SHIFT_Y3
Shift for the image in the Y axis (double)
Definition: metadata_label.h:437
MDL_ANGLE_ROT2
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:56
MDL_SPH_DEFORMATION
Deformation in voxels.
Definition: metadata_label.h:449
MDL_COST
Cost for the image (double)
Definition: metadata_label.h:90
ProgForwardZernikeImages::old_rot
std::vector< double > old_rot
Definition: forward_zernike_images.h:114
MDL_ANGLE_PSI2
Psi angle of an image (double,degrees)
Definition: metadata_label.h:52
j
#define j
Definition: numerical_recipes.cpp:2493
MDRow::setValue
void setValue(MDLabel label, const T &d, bool addLabel=true)
Definition: metadata_row_base.h:151
MDL_SHIFT_X2
Shift for the image in the X axis (double)
Definition: metadata_label.h:432
ProgForwardZernikeImages::old_shiftX
std::vector< double > old_shiftX
Definition: forward_zernike_images.h:116
MDL_SPH_COEFFICIENTS
Deformation coefficients.
Definition: metadata_label.h:448
MDL_ANGLE_ROT3
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:57
MDL_SHIFT_Y
Shift for the image in the Y axis (double)
Definition: metadata_label.h:435
MDL_ANGLE_PSI3
Psi angle of an image (double,degrees)
Definition: metadata_label.h:53

Member Data Documentation

◆ A1

Matrix2D<double> ProgForwardZernikeImages::A1

Definition at line 112 of file forward_zernike_images.h.

◆ A2

Matrix2D<double> ProgForwardZernikeImages::A2

Definition at line 112 of file forward_zernike_images.h.

◆ A3

Matrix2D<double> ProgForwardZernikeImages::A3

Definition at line 112 of file forward_zernike_images.h.

◆ algn_params

int ProgForwardZernikeImages::algn_params

Definition at line 85 of file forward_zernike_images.h.

◆ blob

struct blobtype ProgForwardZernikeImages::blob

Definition at line 145 of file forward_zernike_images.h.

◆ blob_r

double ProgForwardZernikeImages::blob_r

Definition at line 146 of file forward_zernike_images.h.

◆ clnm

Matrix1D<double> ProgForwardZernikeImages::clnm

Definition at line 132 of file forward_zernike_images.h.

◆ correlation

double ProgForwardZernikeImages::correlation

Definition at line 141 of file forward_zernike_images.h.

◆ ctf_params

int ProgForwardZernikeImages::ctf_params

Definition at line 87 of file forward_zernike_images.h.

◆ currentAngle

std::vector<double> ProgForwardZernikeImages::currentAngle

Definition at line 124 of file forward_zernike_images.h.

◆ currentDefocusU

std::vector<double> ProgForwardZernikeImages::currentDefocusU

Definition at line 124 of file forward_zernike_images.h.

◆ currentDefocusV

std::vector<double> ProgForwardZernikeImages::currentDefocusV

Definition at line 124 of file forward_zernike_images.h.

◆ deltaDefocusAngle

std::vector<double> ProgForwardZernikeImages::deltaDefocusAngle

Definition at line 122 of file forward_zernike_images.h.

◆ deltaDefocusU

std::vector<double> ProgForwardZernikeImages::deltaDefocusU

Definition at line 122 of file forward_zernike_images.h.

◆ deltaDefocusV

std::vector<double> ProgForwardZernikeImages::deltaDefocusV

Definition at line 122 of file forward_zernike_images.h.

◆ deltaPsi

std::vector<double> ProgForwardZernikeImages::deltaPsi

Definition at line 114 of file forward_zernike_images.h.

◆ deltaRot

std::vector<double> ProgForwardZernikeImages::deltaRot

Definition at line 114 of file forward_zernike_images.h.

◆ deltaTilt

std::vector<double> ProgForwardZernikeImages::deltaTilt

Definition at line 114 of file forward_zernike_images.h.

◆ deltaX

std::vector<double> ProgForwardZernikeImages::deltaX

Definition at line 116 of file forward_zernike_images.h.

◆ deltaY

std::vector<double> ProgForwardZernikeImages::deltaY

Definition at line 116 of file forward_zernike_images.h.

◆ df

MultidimArray<double> ProgForwardZernikeImages::df

Definition at line 149 of file forward_zernike_images.h.

◆ filter

FourierFilter ProgForwardZernikeImages::filter

Definition at line 110 of file forward_zernike_images.h.

◆ FilterCTF1

FourierFilter ProgForwardZernikeImages::FilterCTF1

Definition at line 126 of file forward_zernike_images.h.

◆ FilterCTF2

FourierFilter ProgForwardZernikeImages::FilterCTF2

Definition at line 127 of file forward_zernike_images.h.

◆ FilterCTF3

FourierFilter ProgForwardZernikeImages::FilterCTF3

Definition at line 128 of file forward_zernike_images.h.

◆ filterp

FourierFilter ProgForwardZernikeImages::filterp

Definition at line 110 of file forward_zernike_images.h.

◆ flagEnabled

int ProgForwardZernikeImages::flagEnabled

Definition at line 89 of file forward_zernike_images.h.

◆ fnImage

std::vector<FileName> ProgForwardZernikeImages::fnImage

Definition at line 101 of file forward_zernike_images.h.

◆ fnMaskR

FileName ProgForwardZernikeImages::fnMaskR

Filename of the reference volume mask

Definition at line 49 of file forward_zernike_images.h.

◆ fnOutDir

FileName ProgForwardZernikeImages::fnOutDir

Output directory.

Definition at line 51 of file forward_zernike_images.h.

◆ fnVolR

FileName ProgForwardZernikeImages::fnVolR

Filename of the reference volume

Definition at line 47 of file forward_zernike_images.h.

◆ hasCTF

bool ProgForwardZernikeImages::hasCTF

Definition at line 120 of file forward_zernike_images.h.

◆ I

std::vector<Image<double> > ProgForwardZernikeImages::I

Definition at line 104 of file forward_zernike_images.h.

◆ idx_z_clnm

std::vector<size_t> ProgForwardZernikeImages::idx_z_clnm

Definition at line 150 of file forward_zernike_images.h.

◆ Ifiltered

std::vector<Image<double> > ProgForwardZernikeImages::Ifiltered

Definition at line 105 of file forward_zernike_images.h.

◆ Ifilteredp

std::vector<Image<double> > ProgForwardZernikeImages::Ifilteredp

Definition at line 106 of file forward_zernike_images.h.

◆ ignoreCTF

bool ProgForwardZernikeImages::ignoreCTF

Definition at line 73 of file forward_zernike_images.h.

◆ image_mode

int ProgForwardZernikeImages::image_mode

Definition at line 90 of file forward_zernike_images.h.

◆ L1

int ProgForwardZernikeImages::L1

Degrees of Zernike polynomials and spherical harmonics

Definition at line 53 of file forward_zernike_images.h.

◆ L2

int ProgForwardZernikeImages::L2

Definition at line 53 of file forward_zernike_images.h.

◆ lambda

double ProgForwardZernikeImages::lambda

Definition at line 79 of file forward_zernike_images.h.

◆ loop_step

int ProgForwardZernikeImages::loop_step

Definition at line 143 of file forward_zernike_images.h.

◆ mask2D

MultidimArray<int> ProgForwardZernikeImages::mask2D

Definition at line 97 of file forward_zernike_images.h.

◆ maxAngularChange

double ProgForwardZernikeImages::maxAngularChange

Maximum angular change allowed

Definition at line 59 of file forward_zernike_images.h.

◆ maxResol

double ProgForwardZernikeImages::maxResol

Maximum frequency (A)

Definition at line 61 of file forward_zernike_images.h.

◆ maxShift

double ProgForwardZernikeImages::maxShift

Maximum shift allowed

Definition at line 57 of file forward_zernike_images.h.

◆ num_images

int ProgForwardZernikeImages::num_images

Definition at line 83 of file forward_zernike_images.h.

◆ old_defocusAngle

std::vector<double> ProgForwardZernikeImages::old_defocusAngle

Definition at line 122 of file forward_zernike_images.h.

◆ old_defocusU

std::vector<double> ProgForwardZernikeImages::old_defocusU

Definition at line 122 of file forward_zernike_images.h.

◆ old_defocusV

std::vector<double> ProgForwardZernikeImages::old_defocusV

Definition at line 122 of file forward_zernike_images.h.

◆ old_flip

bool ProgForwardZernikeImages::old_flip

Definition at line 118 of file forward_zernike_images.h.

◆ old_psi

std::vector<double> ProgForwardZernikeImages::old_psi

Definition at line 114 of file forward_zernike_images.h.

◆ old_rot

std::vector<double> ProgForwardZernikeImages::old_rot

Definition at line 114 of file forward_zernike_images.h.

◆ old_shiftX

std::vector<double> ProgForwardZernikeImages::old_shiftX

Definition at line 116 of file forward_zernike_images.h.

◆ old_shiftY

std::vector<double> ProgForwardZernikeImages::old_shiftY

Definition at line 116 of file forward_zernike_images.h.

◆ old_tilt

std::vector<double> ProgForwardZernikeImages::old_tilt

Definition at line 114 of file forward_zernike_images.h.

◆ optimizeAlignment

bool ProgForwardZernikeImages::optimizeAlignment

Definition at line 67 of file forward_zernike_images.h.

◆ optimizeDefocus

bool ProgForwardZernikeImages::optimizeDefocus

Definition at line 71 of file forward_zernike_images.h.

◆ optimizeDeformation

bool ProgForwardZernikeImages::optimizeDeformation

Definition at line 69 of file forward_zernike_images.h.

◆ optimizeRadius

bool ProgForwardZernikeImages::optimizeRadius

Definition at line 75 of file forward_zernike_images.h.

◆ p

Matrix1D<double> ProgForwardZernikeImages::p

Definition at line 88 of file forward_zernike_images.h.

◆ P

std::vector<Image<double> > ProgForwardZernikeImages::P

Definition at line 108 of file forward_zernike_images.h.

◆ phaseFlipped

bool ProgForwardZernikeImages::phaseFlipped

Definition at line 77 of file forward_zernike_images.h.

◆ prev_clnm

Matrix1D<double> ProgForwardZernikeImages::prev_clnm

Definition at line 132 of file forward_zernike_images.h.

◆ prior_deformation

double ProgForwardZernikeImages::prior_deformation

Definition at line 136 of file forward_zernike_images.h.

◆ R

Matrix2D<double> ProgForwardZernikeImages::R

Definition at line 152 of file forward_zernike_images.h.

◆ resume

bool ProgForwardZernikeImages::resume

Resume computations

Definition at line 95 of file forward_zernike_images.h.

◆ Rmax

int ProgForwardZernikeImages::Rmax

Maximum radius

Definition at line 65 of file forward_zernike_images.h.

◆ RmaxDef

int ProgForwardZernikeImages::RmaxDef

Definition at line 81 of file forward_zernike_images.h.

◆ showOptimization

bool ProgForwardZernikeImages::showOptimization

Definition at line 139 of file forward_zernike_images.h.

◆ steps_cp

Matrix1D<double> ProgForwardZernikeImages::steps_cp

Definition at line 134 of file forward_zernike_images.h.

◆ sumV

int ProgForwardZernikeImages::sumV

Definition at line 137 of file forward_zernike_images.h.

◆ totalDeformation

double ProgForwardZernikeImages::totalDeformation

Definition at line 136 of file forward_zernike_images.h.

◆ Ts

double ProgForwardZernikeImages::Ts

Sampling rate

Definition at line 63 of file forward_zernike_images.h.

◆ useCTF

bool ProgForwardZernikeImages::useCTF

Definition at line 91 of file forward_zernike_images.h.

◆ V

Image<double> ProgForwardZernikeImages::V

Definition at line 103 of file forward_zernike_images.h.

◆ V_mask

MultidimArray<int> ProgForwardZernikeImages::V_mask

Definition at line 97 of file forward_zernike_images.h.

◆ Vdeformed

Image<double> ProgForwardZernikeImages::Vdeformed

Definition at line 103 of file forward_zernike_images.h.

◆ vecSize

int ProgForwardZernikeImages::vecSize

Definition at line 130 of file forward_zernike_images.h.

◆ vL1

Matrix1D<int> ProgForwardZernikeImages::vL1

Zernike and SPH coefficients vectors

Definition at line 55 of file forward_zernike_images.h.

◆ vL2

Matrix1D<int> ProgForwardZernikeImages::vL2

Definition at line 55 of file forward_zernike_images.h.

◆ vM

Matrix1D<int> ProgForwardZernikeImages::vM

Definition at line 55 of file forward_zernike_images.h.

◆ vN

Matrix1D<int> ProgForwardZernikeImages::vN

Definition at line 55 of file forward_zernike_images.h.

◆ vpos

MultidimArray<double> ProgForwardZernikeImages::vpos

Definition at line 149 of file forward_zernike_images.h.

◆ Xdim

size_t ProgForwardZernikeImages::Xdim

Definition at line 99 of file forward_zernike_images.h.

◆ z_clnm_diff

std::vector<double> ProgForwardZernikeImages::z_clnm_diff

Definition at line 151 of file forward_zernike_images.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.13