Xmipp  v3.23.11-Nereus
forward_zernike_images.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Carlos Oscar Sanchez Sorzano (coss@cnb.csic.es)
4  * David Herreros Calero (dherreros@cnb.csic.es)
5  *
6  * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
21  * 02111-1307 USA
22  *
23  * All comments concerning this program package may be sent to the
24  * e-mail address 'xmipp@cnb.csic.es'
25  ***************************************************************************/
26 
27 #include "forward_zernike_images.h"
28 #include "core/transformations.h"
29 #include "core/xmipp_image_extension.h"
30 #include "core/xmipp_image_generic.h"
31 #include "data/projection.h"
32 #include "data/mask.h"
33 
34 // Empty constructor =======================================================
35 ProgForwardZernikeImages::ProgForwardZernikeImages() : Rerunable("")
36 {
37  resume = false;
38  produces_a_metadata = true;
39  each_image_produces_an_output = false;
40  showOptimization = false;
41 }
42 
43 ProgForwardZernikeImages::~ProgForwardZernikeImages() = default;
44 
45 // Read arguments ==========================================================
46 void ProgForwardZernikeImages::readParams()
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 }
73 
74 // Show ====================================================================
75 void ProgForwardZernikeImages::show()
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 }
103 
104 // usage ===================================================================
105 void ProgForwardZernikeImages::defineParams()
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 }
137 
138 void ProgForwardZernikeImages::preProcess()
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 }
284 
285 void ProgForwardZernikeImages::finishProcessing() {
286  XmippMetadataProgram::finishProcessing();
287  rename(Rerunable::getFileName().c_str(), (fnOutDir + fn_out).c_str());
288 }
289 
290 // #define DEBUG
291 double ProgForwardZernikeImages::transformImageSph(double *pclnm)
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 }
519 
520 double continuousZernikeCost(double *x, void *_prm)
521 {
522  ProgForwardZernikeImages *prm=(ProgForwardZernikeImages *)_prm;
523  int idx = 3*(prm->vecSize);
524  // TODO: Optimize parameters for each image (not sharing)
525  // prm->deltaDefocusU[0]=x[idx+6]; prm->deltaDefocusU[1]=x[idx+6];
526  // prm->deltaDefocusV[0]=x[idx+7]; prm->deltaDefocusV[1]=x[idx+7];
527  // prm->deltaDefocusAngle[0]=x[idx+8]; prm->deltaDefocusAngle[1]=x[idx+8];
528 
529  switch (prm->num_images)
530  {
531  case 2:
532  prm->deltaX[0] = x[idx + 1];
533  prm->deltaY[0] = x[idx + 3];
534  prm->deltaRot[0] = x[idx + 5];
535  prm->deltaTilt[0] = x[idx + 7];
536  prm->deltaPsi[0] = x[idx + 9];
537  // prm->deltaDefocusU[0]=x[idx + 11];
538  // prm->deltaDefocusV[0]=x[idx + 13];
539  // prm->deltaDefocusAngle[0]=x[idx + 15];
540 
541  prm->deltaX[1] = x[idx + 2];
542  prm->deltaY[1] = x[idx + 4];
543  prm->deltaRot[1] = x[idx + 6];
544  prm->deltaTilt[1] = x[idx + 8];
545  prm->deltaPsi[1] = x[idx + 10];
546  // prm->deltaDefocusU[1]=x[idx + 12];
547  // prm->deltaDefocusV[1]=x[idx + 14];
548  // prm->deltaDefocusAngle[1]=x[idx + 16];
549 
550  MAT_ELEM(prm->A1, 0, 2) = prm->old_shiftX[0] + prm->deltaX[0];
551  MAT_ELEM(prm->A1, 1, 2) = prm->old_shiftY[0] + prm->deltaY[0];
552  MAT_ELEM(prm->A1, 0, 0) = 1;
553  MAT_ELEM(prm->A1, 0, 1) = 0;
554  MAT_ELEM(prm->A1, 1, 0) = 0;
555  MAT_ELEM(prm->A1, 1, 1) = 1;
556 
557  MAT_ELEM(prm->A2, 0, 2) = prm->old_shiftX[1] + prm->deltaX[1];
558  MAT_ELEM(prm->A2, 1, 2) = prm->old_shiftY[1] + prm->deltaY[1];
559  MAT_ELEM(prm->A2, 0, 0) = 1;
560  MAT_ELEM(prm->A2, 0, 1) = 0;
561  MAT_ELEM(prm->A2, 1, 0) = 0;
562  MAT_ELEM(prm->A2, 1, 1) = 1;
563  break;
564 
565  case 3:
566  prm->deltaX[0] = x[idx + 1];
567  prm->deltaY[0] = x[idx + 4];
568  prm->deltaRot[0] = x[idx + 7];
569  prm->deltaTilt[0] = x[idx + 10];
570  prm->deltaPsi[0] = x[idx + 13];
571  // prm->deltaDefocusU[0]=x[idx + 16];
572  // prm->deltaDefocusV[0]=x[idx + 19];
573  // prm->deltaDefocusAngle[0]=x[idx + 22];
574 
575  prm->deltaX[1] = x[idx + 2];
576  prm->deltaY[1] = x[idx + 5];
577  prm->deltaRot[1] = x[idx + 8];
578  prm->deltaTilt[1] = x[idx + 11];
579  prm->deltaPsi[1] = x[idx + 14];
580  // prm->deltaDefocusU[1]=x[idx + 17];
581  // prm->deltaDefocusV[1]=x[idx + 20];
582  // prm->deltaDefocusAngle[1]=x[idx + 23];
583 
584  prm->deltaX[2] = x[idx + 3];
585  prm->deltaY[2] = x[idx + 6];
586  prm->deltaRot[2] = x[idx + 9];
587  prm->deltaTilt[2] = x[idx + 12];
588  prm->deltaPsi[2] = x[idx + 15];
589  // prm->deltaDefocusU[2]=x[idx + 18];
590  // prm->deltaDefocusV[2]=x[idx + 21];
591  // prm->deltaDefocusAngle[2]=x[idx + 24];
592 
593  MAT_ELEM(prm->A1, 0, 2) = prm->old_shiftX[0] + prm->deltaX[0];
594  MAT_ELEM(prm->A1, 1, 2) = prm->old_shiftY[0] + prm->deltaY[0];
595  MAT_ELEM(prm->A1, 0, 0) = 1;
596  MAT_ELEM(prm->A1, 0, 1) = 0;
597  MAT_ELEM(prm->A1, 1, 0) = 0;
598  MAT_ELEM(prm->A1, 1, 1) = 1;
599 
600  MAT_ELEM(prm->A2, 0, 2) = prm->old_shiftX[1] + prm->deltaX[1];
601  MAT_ELEM(prm->A2, 1, 2) = prm->old_shiftY[1] + prm->deltaY[1];
602  MAT_ELEM(prm->A2, 0, 0) = 1;
603  MAT_ELEM(prm->A2, 0, 1) = 0;
604  MAT_ELEM(prm->A2, 1, 0) = 0;
605  MAT_ELEM(prm->A2, 1, 1) = 1;
606 
607  MAT_ELEM(prm->A3, 0, 2) = prm->old_shiftX[2] + prm->deltaX[2];
608  MAT_ELEM(prm->A3, 1, 2) = prm->old_shiftY[2] + prm->deltaY[2];
609  MAT_ELEM(prm->A3, 0, 0) = 1;
610  MAT_ELEM(prm->A3, 0, 1) = 0;
611  MAT_ELEM(prm->A3, 1, 0) = 0;
612  MAT_ELEM(prm->A3, 1, 1) = 1;
613  break;
614 
615 
616  default:
617  prm->deltaX[0] = x[idx + 1];
618  prm->deltaY[0] = x[idx + 2];
619  prm->deltaRot[0] = x[idx + 3];
620  prm->deltaTilt[0] = x[idx + 4];
621  prm->deltaPsi[0] = x[idx + 5];
622  prm->deltaDefocusU[0]=x[idx + 6];
623  prm->deltaDefocusV[0]=x[idx + 7];
624  prm->deltaDefocusAngle[0]=x[idx + 8];
625 
626  MAT_ELEM(prm->A1, 0, 2) = prm->old_shiftX[0] + prm->deltaX[0];
627  MAT_ELEM(prm->A1, 1, 2) = prm->old_shiftY[0] + prm->deltaY[0];
628  MAT_ELEM(prm->A1, 0, 0) = 1;
629  MAT_ELEM(prm->A1, 0, 1) = 0;
630  MAT_ELEM(prm->A1, 1, 0) = 0;
631  MAT_ELEM(prm->A1, 1, 1) = 1;
632  break;
633  }
634 
635  return prm->transformImageSph(x);
636 
637 }
638 
639 // Predict =================================================================
640 //#define DEBUG
641 void ProgForwardZernikeImages::processImage(const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut)
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 }
871 #undef DEBUG
872 
873 void ProgForwardZernikeImages::writeImageParameters(MDRow &row) {
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 }
933 
934 void ProgForwardZernikeImages::checkPoint() {
935  MDRowVec rowAppend;
936  MetaDataVec checkPoint;
937  getOutputMd().getRow(rowAppend, getOutputMd().lastRowId());
938  checkPoint.addRow(rowAppend);
939  checkPoint.append(Rerunable::getFileName());
940 }
941 
942 void ProgForwardZernikeImages::numCoefficients(int l1, int l2, int &vecSize)
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 }
957 
958 void ProgForwardZernikeImages::minimizepos(int L1, int l2, Matrix1D<double> &steps)
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 }
990 
991 void ProgForwardZernikeImages::fillVectorTerms(int l1, int l2, Matrix1D<int> &vL1, Matrix1D<int> &vN,
992  Matrix1D<int> &vL2, Matrix1D<int> &vM)
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 }
1013 
1014 void ProgForwardZernikeImages::updateCTFImage(double defocusU, double defocusV, double angle)
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 }
1022 
1023 template<ProgForwardZernikeImages::Direction DIRECTION>
1024 void ProgForwardZernikeImages::rotateCoefficients() {
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 }
1046 
1047 void ProgForwardZernikeImages::deformVol(MultidimArray<double> &mP, const MultidimArray<double> &mV, double &def,
1048  double rot, double tilt, double psi)
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 }
1142 
1143 Matrix1D<double> ProgForwardZernikeImages::weightsInterpolation3D(double x, double y, double z) {
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 }
1173 
1174 void ProgForwardZernikeImages::splattingAtPos(std::array<double, 3> r, double weight, MultidimArray<double> &mP, const MultidimArray<double> &mV) {
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 }
1193 
1194 double ProgForwardZernikeImages::bspline1(double x)
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 }
1204 
1205 double ProgForwardZernikeImages::bspline3(double x)
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 }
1223 
1224 void ProgForwardZernikeImages::removePixels()
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 }
1256 
1257 void ProgForwardZernikeImages::rotatePositions(double rot, double tilt, double psi)
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 }
1298 
1299 
1300 void ProgForwardZernikeImages::preComputeDF()
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 }
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::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
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
blobtype::alpha
double alpha
Smoothness parameter.
Definition: blobs.h:121
MDL_CTF_DEFOCUSU
Defocus U (Angstroms)
Definition: metadata_label.h:135
XMIPP_MAX
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
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
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
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
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
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
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
ProgForwardZernikeImages::rotatePositions
void rotatePositions(double rot, double tilt, double psi)
Definition: forward_zernike_images.cpp:1257
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
ProgForwardZernikeImages::fnOutDir
FileName fnOutDir
Output directory.
Definition: forward_zernike_images.h:51
ProgForwardZernikeImages::old_defocusV
std::vector< double > old_defocusV
Definition: forward_zernike_images.h:122
Mask
Definition: mask.h:360
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ProgForwardZernikeImages::deltaPsi
std::vector< double > deltaPsi
Definition: forward_zernike_images.h:114
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
y
static double * y
Definition: numerical_recipes.cpp:8487
ProgForwardZernikeImages::I
std::vector< Image< double > > I
Definition: forward_zernike_images.h:104
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
ProgForwardZernikeImages::deltaTilt
std::vector< double > deltaTilt
Definition: forward_zernike_images.h:114
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::~ProgForwardZernikeImages
~ProgForwardZernikeImages()
Destructor.
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::weightsInterpolation3D
Matrix1D< double > weightsInterpolation3D(double x, double y, double z)
Definition: forward_zernike_images.cpp:1143
MDL_ANGLE_TILT3
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:61
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
ProgForwardZernikeImages::Ifilteredp
std::vector< Image< double > > Ifilteredp
Definition: forward_zernike_images.h:106
MDL_ANGLE_PSI
Special label to be used when gathering MDs in MpiMetadataPrograms.
Definition: metadata_label.h:51
w
doublereal * w
Definition: numerical_recipes.cpp:2477
ProgForwardZernikeImages::deltaDefocusAngle
std::vector< double > deltaDefocusAngle
Definition: forward_zernike_images.h:122
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
Matrix2D::mdata
T * mdata
Definition: matrix2d.h:395
MDL_CTF_MODEL
Name for the CTF Model (std::string)
Definition: metadata_label.h:129
ProgForwardZernikeImages::deltaDefocusU
std::vector< double > deltaDefocusU
Definition: forward_zernike_images.h:122
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
ProgForwardZernikeImages::ProgForwardZernikeImages
ProgForwardZernikeImages()
Empty constructor.
Definition: forward_zernike_images.cpp:35
xmipp_image_extension.h
Rerunable::setFileName
void setFileName(const FileName &fn)
Definition: rerunable_program.h:82
z0
#define z0
ProgForwardZernikeImages::updateCTFImage
void updateCTFImage(double defocusU, double defocusV, double angle)
Definition: forward_zernike_images.cpp:1014
FINISHINGZ
#define FINISHINGZ(v)
Definition: multidim_array_base.h:71
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
ProgForwardZernikeImages::deltaX
std::vector< double > deltaX
Definition: forward_zernike_images.h:116
MetaDataVec::getRow
std::unique_ptr< MDRow > getRow(size_t id) override
Definition: metadata_vec.cpp:310
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
ProgForwardZernikeImages::writeImageParameters
virtual void writeImageParameters(MDRow &row)
Definition: forward_zernike_images.cpp:873
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
MDL_ENABLED
Is this image enabled? (int [-1 or 1])
Definition: metadata_label.h:231
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
MetaDataVec::addRow
size_t addRow(const MDRow &row) override
Definition: metadata_vec.cpp:225
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
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
continuousZernikeCost
double continuousZernikeCost(double *x, void *_prm)
Definition: forward_zernike_images.cpp:520
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
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
FOR_ALL_ELEMENTS_IN_MATRIX1D
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
ProgForwardZernikeImages::deltaY
std::vector< double > deltaY
Definition: forward_zernike_images.h:116
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
Rerunable::getFileName
const FileName & getFileName() const
Definition: rerunable_program.h:80
MDRow::getValue
T & getValue(MDLabel label)
Definition: metadata_row_base.h:97
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
ProgForwardZernikeImages::currentAngle
std::vector< double > currentAngle
Definition: forward_zernike_images.h:124
y0
#define y0
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
x0
#define x0
XX
#define XX(v)
Definition: matrix1d.h:85
xmipp_image_generic.h
CEIL
#define CEIL(x)
Definition: xmipp_macros.h:225
MDL_SPH_DEFORMATION
Deformation in voxels.
Definition: metadata_label.h:449
CTF
#define CTF
Definition: fourier_filter.h:81
MDL_COST
Cost for the image (double)
Definition: metadata_label.h:90
ProgForwardZernikeImages::defineParams
void defineParams()
Define parameters.
Definition: forward_zernike_images.cpp:105
ProgForwardZernikeImages::old_rot
std::vector< double > old_rot
Definition: forward_zernike_images.h:114
Mask::R1
double R1
Definition: mask.h:413
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
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
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
ProgForwardZernikeImages::processImage
void processImage(const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut)
Definition: forward_zernike_images.cpp:641
CTFDescription1D::Tm
double Tm
Sampling rate (A/pixel)
Definition: ctf.h:240
Mask::type
int type
Definition: mask.h:402
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
z
double z
Definition: numerical_recipes.cpp:8490
XmippProgram::addExampleLine
void addExampleLine(const char *example, bool verbatim=true)
Definition: xmipp_program.cpp:286
MDL_IMAGE1
Image associated to this object (std::string)
Definition: metadata_label.h:255
ProgForwardZernikeImages::z_clnm_diff
std::vector< double > z_clnm_diff
Definition: forward_zernike_images.h:151
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
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
j
#define j
Definition: numerical_recipes.cpp:2493
steps
double steps
Definition: numerical_recipes.cpp:5166
YY
#define YY(v)
Definition: matrix1d.h:93
m
int m
Definition: image_find_center.cpp:57
MDRow::getValueOrDefault
const T & getValueOrDefault(MDLabel label, const T &def) const
Definition: metadata_row_base.h:129
ProgForwardZernikeImages::numCoefficients
void numCoefficients(int l1, int l2, int &vecSize)
Length of coefficients vector.
Definition: forward_zernike_images.cpp:942
CTFDescription1D::K
double K
Global gain. By default, 1.
Definition: ctf.h:238
MDRow::setValue
void setValue(MDLabel label, const T &d, bool addLabel=true)
Definition: metadata_row_base.h:151
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
XmippMetadataProgram::show
void show() const override
Definition: xmipp_metadata_program.cpp:232
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
Mask::generate_mask
void generate_mask(bool apply_geo=false)
Definition: mask.cpp:1577
ProgForwardZernikeImages::readParams
void readParams()
Read argument from command line.
Definition: forward_zernike_images.cpp:46
MDRow::containsLabel
virtual bool containsLabel(MDLabel label) const =0
typeCast
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
MetaData::append
void append(const FileName &outFile) const
Definition: metadata_base.cpp:647
BINARY_CIRCULAR_MASK
#define BINARY_CIRCULAR_MASK
Definition: mask.h:365
ZernikeSphericalHarmonics
double ZernikeSphericalHarmonics(int n, int m, double xr, double yr, double zr, double r)
Definition: numerical_tools.cpp:1054
CTFDescription::produceSideInfo
void produceSideInfo()
Produce Side information.
Definition: ctf.cpp:1392
alglib::psi
double psi(const double x)
Definition: specialfunctions.cpp:2492
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
blobtype::order
int order
Derivation order and Bessel function order.
Definition: blobs.h:118
ProgForwardZernikeImages::splattingAtPos
void splattingAtPos(std::array< double, 3 > r, double weight, MultidimArray< double > &mP, const MultidimArray< double > &mV)
Definition: forward_zernike_images.cpp:1174
XmippMetadataProgram::keep_input_columns
bool keep_input_columns
Keep input metadata columns.
Definition: xmipp_metadata_program.h:105
ProgForwardZernikeImages::idx_z_clnm
std::vector< size_t > idx_z_clnm
Definition: forward_zernike_images.h:150
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
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
forward_zernike_images.h
r2
float r2
Definition: image_find_center.cpp:54
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
REALGAUSSIAN
#define REALGAUSSIAN
Definition: fourier_filter.h:87
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
ProgForwardZernikeImages::checkPoint
virtual void checkPoint()
For very long programs, it may be needed to write checkpoints.
Definition: forward_zernike_images.cpp:934
MDL_SHIFT_Y
Shift for the image in the Y axis (double)
Definition: metadata_label.h:435
XmippProgram::addUsageLine
void addUsageLine(const char *line, bool verbatim=false)
Definition: xmipp_program.cpp:282
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
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
Matrix1D::adaptForNumericalRecipes
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331
FourierFilter::generateMask
void generateMask(MultidimArray< double > &v)
Definition: fourier_filter.cpp:605
MDL_ANGLE_PSI3
Psi angle of an image (double,degrees)
Definition: metadata_label.h:53
ProgForwardZernikeImages::deltaDefocusV
std::vector< double > deltaDefocusV
Definition: forward_zernike_images.h:122
ProgForwardZernikeImages::transformImageSph
double transformImageSph(double *pclnm)
Definition: forward_zernike_images.cpp:291
STARTINGZ
#define STARTINGZ(v)
Definition: multidim_array_base.h:67
blobtype::radius
double radius
Spatial radius in Universal System units.
Definition: blobs.h:115
n
int * n
Definition: numerical_recipes.cpp:2229
MDL_IMAGE
Name of an image (std::string)
Definition: metadata_label.h:246
ZZ
#define ZZ(v)
Definition: matrix1d.h:101
XmippMetadataProgram::produces_a_metadata
bool produces_a_metadata
Indicate that the unique final output file is a Metadata.
Definition: xmipp_metadata_program.h:88
LOWPASS
#define LOWPASS
Definition: fourier_filter.h:77
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305
ProgForwardZernikeImages::df
MultidimArray< double > df
Definition: forward_zernike_images.h:149
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
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710
ProgForwardZernikeImages::deltaRot
std::vector< double > deltaRot
Definition: forward_zernike_images.h:114
XmippProgram::defaultComments
std::map< String, CommentList > defaultComments
Definition: xmipp_program.h:83
INNER_MASK
constexpr int INNER_MASK
Definition: mask.h:47
Generated by   doxygen 1.8.13