Xmipp  v3.23.11-Nereus
forward_zernike_images_priors.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_priors.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 #include "fstream"
34 
35 // Empty constructor =======================================================
36 ProgForwardZernikeImagesPriors::ProgForwardZernikeImagesPriors() : Rerunable("")
37 {
38  resume = false;
39  produces_a_metadata = true;
40  each_image_produces_an_output = false;
41  showOptimization = false;
42 }
43 
44 ProgForwardZernikeImagesPriors::~ProgForwardZernikeImagesPriors() = default;
45 
46 // Read arguments ==========================================================
47 void ProgForwardZernikeImagesPriors::readParams()
48 {
49  XmippMetadataProgram::readParams();
50  fnVolR = getParam("--ref");
51  fnMaskR = getParam("--mask");
52  fnPriors = getParam("--priors");
53  fnOutDir = getParam("--odir");
54  maxShift = getDoubleParam("--max_shift");
55  maxAngularChange = getDoubleParam("--max_angular_change");
56  maxResol = getDoubleParam("--max_resolution");
57  Ts = getDoubleParam("--sampling");
58  Rmax = getIntParam("--Rmax");
59  RmaxDef = getIntParam("--RDef");
60  optimizeAlignment = checkParam("--optimizeAlignment");
61  optimizeDefocus = checkParam("--optimizeDefocus");
62  phaseFlipped = checkParam("--phaseFlipped");
63  useCTF = checkParam("--useCTF");
64  L1 = getIntParam("--l1");
65  L2 = getIntParam("--l2");
66  loop_step = getIntParam("--step");
67  image_mode = getIntParam("--image_mode");
68  resume = checkParam("--resume");
69  lambda = getDoubleParam("--regularization");
70  Rerunable::setFileName(fnOutDir + "/sphDone.xmd");
71  blob_r = getDoubleParam("--blobr");
72  keep_input_columns = true;
73 }
74 
75 // Show ====================================================================
76 void ProgForwardZernikeImagesPriors::show()
77 {
78  if (!verbose)
79  return;
80  XmippMetadataProgram::show();
81  std::cout
82  << "Output directory: " << fnOutDir << std::endl
83  << "Reference volume: " << fnVolR << std::endl
84  << "Reference mask: " << fnMaskR << std::endl
85  << "Max. Shift: " << maxShift << std::endl
86  << "Max. Angular Change: " << maxAngularChange << std::endl
87  << "Max. Resolution: " << maxResol << std::endl
88  << "Sampling: " << Ts << std::endl
89  << "Max. Radius: " << Rmax << std::endl
90  << "Max. Radius Deform. " << RmaxDef << std::endl
91  << "Zernike Degree: " << L1 << std::endl
92  << "SH Degree: " << L2 << std::endl
93  << "Step: " << loop_step << std::endl
94  << "Correct CTF: " << useCTF << std::endl
95  << "Optimize alignment: " << optimizeAlignment << std::endl
96  << "Optimize defocus: " << optimizeDefocus << std::endl
97  << "Regularization: " << lambda << std::endl
98  << "Phase flipped: " << phaseFlipped << std::endl
99  << "Blob radius: " << blob_r << std::endl
100  << "Image mode: " << image_mode << std::endl
101  ;
102 }
103 
104 // usage ===================================================================
105 void ProgForwardZernikeImagesPriors::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(" --priors <metadata_file> : List of priors (deformation coefficients)");
115  addParamsLine(" [--mask <m=\"\">] : Reference volume");
116  addParamsLine(" [--odir <outputDir=\".\">] : Output directory");
117  addParamsLine(" [--max_shift <s=-1>] : Maximum shift allowed in pixels");
118  addParamsLine(" [--max_angular_change <a=5>] : Maximum angular change allowed (in degrees)");
119  addParamsLine(" [--max_resolution <f=4>] : Maximum resolution (A)");
120  addParamsLine(" [--sampling <Ts=1>] : Sampling rate (A/pixel)");
121  addParamsLine(" [--Rmax <R=-1>] : Maximum radius (px). -1=Half of volume size");
122  addParamsLine(" [--RDef <r=-1>] : Maximum radius of the deformation (px). -1=Half of volume size");
123  addParamsLine(" [--l1 <l1=3>] : Degree Zernike Polynomials=1,2,3,...");
124  addParamsLine(" [--l2 <l2=2>] : Harmonical depth of the deformation=1,2,3,...");
125  addParamsLine(" [--regularization <l=0.005>] : Regularization weight");
126  addParamsLine(" [--step <step=1>] : Voxel index step");
127  addParamsLine(" [--useCTF] : Correct CTF");
128  addParamsLine(" [--optimizeAlignment] : Optimize alignment");
129  addParamsLine(" [--optimizeDefocus] : Optimize defocus");
130  addParamsLine(" [--phaseFlipped] : Input images have been phase flipped");
131  addParamsLine(" [--blobr <b=1.5>] : 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_angular_sph_alignment -i anglesFromContinuousAssignment.xmd --ref reference.vol -o assigned_anglesAndDeformations.xmd --optimizeAlignment --optimizeDeformation --depth 1");
136 }
137 
138 void ProgForwardZernikeImagesPriors::preProcess()
139 {
140  V.read(fnVolR);
141  V().setXmippOrigin();
142  Xdim=XSIZE(V());
143  Vdeformed().initZeros(V());
144  Vdeformed().setXmippOrigin();
145 
146  std::string line;
147  std::ifstream priorsFile(fnPriors.getString());
148  bool first_line = true;
149  while (std::getline(priorsFile, line))
150  {
151  if (first_line)
152  {
153  std::vector<double> basisParams = string2vector(line);
154  L1 = (int)basisParams[0]; L2 = (int)basisParams[1]; RmaxDef = (int)basisParams[2];
155  first_line = false;
156  for (int idx=3; idx < basisParams.size(); idx++)
157  prior_deformations.push_back(basisParams[idx]);
158  }
159  else
160  {
161  priors.push_back(string2vector(line));
162  }
163  }
164 
165  // Check execution mode (single, pair, or triplet)
166  if (image_mode > 0)
167  {
168  num_images = image_mode;
169  }
170  else{
171  num_images = 1; // Single image
172  if (getInputMd()->containsLabel(MDL_IMAGE1) && !getInputMd()->containsLabel(MDL_IMAGE2))
173  {
174  num_images = 2; // Pair
175  }
176  else if (getInputMd()->containsLabel(MDL_IMAGE1) && getInputMd()->containsLabel(MDL_IMAGE2))
177  {
178  num_images = 3; // Triplet
179  }
180  }
181 
182  // Preallocate vectors (Size depends on image number)
183  fnImage.resize(num_images, "");
184  I.resize(num_images); Ifiltered.resize(num_images); Ifilteredp.resize(num_images);
185  P.resize(num_images);
186  old_rot.resize(num_images, 0.); old_tilt.resize(num_images, 0.); old_psi.resize(num_images, 0.);
187  deltaRot.resize(num_images, 0.); deltaTilt.resize(num_images, 0.); deltaPsi.resize(num_images, 0.);
188  old_shiftX.resize(num_images, 0.); old_shiftY.resize(num_images, 0.);
189  deltaX.resize(num_images, 0.); deltaY.resize(num_images, 0.);
190  old_defocusU.resize(num_images, 0.); old_defocusV.resize(num_images, 0.); old_defocusAngle.resize(num_images, 0.);
191  deltaDefocusU.resize(num_images, 0.); deltaDefocusV.resize(num_images, 0.); deltaDefocusAngle.resize(num_images, 0.);
192  currentDefocusU.resize(num_images, 0.); currentDefocusV.resize(num_images, 0.); currentAngle.resize(num_images, 0.);
193 
194  switch (num_images)
195  {
196  case 2:
197  Ifilteredp[0]().initZeros(Xdim,Xdim); Ifilteredp[1]().initZeros(Xdim,Xdim);
198  Ifilteredp[0]().setXmippOrigin(); Ifilteredp[1]().setXmippOrigin();
199  P[0]().initZeros(Xdim,Xdim); P[1]().initZeros(Xdim,Xdim);
200  break;
201  case 3:
202  Ifilteredp[0]().initZeros(Xdim,Xdim); Ifilteredp[1]().initZeros(Xdim,Xdim); Ifilteredp[2]().initZeros(Xdim,Xdim);
203  Ifilteredp[0]().setXmippOrigin(); Ifilteredp[1]().setXmippOrigin(); Ifilteredp[2]().setXmippOrigin();
204  P[0]().initZeros(Xdim,Xdim); P[1]().initZeros(Xdim,Xdim); P[2]().initZeros(Xdim,Xdim);
205  break;
206 
207  default:
208  Ifilteredp[0]().initZeros(Xdim,Xdim);
209  Ifilteredp[0]().setXmippOrigin();
210  P[0]().initZeros(Xdim,Xdim);
211  break;
212  }
213 
214  if (RmaxDef<0)
215  RmaxDef = Xdim/2;
216 
217  // Read Reference mask if avalaible (otherwise sphere of radius RmaxDef is used)
218  Mask mask;
219  mask.type = BINARY_CIRCULAR_MASK;
220  mask.mode = INNER_MASK;
221  if (fnMaskR != "") {
222  Image<double> aux;
223  aux.read(fnMaskR);
224  typeCast(aux(), V_mask);
225  V_mask.setXmippOrigin();
226  double Rmax2 = RmaxDef*RmaxDef;
227  for (int k=STARTINGZ(V_mask); k<=FINISHINGZ(V_mask); k++) {
228  for (int i=STARTINGY(V_mask); i<=FINISHINGY(V_mask); i++) {
229  for (int j=STARTINGX(V_mask); j<=FINISHINGX(V_mask); j++) {
230  double r2 = k*k + i*i + j*j;
231  if (r2>=Rmax2)
232  A3D_ELEM(V_mask,k,i,j) = 0;
233  }
234  }
235  }
236  }
237  else {
238  mask.R1 = RmaxDef;
239  mask.generate_mask(V());
240  V_mask = mask.get_binary_mask();
241  V_mask.setXmippOrigin();
242  }
243 
244  // Construct mask
245  if (Rmax<0)
246  Rmax=Xdim/2;
247  mask.R1 = Rmax;
248  mask.generate_mask(Xdim,Xdim);
249  mask2D=mask.get_binary_mask();
250 
251  // Low pass filter
252  filter.FilterBand=LOWPASS;
253  filter.w1=Ts/maxResol;
254  filter.raised_w=0.02;
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 = 3.6; // Smoothness parameter
282 
283 }
284 
285 void ProgForwardZernikeImagesPriors::finishProcessing() {
286  XmippMetadataProgram::finishProcessing();
287  rename(Rerunable::getFileName().c_str(), (fnOutDir + fn_out).c_str());
288 }
289 
290 double ProgForwardZernikeImagesPriors::transformImageSph(double *pc_priors)
291 {
292  const MultidimArray<double> &mV=V();
293  FOR_ALL_ELEMENTS_IN_MATRIX1D(c_priors)
294  VEC_ELEM(c_priors,i)=pc_priors[i+1];
295 
296  linearCombinationClnm();
297 
298  double deformation=0.0;
299  totalDeformation=0.0;
300 
301  P[0]().initZeros(Xdim,Xdim);
302  P[0]().setXmippOrigin();
303  double currentRot=old_rot[0] + deltaRot[0];
304  double currentTilt=old_tilt[0] + deltaTilt[0];
305  double currentPsi=old_psi[0] + deltaPsi[0];
306  deformVol(P[0](), mV, deformation, currentRot, currentTilt, currentPsi);
307 
308  double cost=0.0;
309  const MultidimArray<int> &mMask2D=mask2D;
310  double corr2 = 0.0;
311  double corr3 = 0.0;
312 
313  switch (num_images)
314  {
315  case 2:
316  {
317  P[1]().initZeros(Xdim,Xdim);
318  P[1]().setXmippOrigin();
319  currentRot = old_rot[1] + deltaRot[1];
320  currentTilt = old_tilt[1] + deltaTilt[1];
321  currentPsi = old_psi[1] + deltaPsi[1];
322  deformVol(P[1](), mV, deformation, currentRot, currentTilt, currentPsi);
323 
324  if (old_flip)
325  {
326  MAT_ELEM(A1, 0, 0) *= -1;
327  MAT_ELEM(A1, 0, 1) *= -1;
328  MAT_ELEM(A1, 0, 2) *= -1;
329  MAT_ELEM(A2, 0, 0) *= -1;
330  MAT_ELEM(A2, 0, 1) *= -1;
331  MAT_ELEM(A2, 0, 2) *= -1;
332  }
333  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[0](), Ifiltered[0](), A1,
334  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
335  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[1](), Ifiltered[1](), A2,
336  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
337  // filter.applyMaskSpace(P[1]());
338  const MultidimArray<double> mP2 = P[1]();
339  MultidimArray<double> &mI2filteredp = Ifilteredp[1]();
340  corr2 = correlationIndex(mI2filteredp, mP2, &mMask2D);
341  }
342  break;
343 
344  case 3:
345  {
346  P[1]().initZeros(Xdim,Xdim);
347  P[1]().setXmippOrigin();
348  currentRot = old_rot[1] + deltaRot[1];
349  currentTilt = old_tilt[1] + deltaTilt[1];
350  currentPsi = old_psi[1] + deltaPsi[1];
351  deformVol(P[1](), mV, deformation, currentRot, currentTilt, currentPsi);
352 
353  P[2]().initZeros(Xdim,Xdim);
354  P[2]().setXmippOrigin();
355  currentRot = old_rot[2] + deltaRot[2];
356  currentTilt = old_tilt[2] + deltaTilt[2];
357  currentPsi = old_psi[2] + deltaPsi[2];
358  deformVol(P[2](), mV, deformation, currentRot, currentTilt, currentPsi);
359 
360  if (old_flip)
361  {
362  MAT_ELEM(A1, 0, 0) *= -1;
363  MAT_ELEM(A1, 0, 1) *= -1;
364  MAT_ELEM(A1, 0, 2) *= -1;
365  MAT_ELEM(A2, 0, 0) *= -1;
366  MAT_ELEM(A2, 0, 1) *= -1;
367  MAT_ELEM(A2, 0, 2) *= -1;
368  MAT_ELEM(A3, 0, 0) *= -1;
369  MAT_ELEM(A3, 0, 1) *= -1;
370  MAT_ELEM(A3, 0, 2) *= -1;
371  }
372  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[0](), Ifiltered[0](), A1,
373  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
374  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[1](), Ifiltered[1](), A2,
375  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
376  applyGeometry(xmipp_transformation::LINEAR, Ifilteredp[2](), Ifiltered[2](), A3,
377  xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, 0.);
378  // filter.applyMaskSpace(P[1]());
379  // filter.applyMaskSpace(P[2]());
380  const MultidimArray<double> mP2 = P[1]();
381  const MultidimArray<double> mP3 = P[2]();
382  MultidimArray<double> &mI2filteredp = Ifilteredp[1]();
383  MultidimArray<double> &mI3filteredp = Ifilteredp[2]();
384  corr2 = correlationIndex(mI2filteredp, mP2, &mMask2D);
385  corr3 = correlationIndex(mI3filteredp, mP3, &mMask2D);
386  }
387  break;
388 
389  default:
390  if (old_flip)
391  {
392  MAT_ELEM(A1, 0, 0) *= -1;
393  MAT_ELEM(A1, 0, 1) *= -1;
394  MAT_ELEM(A1, 0, 2) *= -1;
395  }
396  applyGeometry(xmipp_transformation::LINEAR,Ifilteredp[0](),Ifiltered[0](),A1,
397  xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP,0.);
398  break;
399  }
400 
401  if (hasCTF)
402  {
403  double defocusU=old_defocusU[0]+deltaDefocusU[0];
404  double defocusV=old_defocusV[0]+deltaDefocusV[0];
405  double angle=old_defocusAngle[0]+deltaDefocusAngle[0];
406  if (defocusU!=currentDefocusU[0] || defocusV!=currentDefocusV[0] || angle!=currentAngle[0]) {
407  updateCTFImage(defocusU,defocusV,angle);
408  }
409  FilterCTF1.generateMask(P[0]());
410  if (phaseFlipped)
411  FilterCTF1.correctPhase();
412  FilterCTF1.applyMaskSpace(P[0]());
413  }
414  filter.applyMaskSpace(P[0]());
415 
416 
417  const MultidimArray<double> mP1=P[0]();
418  MultidimArray<double> &mI1filteredp=Ifilteredp[0]();
419  double corr1=correlationIndex(mI1filteredp,mP1,&mMask2D);
420 
421  switch (num_images)
422  {
423  case 2:
424  cost=-(corr1+corr2+corr3) / 2.0;
425  break;
426  case 3:
427  cost=-(corr1+corr2+corr3) / 3.0;
428  break;
429  default:
430  cost=-(corr1+corr2+corr3);
431  break;
432  }
433 
434 #ifdef DEBUG
435  std::cout << "A=" << A << std::endl;
436  Image<double> save;
437  save()=P();
438  save.write("PPPtheo.xmp");
439  save()=Ifilteredp();
440  save.write("PPPfilteredp.xmp");
441  save()=Ifiltered();
442  save.write("PPPfiltered.xmp");
443  // Vdeformed.write("PPPVdeformed.vol");
444  std::cout << "Cost=" << cost << " deformation=" << deformation << std::endl;
445  std::cout << "Press any key" << std::endl;
446  char c; std::cin >> c;
447 #endif
448 
449  if (showOptimization)
450  {
451  std::cout << "A1=" << A1 << std::endl;
452  Image<double> save;
453  save()=P[0]();
454  save.write("PPPtheo1.xmp");
455  save()=Ifilteredp[0]();
456  save.write("PPPfilteredp1.xmp");
457  save()=Ifiltered[0]();
458  save.write("PPPfiltered1.xmp");
459 
460  switch (num_images)
461  {
462  case 2:
463  {
464  std::cout << "A2=" << A2 << std::endl;
465  save()=P[1]();
466  save.write("PPPtheo2.xmp");
467  save()=Ifilteredp[1]();
468  save.write("PPPfilteredp2.xmp");
469  save()=Ifiltered[1]();
470  save.write("PPPfiltered2.xmp");
471  }
472  break;
473 
474  case 3:
475  {
476  std::cout << "A2=" << A2 << std::endl;
477  save()=P[1]();
478  save.write("PPPtheo2.xmp");
479  save()=Ifilteredp[1]();
480  save.write("PPPfilteredp2.xmp");
481  save()=Ifiltered[1]();
482  save.write("PPPfiltered2.xmp");
483 
484  std::cout << "A3=" << A3 << std::endl;
485  save()=P[2]();
486  save.write("PPPtheo3.xmp");
487  save()=Ifilteredp[2]();
488  save.write("PPPfilteredp3.xmp");
489  save()=Ifiltered[2]();
490  save.write("PPPfiltered3.xmp");
491  }
492  break;
493 
494  default:
495  break;
496  }
497  Vdeformed.write("PPPVdeformed.vol");
498  std::cout << "Deformation=" << totalDeformation << std::endl;
499  std::cout << "Press any key" << std::endl;
500  char c; std::cin >> c;
501  }
502 
503  double bound = 1.0;
504  double sum = 0.0;
505  for (int i=0; i < priors.size(); i++) {
506  sum += c_priors[i];
507  if (c_priors[i] < 0.0)
508  {
509  bound = 0.0;
510  break;
511  }
512  }
513 
514  if (sum > 1.0 && bound == 1.0)
515  bound = 0.0;
516 
517  if (showOptimization)
518  std::cout << cost << " " << deformation << std::endl;
519  return bound * cost + lambda * abs(deformation - prior_deformation);
520 }
521 
522 double continuousZernikePriorsCost(double *x, void *_prm)
523 {
524  ProgForwardZernikeImagesPriors *prm=(ProgForwardZernikeImagesPriors *)_prm;
525  int idx = prm->priors.size();
526  // TODO: Optimize parameters for each image (not sharing)
527  // prm->deltaDefocusU[0]=x[idx+6]; prm->deltaDefocusU[1]=x[idx+6];
528  // prm->deltaDefocusV[0]=x[idx+7]; prm->deltaDefocusV[1]=x[idx+7];
529  // prm->deltaDefocusAngle[0]=x[idx+8]; prm->deltaDefocusAngle[1]=x[idx+8];
530 
531  switch (prm->num_images)
532  {
533  case 2:
534  prm->deltaX[0] = x[idx + 1];
535  prm->deltaY[0] = x[idx + 3];
536  prm->deltaRot[0] = x[idx + 5];
537  prm->deltaTilt[0] = x[idx + 7];
538  prm->deltaPsi[0] = x[idx + 9];
539  // prm->deltaDefocusU[0]=x[idx + 11];
540  // prm->deltaDefocusV[0]=x[idx + 13];
541  // prm->deltaDefocusAngle[0]=x[idx + 15];
542 
543  prm->deltaX[1] = x[idx + 2];
544  prm->deltaY[1] = x[idx + 4];
545  prm->deltaRot[1] = x[idx + 6];
546  prm->deltaTilt[1] = x[idx + 8];
547  prm->deltaPsi[1] = x[idx + 10];
548  // prm->deltaDefocusU[1]=x[idx + 12];
549  // prm->deltaDefocusV[1]=x[idx + 14];
550  // prm->deltaDefocusAngle[1]=x[idx + 16];
551 
552  MAT_ELEM(prm->A1, 0, 2) = prm->old_shiftX[0] + prm->deltaX[0];
553  MAT_ELEM(prm->A1, 1, 2) = prm->old_shiftY[0] + prm->deltaY[0];
554  MAT_ELEM(prm->A1, 0, 0) = 1;
555  MAT_ELEM(prm->A1, 0, 1) = 0;
556  MAT_ELEM(prm->A1, 1, 0) = 0;
557  MAT_ELEM(prm->A1, 1, 1) = 1;
558 
559  MAT_ELEM(prm->A2, 0, 2) = prm->old_shiftX[1] + prm->deltaX[1];
560  MAT_ELEM(prm->A2, 1, 2) = prm->old_shiftY[1] + prm->deltaY[1];
561  MAT_ELEM(prm->A2, 0, 0) = 1;
562  MAT_ELEM(prm->A2, 0, 1) = 0;
563  MAT_ELEM(prm->A2, 1, 0) = 0;
564  MAT_ELEM(prm->A2, 1, 1) = 1;
565  break;
566 
567  case 3:
568  prm->deltaX[0] = x[idx + 1];
569  prm->deltaY[0] = x[idx + 4];
570  prm->deltaRot[0] = x[idx + 7];
571  prm->deltaTilt[0] = x[idx + 10];
572  prm->deltaPsi[0] = x[idx + 13];
573  // prm->deltaDefocusU[0]=x[idx + 16];
574  // prm->deltaDefocusV[0]=x[idx + 19];
575  // prm->deltaDefocusAngle[0]=x[idx + 22];
576 
577  prm->deltaX[1] = x[idx + 2];
578  prm->deltaY[1] = x[idx + 5];
579  prm->deltaRot[1] = x[idx + 8];
580  prm->deltaTilt[1] = x[idx + 11];
581  prm->deltaPsi[1] = x[idx + 14];
582  // prm->deltaDefocusU[1]=x[idx + 17];
583  // prm->deltaDefocusV[1]=x[idx + 20];
584  // prm->deltaDefocusAngle[1]=x[idx + 23];
585 
586  prm->deltaX[2] = x[idx + 3];
587  prm->deltaY[2] = x[idx + 6];
588  prm->deltaRot[2] = x[idx + 9];
589  prm->deltaTilt[2] = x[idx + 12];
590  prm->deltaPsi[2] = x[idx + 15];
591  // prm->deltaDefocusU[2]=x[idx + 18];
592  // prm->deltaDefocusV[2]=x[idx + 21];
593  // prm->deltaDefocusAngle[2]=x[idx + 24];
594 
595  MAT_ELEM(prm->A1, 0, 2) = prm->old_shiftX[0] + prm->deltaX[0];
596  MAT_ELEM(prm->A1, 1, 2) = prm->old_shiftY[0] + prm->deltaY[0];
597  MAT_ELEM(prm->A1, 0, 0) = 1;
598  MAT_ELEM(prm->A1, 0, 1) = 0;
599  MAT_ELEM(prm->A1, 1, 0) = 0;
600  MAT_ELEM(prm->A1, 1, 1) = 1;
601 
602  MAT_ELEM(prm->A2, 0, 2) = prm->old_shiftX[1] + prm->deltaX[1];
603  MAT_ELEM(prm->A2, 1, 2) = prm->old_shiftY[1] + prm->deltaY[1];
604  MAT_ELEM(prm->A2, 0, 0) = 1;
605  MAT_ELEM(prm->A2, 0, 1) = 0;
606  MAT_ELEM(prm->A2, 1, 0) = 0;
607  MAT_ELEM(prm->A2, 1, 1) = 1;
608 
609  MAT_ELEM(prm->A3, 0, 2) = prm->old_shiftX[2] + prm->deltaX[2];
610  MAT_ELEM(prm->A3, 1, 2) = prm->old_shiftY[2] + prm->deltaY[2];
611  MAT_ELEM(prm->A3, 0, 0) = 1;
612  MAT_ELEM(prm->A3, 0, 1) = 0;
613  MAT_ELEM(prm->A3, 1, 0) = 0;
614  MAT_ELEM(prm->A3, 1, 1) = 1;
615  break;
616 
617 
618  default:
619  prm->deltaX[0] = x[idx + 1];
620  prm->deltaY[0] = x[idx + 2];
621  prm->deltaRot[0] = x[idx + 3];
622  prm->deltaTilt[0] = x[idx + 4];
623  prm->deltaPsi[0] = x[idx + 5];
624  prm->deltaDefocusU[0]=x[idx + 6];
625  prm->deltaDefocusV[0]=x[idx + 7];
626  prm->deltaDefocusAngle[0]=x[idx + 8];
627 
628  MAT_ELEM(prm->A1, 0, 2) = prm->old_shiftX[0] + prm->deltaX[0];
629  MAT_ELEM(prm->A1, 1, 2) = prm->old_shiftY[0] + prm->deltaY[0];
630  MAT_ELEM(prm->A1, 0, 0) = 1;
631  MAT_ELEM(prm->A1, 0, 1) = 0;
632  MAT_ELEM(prm->A1, 1, 0) = 0;
633  MAT_ELEM(prm->A1, 1, 1) = 1;
634  break;
635  }
636 
637  return prm->transformImageSph(x);
638 
639 }
640 
641 // Predict =================================================================
642 void ProgForwardZernikeImagesPriors::processImage(const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut)
643 {
644  Matrix1D<double> steps;
645  // totalSize = 3*num_Z_coeff + num_images * 2 shifts + num_images * 3 angles + num_images * 3 CTF
646  algn_params = num_images * 5;
647  ctf_params = num_images * 3;
648  int totalSize = priors.size() + algn_params + ctf_params;
649  p.initZeros(totalSize);
650  c_priors.initZeros(totalSize);
651  clnm.initZeros(3*vecSize);
652 
653  flagEnabled=1;
654 
655  rowIn.getValueOrDefault(MDL_IMAGE, fnImage[0], "");
656  rowIn.getValueOrDefault(MDL_ANGLE_ROT, old_rot[0], 0.0);
657  rowIn.getValueOrDefault(MDL_ANGLE_TILT, old_tilt[0], 0.0);
658  rowIn.getValueOrDefault(MDL_ANGLE_PSI, old_psi[0], 0.0);
659  rowIn.getValueOrDefault(MDL_SHIFT_X, old_shiftX[0], 0.0);
660  rowIn.getValueOrDefault(MDL_SHIFT_Y, old_shiftY[0], 0.0);
661  I[0].read(fnImage[0]);
662  I[0]().setXmippOrigin();
663  Ifiltered[0]() = I[0]();
664  filter.applyMaskSpace(Ifiltered[0]());
665 
666  switch (num_images)
667  {
668  case 2:
669  rowIn.getValueOrDefault(MDL_IMAGE1, fnImage[1], "");
670  rowIn.getValueOrDefault(MDL_ANGLE_ROT2, old_rot[1], 0.0);
671  rowIn.getValueOrDefault(MDL_ANGLE_TILT2, old_tilt[1], 0.0);
672  rowIn.getValueOrDefault(MDL_ANGLE_PSI2, old_psi[1], 0.0);
673  rowIn.getValueOrDefault(MDL_SHIFT_X2, old_shiftX[1], 0.0);
674  rowIn.getValueOrDefault(MDL_SHIFT_Y2, old_shiftY[1], 0.0);
675  I[1].read(fnImage[1]);
676  I[1]().setXmippOrigin();
677  Ifiltered[1]() = I[1]();
678  filter.applyMaskSpace(Ifiltered[1]());
679 
680  if (verbose >= 2)
681  std::cout << "Processing Pair (" << fnImage[0] << "," << fnImage[1] << ")" << std::endl;
682  break;
683 
684  case 3:
685  rowIn.getValueOrDefault(MDL_IMAGE1, fnImage[1], "");
686  rowIn.getValueOrDefault(MDL_ANGLE_ROT2, old_rot[1], 0.0);
687  rowIn.getValueOrDefault(MDL_ANGLE_TILT2, old_tilt[1], 0.0);
688  rowIn.getValueOrDefault(MDL_ANGLE_PSI2, old_psi[1], 0.0);
689  rowIn.getValueOrDefault(MDL_SHIFT_X2, old_shiftX[1], 0.0);
690  rowIn.getValueOrDefault(MDL_SHIFT_Y2, old_shiftY[1], 0.0);
691  I[1].read(fnImage[1]);
692  I[1]().setXmippOrigin();
693  Ifiltered[1]() = I[1]();
694  filter.applyMaskSpace(Ifiltered[1]());
695 
696  rowIn.getValueOrDefault(MDL_IMAGE2, fnImage[2], "");
697  rowIn.getValueOrDefault(MDL_ANGLE_ROT3, old_rot[2], 0.0);
698  rowIn.getValueOrDefault(MDL_ANGLE_TILT3, old_tilt[2], 0.0);
699  rowIn.getValueOrDefault(MDL_ANGLE_PSI3, old_psi[2], 0.0);
700  rowIn.getValueOrDefault(MDL_SHIFT_X3, old_shiftX[2], 0.0);
701  rowIn.getValueOrDefault(MDL_SHIFT_Y3, old_shiftY[2], 0.0);
702  I[2].read(fnImage[2]);
703  I[2]().setXmippOrigin();
704  Ifiltered[2]() = I[2]();
705  filter.applyMaskSpace(Ifiltered[2]());
706 
707  if (verbose >= 2)
708  std::cout << "Processing Triplet (" << fnImage[0] << "," << fnImage[1] << "," << fnImage[2] << ")" << std::endl;
709  break;
710 
711  default:
712  if (verbose >= 2)
713  std::cout << "Processing Image (" << fnImage[0] << ")" << std::endl;
714  break;
715  }
716 
717  if (rowIn.containsLabel(MDL_FLIP))
718  rowIn.getValue(MDL_FLIP,old_flip);
719  else
720  old_flip = false;
721 
722  // FIXME: Add defocus per image and make CTF correction available
723  if ((rowIn.containsLabel(MDL_CTF_DEFOCUSU) || rowIn.containsLabel(MDL_CTF_MODEL)) && useCTF)
724  {
725  hasCTF=true;
726  FilterCTF1.ctf.readFromMdRow(rowIn);
727  FilterCTF1.ctf.Tm = Ts;
728  FilterCTF1.ctf.produceSideInfo();
729  old_defocusU[0]=FilterCTF1.ctf.DeltafU;
730  old_defocusV[0]=FilterCTF1.ctf.DeltafV;
731  old_defocusAngle[0]=FilterCTF1.ctf.azimuthal_angle;
732  }
733  else
734  hasCTF=false;
735 
736  prior_deformation = 0.0;
737  optimalPowellOrder();
738  p.initZeros(totalSize);
739  c_priors.initZeros(totalSize);
740  clnm.initZeros(3*vecSize);
741  if (!do_not_init)
742  {
743  p[0] = 1.0;
744  c_priors[0] = 1.0;
745  clnm = priors[0];
746  prior_deformation = prior_deformations[0];
747  }
748 
749  if (verbose >= 2)
750  {
751  std::cout << std::endl;
752  std::cout << "------------------------------ Fiding priors linear combination ----------------------------" << std::endl;
753  }
754  steps.clear();
755  steps.initZeros(totalSize);
756 
757  // Optimize
758  double cost = -1;
759  try
760  {
761  cost = 1e38;
762  int iter;
763  if (optimizeAlignment)
764  {
765  int init = steps.size() - algn_params - ctf_params;
766  int end = steps.size() - ctf_params;
767  for (int i = init; i < end; i++)
768  steps(i) = 1.;
769  }
770  if (optimizeDefocus)
771  {
772  int init = steps.size() - ctf_params;
773  int end = steps.size();
774  for (int i = init; i < end; i++)
775  steps(i) = 1.;
776  }
777  minimizepos(steps);
778  steps_cp = steps;
779  powellOptimizer(p, 1, totalSize, &continuousZernikePriorsCost, this, 0.01, cost, iter, steps, verbose >= 2);
780 
781  if (verbose >= 3)
782  {
783  showOptimization = true;
784  continuousZernikePriorsCost(p.adaptForNumericalRecipes(), this);
785  showOptimization = false;
786  }
787 
788  if (cost > 0)
789  {
790  flagEnabled = -1;
791  p.initZeros();
792  }
793  cost = -cost;
794  correlation = cost;
795  if (verbose >= 2)
796  {
797  std::cout << std::endl;
798  for (int j = 1; j < 4; j++)
799  {
800  switch (j)
801  {
802  case 1:
803  std::cout << "X Coefficients=(";
804  break;
805  case 2:
806  std::cout << "Y Coefficients=(";
807  break;
808  case 3:
809  std::cout << "Z Coefficients=(";
810  break;
811  default:
812  // Default case will be never reached
813  break;
814  }
815  for (int i = (j - 1) * vecSize; i < j * vecSize; i++)
816  {
817  std::cout << clnm(i);
818  if (i < j * vecSize - 1)
819  std::cout << ",";
820  }
821  std::cout << ")" << std::endl;
822  }
823  std::cout << "Radius=" << RmaxDef << std::endl;
824  std::cout << " Dshift=(" << p(totalSize - 5) << "," << p(totalSize - 4) << ") "
825  << "Drot=" << p(totalSize - 3) << " Dtilt=" << p(totalSize - 2)
826  << " Dpsi=" << p(totalSize - 1) << std::endl;
827  std::cout << " Total deformation=" << totalDeformation << std::endl;
828  std::cout << std::endl;
829  }
830  }
831  catch (XmippError &XE)
832  {
833  std::cerr << XE.what() << std::endl;
834  std::cerr << "Warning: Cannot refine " << fnImg << std::endl;
835  flagEnabled = -1;
836  }
837 
838  //AJ NEW
839  writeImageParameters(rowOut);
840  //END AJ
841 
842 }
843 
844 void ProgForwardZernikeImagesPriors::writeImageParameters(MDRow &row) {
845  int pos = priors.size();
846  if (flagEnabled==1) {
847  row.setValue(MDL_ENABLED, 1);
848  }
849  else {
850  row.setValue(MDL_ENABLED, -1);
851  }
852 
853  switch (num_images)
854  {
855  case 2:
856  row.setValue(MDL_ANGLE_ROT, old_rot[0]+p(pos+4));
857  row.setValue(MDL_ANGLE_ROT2, old_rot[1]+p(pos+5));
858  row.setValue(MDL_ANGLE_TILT, old_tilt[0]+p(pos+6));
859  row.setValue(MDL_ANGLE_TILT2, old_tilt[1]+p(pos+7));
860  row.setValue(MDL_ANGLE_PSI, old_psi[0]+p(pos+8));
861  row.setValue(MDL_ANGLE_PSI2, old_psi[1]+p(pos+9));
862  row.setValue(MDL_SHIFT_X, old_shiftX[0]+p(pos));
863  row.setValue(MDL_SHIFT_X2, old_shiftX[1]+p(pos+1));
864  row.setValue(MDL_SHIFT_Y, old_shiftY[0]+p(pos+2));
865  row.setValue(MDL_SHIFT_Y2, old_shiftY[1]+p(pos+3));
866  break;
867 
868  case 3:
869  row.setValue(MDL_ANGLE_ROT, old_rot[0]+p(pos+6));
870  row.setValue(MDL_ANGLE_ROT2, old_rot[1]+p(pos+7));
871  row.setValue(MDL_ANGLE_ROT3, old_rot[2]+p(pos+8));
872  row.setValue(MDL_ANGLE_TILT, old_tilt[0]+p(pos+9));
873  row.setValue(MDL_ANGLE_TILT2, old_tilt[1]+p(pos+10));
874  row.setValue(MDL_ANGLE_TILT3, old_tilt[2]+p(pos+11));
875  row.setValue(MDL_ANGLE_PSI, old_psi[0]+p(pos+12));
876  row.setValue(MDL_ANGLE_PSI2, old_psi[1]+p(pos+13));
877  row.setValue(MDL_ANGLE_PSI3, old_psi[2]+p(pos+14));
878  row.setValue(MDL_SHIFT_X, old_shiftX[0]+p(pos));
879  row.setValue(MDL_SHIFT_X2, old_shiftX[1]+p(pos+1));
880  row.setValue(MDL_SHIFT_X3, old_shiftX[2]+p(pos+2));
881  row.setValue(MDL_SHIFT_Y, old_shiftY[0]+p(pos+3));
882  row.setValue(MDL_SHIFT_Y2, old_shiftY[1]+p(pos+4));
883  row.setValue(MDL_SHIFT_Y3, old_shiftY[2]+p(pos+5));
884  break;
885 
886  default:
887  row.setValue(MDL_ANGLE_ROT, old_rot[0]+p(pos+2));
888  row.setValue(MDL_ANGLE_TILT, old_tilt[0]+p(pos+3));
889  row.setValue(MDL_ANGLE_PSI, old_psi[0]+p(pos+4));
890  row.setValue(MDL_SHIFT_X, old_shiftX[0]+p(pos));
891  row.setValue(MDL_SHIFT_Y, old_shiftY[0]+p(pos+1));
892  break;
893  }
894 
895  row.setValue(MDL_SPH_DEFORMATION, totalDeformation);
896  std::vector<double> vectortemp;
897  size_t end_clnm = VEC_XSIZE(clnm);
898  for (int j = 0; j < end_clnm; j++) {
899  vectortemp.push_back(clnm(j));
900  }
901  row.setValue(MDL_SPH_COEFFICIENTS, vectortemp);
902  row.setValue(MDL_COST, correlation);
903 }
904 
905 void ProgForwardZernikeImagesPriors::checkPoint() {
906  MDRowVec rowAppend;
907  MetaDataVec checkPoint;
908  getOutputMd().getRow(rowAppend, getOutputMd().lastRowId());
909  checkPoint.addRow(rowAppend);
910  checkPoint.append(Rerunable::getFileName());
911 }
912 
913 void ProgForwardZernikeImagesPriors::numCoefficients(int l1, int l2, int &vecSize)
914 {
915  for (int h=0; h<=l2; h++)
916  {
917  int numSPH = 2*h+1;
918  int count=l1-h+1;
919  int numEven=(count>>1)+(count&1 && !(h&1));
920  if (h%2 == 0) {
921  vecSize += numSPH*numEven;
922  }
923  else {
924  vecSize += numSPH*(l1-h+1-numEven);
925  }
926  }
927 }
928 
929 void ProgForwardZernikeImagesPriors::minimizepos(Matrix1D<double> &steps)
930 {
931  int size = priors.size();
932  for (int idx=0; idx<size; idx++) {
933  VEC_ELEM(steps,idx) = 1.;
934  }
935 }
936 
937 void ProgForwardZernikeImagesPriors::fillVectorTerms(int l1, int l2, Matrix1D<int> &vL1, Matrix1D<int> &vN,
938  Matrix1D<int> &vL2, Matrix1D<int> &vM)
939 {
940  int idx = 0;
941  vL1.initZeros(vecSize);
942  vN.initZeros(vecSize);
943  vL2.initZeros(vecSize);
944  vM.initZeros(vecSize);
945  for (int h=0; h<=l2; h++) {
946  int totalSPH = 2*h+1;
947  int aux = std::floor(totalSPH/2);
948  for (int l=h; l<=l1; l+=2) {
949  for (int m=0; m<totalSPH; m++) {
950  VEC_ELEM(vL1,idx) = l;
951  VEC_ELEM(vN,idx) = h;
952  VEC_ELEM(vL2,idx) = h;
953  VEC_ELEM(vM,idx) = m-aux;
954  idx++;
955  }
956  }
957  }
958 }
959 
960 void ProgForwardZernikeImagesPriors::updateCTFImage(double defocusU, double defocusV, double angle)
961 {
962  FilterCTF1.ctf.K=1; // get pure CTF with no envelope
963  currentDefocusU[0]=FilterCTF1.ctf.DeltafU=defocusU;
964  currentDefocusV[0]=FilterCTF1.ctf.DeltafV=defocusV;
965  currentAngle[0]=FilterCTF1.ctf.azimuthal_angle=angle;
966  FilterCTF1.ctf.produceSideInfo();
967 }
968 
969 void ProgForwardZernikeImagesPriors::deformVol(MultidimArray<double> &mP, const MultidimArray<double> &mV, double &def,
970  double rot, double tilt, double psi)
971 {
972  size_t idxY0=vecSize;
973  double Ncount=0.0;
974  double modg=0.0;
975  double diff2=0.0;
976 
977  def=0.0;
978  size_t idxZ0=2*idxY0;
979  double RmaxF=RmaxDef;
980  double RmaxF2=RmaxF*RmaxF;
981  double iRmaxF=1.0/RmaxF;
982  // Rotation Matrix
983  Matrix2D<double> R;
984  R.initIdentity(3);
985  Euler_angles2matrix(rot, tilt, psi, R, false);
986 
987  // TODO: Poner primero i y j en el loop, acumular suma y guardar al final
988  const auto lastZ = FINISHINGZ(mV);
989  const auto lastY = FINISHINGY(mV);
990  const auto lastX = FINISHINGX(mV);
991  for (int k=STARTINGZ(mV); k<=lastZ; k+=loop_step)
992  {
993  for (int i=STARTINGY(mV); i<=lastY; i+=loop_step)
994  {
995  for (int j=STARTINGX(mV); j<=lastX; j+=loop_step)
996  {
997  if (A3D_ELEM(V_mask,k,i,j) == 1) {
998  double gx=0.0, gy=0.0, gz=0.0;
999  double k2=k*k;
1000  double kr=k*iRmaxF;
1001  double k2i2=k2+i*i;
1002  double ir=i*iRmaxF;
1003  double r2=k2i2+j*j;
1004  double jr=j*iRmaxF;
1005  double rr=sqrt(r2)*iRmaxF;
1006  for (size_t idx = 0; idx < idxY0; idx++)
1007  {
1008  auto l1 = VEC_ELEM(vL1,idx);
1009  auto n = VEC_ELEM(vN,idx);
1010  auto l2 = VEC_ELEM(vL2,idx);
1011  auto m = VEC_ELEM(vM,idx);
1012  auto zsph=ZernikeSphericalHarmonics(l1,n,l2,m,jr,ir,kr,rr);
1013  auto c = std::array<double, 3>{};
1014  c[0] = VEC_ELEM(clnm,idx);
1015  c[1] = VEC_ELEM(clnm,idx+idxY0);
1016  c[2] = VEC_ELEM(clnm,idx+idxZ0);
1017  if (rr>0 || l2==0) {
1018  gx += c[0] *(zsph);
1019  gy += c[1] *(zsph);
1020  gz += c[2] *(zsph);
1021  }
1022  }
1023 
1024  auto pos = std::array<double, 3>{};
1025  double r_x = j + gx;
1026  double r_y = i + gy;
1027  double r_z = k + gz;
1028  pos[0] = R.mdata[0] * r_x + R.mdata[1] * r_y + R.mdata[2] * r_z;
1029  pos[1] = R.mdata[3] * r_x + R.mdata[4] * r_y + R.mdata[5] * r_z;
1030  pos[2] = R.mdata[6] * r_x + R.mdata[7] * r_y + R.mdata[8] * r_z;
1031 
1032  double voxel_mV = A3D_ELEM(mV,k,i,j);
1033  splattingAtPos(pos, voxel_mV, mP, mV);
1034  modg += gx*gx+gy*gy+gz*gz;
1035  Ncount++;
1036  }
1037  }
1038  }
1039  }
1040 
1041  def = sqrt(modg/Ncount);
1042  totalDeformation = def;
1043 }
1044 
1045 Matrix1D<double> ProgForwardZernikeImagesPriors::weightsInterpolation3D(double x, double y, double z) {
1046  Matrix1D<double> w;
1047  w.initZeros(8);
1048 
1049  int x0 = FLOOR(x);
1050  double fx0 = x - x0;
1051  int x1 = x0 + 1;
1052  double fx1 = x1 - x;
1053 
1054  int y0 = FLOOR(y);
1055  double fy0 = y - y0;
1056  int y1 = y0 + 1;
1057  double fy1 = y1 - y;
1058 
1059  int z0 = FLOOR(z);
1060  double fz0 = z - z0;
1061  int z1 = z0 + 1;
1062  double fz1 = z1 - z;
1063 
1064  VEC_ELEM(w,0) = fx1 * fy1 * fz1; // w000 (x0,y0,z0)
1065  VEC_ELEM(w,1) = fx1 * fy1 * fz0; // w001 (x0,y0,z1)
1066  VEC_ELEM(w,2) = fx1 * fy0 * fz1; // w010 (x0,y1,z0)
1067  VEC_ELEM(w,3) = fx1 * fy0 * fz0; // w011 (x0,y1,z1)
1068  VEC_ELEM(w,4) = fx0 * fy1 * fz1; // w100 (x1,y0,z0)
1069  VEC_ELEM(w,5) = fx0 * fy1 * fz0; // w101 (x1,y0,z1)
1070  VEC_ELEM(w,6) = fx0 * fy0 * fz1; // w110 (x1,y1,z0)
1071  VEC_ELEM(w,7) = fx0 * fy0 * fz0; // w111 (x1,y1,z1)
1072 
1073  return w;
1074 }
1075 
1076 void ProgForwardZernikeImagesPriors::splattingAtPos(std::array<double, 3> r, double weight, MultidimArray<double> &mP, const MultidimArray<double> &mV) {
1077  // Find the part of the volume that must be updated
1078  double x_pos = r[0];
1079  double y_pos = r[1];
1080  double z_pos = r[2];
1081  int i0 = XMIPP_MAX(FLOOR(y_pos - blob_r), STARTINGY(mV));
1082  int iF = XMIPP_MIN(CEIL(y_pos + blob_r), FINISHINGY(mV));
1083  int j0 = XMIPP_MAX(FLOOR(x_pos - blob_r), STARTINGX(mV));
1084  int jF = XMIPP_MIN(CEIL(x_pos + blob_r), FINISHINGX(mV));
1085  for (int i = i0; i <= iF; i++)
1086  {
1087  double y_val = bspline1(i - y_pos);
1088  for (int j = j0; j <= jF; j++)
1089  {
1090  double x_val = bspline1(j - x_pos);
1091  A2D_ELEM(mP, i, j) += weight * x_val * y_val;
1092  }
1093  }
1094 }
1095 
1096 std::vector<double> ProgForwardZernikeImagesPriors::string2vector(std::string const &s) const
1097 {
1098  std::stringstream iss(s);
1099  double number;
1100  std::vector<double> v;
1101  while (iss >> number)
1102  v.push_back(number);
1103  return v;
1104 }
1105 
1106 void ProgForwardZernikeImagesPriors::linearCombinationClnm()
1107 {
1108  clnm.initZeros(3*vecSize);
1109  for (int idx=0; idx < 3*vecSize; idx++)
1110  {
1111  for (int idc = 0; idc < priors.size(); idc++)
1112  clnm[idx] += c_priors[idc] * priors[idc][idx];
1113  }
1114 }
1115 
1116 void ProgForwardZernikeImagesPriors::optimalPowellOrder()
1117 {
1118  do_not_init = false;
1119  std::vector<double> z(priors[0].size(), 0.0);
1120  priors.push_back(z);
1121  Matrix1D<double> cost_vec;
1122  cost_vec.initZeros(priors.size());
1123  int totalSize = priors.size() + algn_params + ctf_params;
1124 
1125  // Initialize some values
1126  // TODO: Optimize parameters for each image (not sharing)
1127  // prm->deltaDefocusU[0]=x[idx+6]; prm->deltaDefocusU[1]=x[idx+6];
1128  // prm->deltaDefocusV[0]=x[idx+7]; prm->deltaDefocusV[1]=x[idx+7];
1129  // prm->deltaDefocusAngle[0]=x[idx+8]; prm->deltaDefocusAngle[1]=x[idx+8];
1130 
1131  switch (num_images)
1132  {
1133  case 2:
1134  deltaX[0] = 0;
1135  deltaY[0] = 0;
1136  deltaRot[0] = 0;
1137  deltaTilt[0] = 0;
1138  deltaPsi[0] = 0;
1139  // prm->deltaDefocusU[0]=x[idx + 11];
1140  // prm->deltaDefocusV[0]=x[idx + 13];
1141  // prm->deltaDefocusAngle[0]=x[idx + 15];
1142 
1143  deltaX[1] = 0;
1144  deltaY[1] = 0;
1145  deltaRot[1] = 0;
1146  deltaTilt[1] = 0;
1147  deltaPsi[1] = 0;
1148  // prm->deltaDefocusU[1]=x[idx + 12];
1149  // prm->deltaDefocusV[1]=x[idx + 14];
1150  // prm->deltaDefocusAngle[1]=x[idx + 16];
1151 
1152  MAT_ELEM(A1, 0, 2) = old_shiftX[0];
1153  MAT_ELEM(A1, 1, 2) = old_shiftY[0];
1154  MAT_ELEM(A1, 0, 0) = 1;
1155  MAT_ELEM(A1, 0, 1) = 0;
1156  MAT_ELEM(A1, 1, 0) = 0;
1157  MAT_ELEM(A1, 1, 1) = 1;
1158 
1159  MAT_ELEM(A2, 0, 2) = old_shiftX[1];
1160  MAT_ELEM(A2, 1, 2) = old_shiftY[1];
1161  MAT_ELEM(A2, 0, 0) = 1;
1162  MAT_ELEM(A2, 0, 1) = 0;
1163  MAT_ELEM(A2, 1, 0) = 0;
1164  MAT_ELEM(A2, 1, 1) = 1;
1165  break;
1166 
1167  case 3:
1168  deltaX[0] = 0;
1169  deltaY[0] = 0;
1170  deltaRot[0] = 0;
1171  deltaTilt[0] = 0;
1172  deltaPsi[0] = 0;
1173  // prm->deltaDefocusU[0]=x[idx + 16];
1174  // prm->deltaDefocusV[0]=x[idx + 19];
1175  // prm->deltaDefocusAngle[0]=x[idx + 22];
1176 
1177  deltaX[1] = 0;
1178  deltaY[1] = 0;
1179  deltaRot[1] = 0;
1180  deltaTilt[1] = 0;
1181  deltaPsi[1] = 0;
1182  // prm->deltaDefocusU[1]=x[idx + 17];
1183  // prm->deltaDefocusV[1]=x[idx + 20];
1184  // prm->deltaDefocusAngle[1]=x[idx + 23];
1185 
1186  deltaX[2] = 0;
1187  deltaY[2] = 0;
1188  deltaRot[2] = 0;
1189  deltaTilt[2] = 0;
1190  deltaPsi[2] = 0;
1191  // prm->deltaDefocusU[2]=x[idx + 18];
1192  // prm->deltaDefocusV[2]=x[idx + 21];
1193  // prm->deltaDefocusAngle[2]=x[idx + 24];
1194 
1195  MAT_ELEM(A1, 0, 2) = old_shiftX[0];
1196  MAT_ELEM(A1, 1, 2) = old_shiftY[0];
1197  MAT_ELEM(A1, 0, 0) = 1;
1198  MAT_ELEM(A1, 0, 1) = 0;
1199  MAT_ELEM(A1, 1, 0) = 0;
1200  MAT_ELEM(A1, 1, 1) = 1;
1201 
1202  MAT_ELEM(A2, 0, 2) = old_shiftX[1];
1203  MAT_ELEM(A2, 1, 2) = old_shiftY[1];
1204  MAT_ELEM(A2, 0, 0) = 1;
1205  MAT_ELEM(A2, 0, 1) = 0;
1206  MAT_ELEM(A2, 1, 0) = 0;
1207  MAT_ELEM(A2, 1, 1) = 1;
1208 
1209  MAT_ELEM(A3, 0, 2) = old_shiftX[2];
1210  MAT_ELEM(A3, 1, 2) = old_shiftY[2];
1211  MAT_ELEM(A3, 0, 0) = 1;
1212  MAT_ELEM(A3, 0, 1) = 0;
1213  MAT_ELEM(A3, 1, 0) = 0;
1214  MAT_ELEM(A3, 1, 1) = 1;
1215  break;
1216 
1217 
1218  default:
1219  deltaX[0] = 0;
1220  deltaY[0] = 0;
1221  deltaRot[0] = 0;
1222  deltaTilt[0] = 0;
1223  deltaPsi[0] = 0;
1224  deltaDefocusU[0] = 0;
1225  deltaDefocusV[0] = 0;
1226  deltaDefocusAngle[0] = 0;
1227 
1228  MAT_ELEM(A1, 0, 2) = old_shiftX[0];
1229  MAT_ELEM(A1, 1, 2) = old_shiftY[0];
1230  MAT_ELEM(A1, 0, 0) = 1;
1231  MAT_ELEM(A1, 0, 1) = 0;
1232  MAT_ELEM(A1, 1, 0) = 0;
1233  MAT_ELEM(A1, 1, 1) = 1;
1234  break;
1235  }
1236 
1237  auto lambda_aux = lambda;
1238  lambda = 0.0;
1239  for (int idx=0; idx < priors.size(); idx++)
1240  {
1241  clnm.initZeros(3*vecSize);
1242  p.initZeros(totalSize);
1243  p[idx] = 1.0;
1244  double *pptr = p.adaptForNumericalRecipes();
1245  cost_vec[idx] = transformImageSph(pptr);
1246  }
1247  lambda = lambda_aux;
1248 
1249  std::vector<std::size_t> index_vec;
1250  for (std::size_t i = 0; i != priors.size(); ++i)
1251  {
1252  index_vec.push_back(i);
1253  }
1254 
1255  std::sort(
1256  index_vec.begin(), index_vec.end(),
1257  [&](std::size_t a, std::size_t b)
1258  { return cost_vec[a] < cost_vec[b]; });
1259 
1260  auto priors_aux = priors;
1261  auto prior_deformations_aux = prior_deformations;
1262  priors.clear();
1263  prior_deformations.clear();
1264  for (std::size_t i = 0; i != index_vec.size(); ++i)
1265  {
1266  auto v = priors_aux[index_vec[i]];
1267  auto d = prior_deformations_aux[index_vec[i]];
1268  bool zeros = std::all_of(v.begin(), v.end(), [](double j) { return j==0.0; });
1269  if (!zeros)
1270  {
1271  priors.push_back(v);
1272  prior_deformations.push_back(d);
1273  }
1274  else if (zeros && i==0)
1275  do_not_init = true;
1276  }
1277 }
1278 
1279 double ProgForwardZernikeImagesPriors::bspline1(double x)
1280 {
1281  double m = 1 / blob_r;
1282  if (0. < x && x < blob_r)
1283  return m * (blob_r - x);
1284  else if (-blob_r < x && x <= 0.)
1285  return m * (blob_r + x);
1286  else
1287  return 0.;
1288 }
ProgForwardZernikeImagesPriors::splattingAtPos
void splattingAtPos(std::array< double, 3 > r, double weight, MultidimArray< double > &mP, const MultidimArray< double > &mV)
Definition: forward_zernike_images_priors.cpp:1076
ProgForwardZernikeImagesPriors::Ifiltered
std::vector< Image< double > > Ifiltered
Definition: forward_zernike_images_priors.h:107
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
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
ProgForwardZernikeImagesPriors::checkPoint
virtual void checkPoint()
For very long programs, it may be needed to write checkpoints.
Definition: forward_zernike_images_priors.cpp:905
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
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
FourierFilter::ctf
CTFDescription ctf
Definition: fourier_filter.h:134
Matrix1D::clear
void clear()
Definition: matrix1d.cpp:67
ProgForwardZernikeImagesPriors::numCoefficients
void numCoefficients(int l1, int l2, int &vecSize)
Length of coefficients vector.
Definition: forward_zernike_images_priors.cpp:913
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
forward_zernike_images_priors.h
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
Mask
Definition: mask.h:360
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
y
static double * y
Definition: numerical_recipes.cpp:8487
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
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
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
ProgForwardZernikeImagesPriors::priors
std::vector< std::vector< double > > priors
Definition: forward_zernike_images_priors.h:99
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
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
w
doublereal * w
Definition: numerical_recipes.cpp:2477
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
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
xmipp_image_extension.h
Rerunable::setFileName
void setFileName(const FileName &fn)
Definition: rerunable_program.h:82
z0
#define z0
FINISHINGZ
#define FINISHINGZ(v)
Definition: multidim_array_base.h:71
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
ProgForwardZernikeImagesPriors::string2vector
std::vector< double > string2vector(std::string const &s) const
Definition: forward_zernike_images_priors.cpp:1096
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
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
d
doublereal * d
Definition: numerical_recipes.cpp:2230
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
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
b
doublereal * b
Definition: numerical_recipes.cpp:2230
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
y0
#define y0
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
x0
#define x0
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
Mask::R1
double R1
Definition: mask.h:413
MDL_FLIP
Flip the image? (bool)
Definition: metadata_label.h:236
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
ProgForwardZernikeImagesPriors::~ProgForwardZernikeImagesPriors
~ProgForwardZernikeImagesPriors()
Destructor.
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
Mask::type
int type
Definition: mask.h:402
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
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
sort
void sort(struct DCEL_T *dcel)
Definition: sorting.cpp:18
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
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
continuousZernikePriorsCost
double continuousZernikePriorsCost(double *x, void *_prm)
Definition: forward_zernike_images_priors.cpp:522
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
ProgForwardZernikeImagesPriors::readParams
void readParams()
Read argument from command line.
Definition: forward_zernike_images_priors.cpp:47
Mask::generate_mask
void generate_mask(bool apply_geo=false)
Definition: mask.cpp:1577
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
FileName::getString
String getString() const
Definition: xmipp_filename.h:158
blobtype::order
int order
Derivation order and Bessel function order.
Definition: blobs.h:118
ProgForwardZernikeImagesPriors::processImage
void processImage(const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut)
Definition: forward_zernike_images_priors.cpp:642
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
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
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
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
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
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
ProgForwardZernikeImagesPriors::minimizepos
void minimizepos(Matrix1D< double > &steps)
Determine the positions to be minimize of a vector containing spherical harmonic coefficients.
Definition: forward_zernike_images_priors.cpp:929
Matrix1D::adaptForNumericalRecipes
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331
ProgForwardZernikeImagesPriors::Ifilteredp
std::vector< Image< double > > Ifilteredp
Definition: forward_zernike_images_priors.h:108
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
STARTINGZ
#define STARTINGZ(v)
Definition: multidim_array_base.h:67
ProgForwardZernikeImagesPriors::weightsInterpolation3D
Matrix1D< double > weightsInterpolation3D(double x, double y, double z)
Definition: forward_zernike_images_priors.cpp:1045
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
XmippMetadataProgram::produces_a_metadata
bool produces_a_metadata
Indicate that the unique final output file is a Metadata.
Definition: xmipp_metadata_program.h:88
a
doublereal * a
Definition: numerical_recipes.cpp:2230
ProgForwardZernikeImagesPriors::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_priors.cpp:937
ProgForwardZernikeImagesPriors::updateCTFImage
void updateCTFImage(double defocusU, double defocusV, double angle)
Definition: forward_zernike_images_priors.cpp:960
LOWPASS
#define LOWPASS
Definition: fourier_filter.h:77
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305
ir
int ir
Definition: image_find_center.cpp:57
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673
ProgForwardZernikeImagesPriors::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_priors.cpp:969
Mask::mode
int mode
Definition: mask.h:407
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710
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