Xmipp  v3.23.11-Nereus
ml_align2d.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Sjors Scheres scheres@cnb.csic.es (2007)
4  *
5  * Unidad de Bioinformatica del Centro Nacional de Biotecnologia , CSIC
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20  * 02111-1307 USA
21  *
22  * All comments concerning this program package may be sent to the
23  * e-mail address 'xmipp@cnb.uam.es'
24  ***************************************************************************/
25 #include "ml_align2d.h"
26 #include "core/metadata_sql.h"
27 #include "core/transformations.h"
28 //#define DEBUG_JM
29 
30 //Mutex for each thread update sums
31 pthread_mutex_t update_mutex =
32  PTHREAD_MUTEX_INITIALIZER;
33 //Mutex for each thread get next refno
34 pthread_mutex_t refno_mutex =
35  PTHREAD_MUTEX_INITIALIZER;
36 
37 // Constructor
38 ProgML2D::ProgML2D()
39 {
40  do_ML3D = false;
41  refs_per_class = 1;
42 }
43 
44 
45 void ProgML2D::defineParams()
46 {
47  addUsageLine("Perform (multi-reference) 2D-alignment using a maximum-likelihood (ML) target function.");
48  addUsageLine("+Our recommended way of performing ML alignment is to introduce as little bias in the initial reference(s) as possible.");
49  addUsageLine("+This can be done by calculting average images of random subsets of the (unaligned!) input experimental images, using the --nref option.");
50  addUsageLine("+Note that the estimates for the standard deviation in the noise and in the origin offsets are re-estimated every iteration,");
51  addUsageLine("+so that the initial values should not matter too much, as long as they are \"reasonable\". For Xmipp-normalized images,");
52  addUsageLine("+the standard deviation in the noise can be assumed to be 1. For reasonably centered particles the default value of 3 for");
53  addUsageLine("+the offsets should do the job.");
54  addUsageLine("+");
55  addUsageLine("+The output of the program consists of the refined reference images (weighted averages over all experimental images).");
56  addUsageLine("+The experimental images are not altered at all. In terms of the ML approach, optimal transformations and references for");
57  addUsageLine("+each image do not play the same role as in the conventional cross-correlation (or least-sqaures) approach. This program");
58  addUsageLine("+can also be used for reference-free 2D-alignment using only a single reference: just supply =--nref 1= .");
59  addUsageLine("+Although the calculations can be rather time-consuming (especially for many, large experimental images and a large number of references),");
60  addUsageLine("+we strongly recommend to let the calculations converge. In our experience this takes in the order of 10-100 iterations, depending on the");
61  addUsageLine("+number images, the amount of noise, etc. The default stopping criterium has yielded satisfactory results in our experience. A parallel ");
62  addUsageLine("+version of this program has been implemented.");
63  addSeeAlsoLine("mpi_ml_align2d");
64 
65  defineBasicParams(this);
66 
67  defineAdditionalParams(this, "==+ Additional options ==");
68  defineHiddenParams(this);
69 
70  addExampleLine("A typical use of this program is:", false);
71  addExampleLine("xmipp_ml_align2d -i input/images_some.stk --ref input/seeds2.stk --oroot output/ml2d --fast --mirror");
72 }
73 
74 // Read arguments ==========================================================
75 void ProgML2D::readParams()
76 {
77  // Generate new command line for restart procedure
78 
79  cline = "";
80  int argc2 = 0;
81  char ** argv2 = nullptr;
82 
83  double restart_noise=0, restart_offset=0;
84  FileName restart_imgmd, restart_refmd;
85  int restart_iter=0, restart_seed=0;
86 
87  if (!do_ML3D && checkParam("--restart"))
88  {
89  //TODO-------- Think later -----------
90  REPORT_ERROR(ERR_NOT_IMPLEMENTED, "Not implemented restart in ml2d for now");
91  // do_restart = true;
92  // MetaData MDrestart;
93  // char *copy = NULL;
94  //
95  // MDrestart.read(getParam("-restart"));
96  // cline = MDrestart.getComment();
97  // MDrestart.getValue(MDL_SIGMANOISE, restart_noise);
98  // MDrestart.getValue(MDL_SIGMAOFFSET, restart_offset);
99  // MDrestart.getValue(MDL_IMGMD, restart_imgmd);
100  // MDrestart.getValue(MDL_REFMD, restart_refmd);
101  // MDrestart.getValue(MDL_ITER, restart_iter);
102  // MDrestart.getValue(MDL_RANDOMSEED, restart_seed);
103  // generateCommandLine(cline, argc2, argv2, copy);
104  // //Take argument from the restarting command line
105  // progDef->read(argc2, argv2);
106  }
107  else
108  {
109  // no restart, just copy argc to argc2 and argv to argv2
110  do_restart = false;
111  for (int i = 1; i < argc; i++)
112  cline = cline + (String) argv[i] + " ";
113  }
114 
115  factor_nref = getIntParam("--nref");
116  fn_ref = getParam("--ref");
117  fn_img = getParam("-i");
118  fn_root = getParam("--oroot");
119  psi_step = getDoubleParam("--psi_step");
120  Niter = getIntParam("--iter");
121  //Fixed values for the 3D case
122  if (do_ML3D)
123  {
124  istart = 1;
125  fast_mode = save_mem2 = do_mirror = true;
126  }
127  else
128  {
129  //istart = getIntParam("--istart");
130  istart = 1;
131  do_mirror = checkParam("--mirror");
132  save_mem2 = checkParam("--save_memB");
133  fast_mode = checkParam("--fast");
134  }
135  model.sigma_noise = getDoubleParam("--noise");
136  model.sigma_offset = getDoubleParam("--offset");
137 
138  eps = getDoubleParam("--eps");
139  fn_frac = getParam("--frac");
140  fix_fractions = checkParam("--fix_fractions");
141  fix_sigma_offset = checkParam("--fix_sigma_offset");
142  fix_sigma_noise = checkParam("--fix_sigma_noise");
143 
144  C_fast = getDoubleParam("-C");
145  save_mem1 = checkParam("--save_memA");
146 
147  zero_offsets = checkParam("--zero_offsets");
148  model.do_student = checkParam("--student");
149  df = getDoubleParam("--student");
150  model.do_norm = checkParam("--norm");
151 
152  // Number of threads
153  threads = getIntParam("--thr");
154  //testing the thread load in refno
155  refno_load_param = getIntParam("--load");
156  // Hidden arguments
157  fn_scratch = getParameter(argc2, (const char **)argv2, "--scratch", "");
158  debug = getIntParam("--debug");
159  model.do_student_sigma_trick = !checkParam("--no_sigma_trick");
160  trymindiff_factor = getDoubleParam("--trymindiff_factor");
161  // Random seed to use for image randomization and creation
162  // of initial references and blocks order
163  // could be passed for restart or for debugging
164  seed = getIntParam("--random_seed");
165 
166  if (seed == -1)
167  seed = time(nullptr);
168  else if (!do_restart && verbose)
169  std::cerr << "WARNING: *** Using a non random seed and not in restarting ***" <<std::endl;
170 
171  // Only for interaction with refine3d:
172  search_rot = getDoubleParam("--search_rot");
173  //IEM stuff
174  blocks = getIntParam("--iem");
175 
176  // Now reset some stuff for restart
177  if (do_restart)
178  {
179  fn_img = restart_imgmd;
180  fn_ref = restart_refmd;
181  model.n_ref = 0; // Just to be sure (not strictly necessary)
182  model.sigma_noise = restart_noise;
183  model.sigma_offset = restart_offset;
184  seed = restart_seed;
185  istart = restart_iter + 1;
186  factor_nref = 1;
187  }
188 
189  no_iem = checkParam("--no_iem");
190 
191  //std::cerr << "DEBUG_JM: exiting after readParams..." <<std::endl;
192  //exit(1);
193 }
194 
195 // Show ====================================================================
196 void ProgML2D::show()
197 {
198  if (verbose)
199  {
200  // To screen
201  if (!do_ML3D)
202  {
203  std::cout
204  << " -----------------------------------------------------------------" << std::endl
205  << " | Read more about this program in the following publications: |" << std::endl
206  << " | Scheres ea. (2005) J.Mol.Biol. 348(1), 139-49 |" << std::endl
207  << " | Scheres ea. (2005) Bioinform. 21(suppl.2), ii243-4 (-fast) |" << std::endl
208  << " | |"<< std::endl
209  << " | *** Please cite them if this program is of use to you! *** |"<< std::endl
210  << " -----------------------------------------------------------------"<< std::endl;
211  }
212 
213  std::cout << "--> Maximum-likelihood multi-reference refinement " << std::endl
214  << formatString(" Input images : %s (%lu)\n", fn_img.c_str(), nr_images_global);
215 
216  if (fn_ref != "")
217  std::cout << formatString(" Reference image(s) : %s (%d)\n", fn_ref.c_str(), model.n_ref);
218  if (factor_nref > 1)
219  std::cout << " Reference expanding factor : " << factor_nref << std::endl;
220  std::cout << " Number of references: : " << model.n_ref * factor_nref << std::endl;
221 
222  std::cout
223  << " Output rootname : " << fn_root << std::endl
224  << " Stopping criterium : " << eps << std::endl
225  << " Initial sigma noise : " << model.sigma_noise << std::endl
226  << " Initial sigma offset : " << model.sigma_offset << std::endl
227  << " Psi sampling interval : " << psi_step << std::endl
228  << " Check mirrors : " << (do_mirror ? "true" : "false") << std::endl;
229 
230  if (!fn_frac.empty())
231  std::cout << " Initial model fractions : " << fn_frac << std::endl;
232 
233  if (fast_mode)
234  {
235  std::cout << " -> Use fast, reduced search-space approach with C = " << C_fast << std::endl;
236  if (zero_offsets)
237  std::cout << " + Start from all-zero translations" << std::endl;
238  }
239 
240  if (search_rot < 180.)
241  std::cout << formatString(" + Limit orientational search to +/- %f degrees", search_rot) << std::endl;
242 
243  if (save_mem1)
244  std::cout << " -> Save_memory A: recalculate real-space rotations in -fast" << std::endl;
245 
246  if (save_mem2)
247  std::cout << " -> Save_memory B: limit translations to 3 sigma_offset " << std::endl;
248 
249  if (fix_fractions)
250  std::cout << " -> Do not update estimates of model fractions." << std::endl;
251 
252  if (fix_sigma_offset)
253  std::cout << " -> Do not update sigma-estimate of origin offsets." << std::endl;
254 
255  if (fix_sigma_noise)
256  std::cout << " -> Do not update sigma-estimate of noise." << std::endl;
257 
258  if (model.do_student)
259  {
260  std::cout << " -> Use t-student distribution with df = " << df << std::endl;
261 
262  if (model.do_student_sigma_trick)
263  std::cout << " -> Use sigma-trick for t-student distributions" << std::endl;
264  }
265 
266  if (model.do_norm)
267  std::cout << " -> Refine normalization for each experimental image" << std::endl;
268 
269  if (threads > 1)
270  std::cout << " -> Using " << threads << " parallel threads" << std::endl;
271 
272  if (blocks > 1)
273  std::cout << " -> Doing IEM with " << blocks << " blocks" <<std::endl;
274 
275  std::cout << " -----------------------------------------------------------------" << std::endl;
276 
277  }
278 
279 }
280 
281 void ProgML2D::printModel(const String &msg, const ModelML2D & model)
282 {
283  std::cerr << "================> " << msg << std::endl;
284  model.print();
285 }
286 
287 
288 //#define LOG(str)
289 // Trying to merge produceSideInfo 1 y 2
290 void ProgML2D::produceSideInfo()
291 {
292  LOG(" ProgML2D::produceSideInfo: start");
293  // Read selfile with experimental images
294  // and set some global variables
295  LOG(" ProgML2D::produceSideInfo: reading MDimg");
296  MDimg.read(fn_img);
297  // Remove disabled images
298  MDimg.removeDisabled();
299  nr_images_global = MDimg.size();
300  // By default set myFirst and myLast equal to 0 and N
301  // respectively, this should be changed when using MPI
302  // by calling setWorkingImages before produceSideInfo2
303  myFirstImg = 0;
304  myLastImg = nr_images_global - 1;
305 
306  //Initialize blocks
307  current_block = 0;
308 
309  // Create a vector of objectIDs, which may be randomized later on
310  MDimg.findObjects(img_id);
311  // Get original image size
312  size_t idum, idumLong;
313  LOG(" ProgML2D::produceSideInfo: setting dimensions");
314  getImageSize(MDimg, dim, idum, idum, idumLong);
315  model.dim = dim;
316  hdim = dim / 2;
317  dim2 = dim * dim;
318  ddim2 = (double) dim2;
319  double sigma_noise2 = model.sigma_noise * model.sigma_noise;
320 
321  if (model.do_student)
322  {
323  df2 = -(df + ddim2) / 2.;
324  dfsigma2 = df * sigma_noise2;
325  }
326 
327  if (fn_ref.empty())
328  {
329  //generate an initial reference just by averaging the experimental images
330  LOG(" ProgML2D::produceSideInfo: generate an initial reference just by averaging the experimental images");
331  FileName fn_tmp;
332  Image<double> img, avg(dim, dim);
333  avg().initZeros();
334  avg().setXmippOrigin();
335 
336  for (size_t objId : MDimg.ids())
337  {
338  MDimg.getValue(MDL_IMAGE, fn_tmp, objId);
339  img.read(fn_tmp);
340  img().setXmippOrigin();
341  avg() += img();
342  }
343 
344  avg() /= nr_images_global;
345  model.setNRef(1);
346  model.Iref[0] = avg;
347  fn_ref = fn_root + "_images_average.xmp";
348  avg.write(fn_ref);
349  }
350 
351  // Print some output to screen
352  LOG(" ProgML2D::produceSideInfo show");
353  show();
354  LOG(" ProgML2D::produceSideInfo end");
355 }
356 
357 void ProgML2D::produceSideInfo2()
358 {
359  Image<double> img;
360  FileName fn_tmp;
361 
362  // Read in all reference images in memory
363  MDref.read(fn_ref);
364  MDref.removeDisabled();
365 
366  model.setNRef(MDref.size());
367  int refno = 0;
368  double fraction = (double) 1./ model.n_ref;
369  double weight = fraction * (double) nr_images_global;
370  double rot = 0., tilt = 0.;
371 
372  for (size_t objId : MDref.ids())
373  {
374  MDref.getValue(MDL_IMAGE, fn_tmp, objId);
375  MDref.getValue(MDL_ANGLE_ROT, rot, objId);
376  MDref.getValue(MDL_ANGLE_TILT, rot, objId);
377  img.read(fn_tmp);
378  img.setEulerAngles(rot, tilt, 0.);
379  img().setXmippOrigin();
380  img.setWeight(weight);
381  model.Iref[refno] = img;
382  //model.WsumMref[refno] = img;
383  // Default start is all equal model fractions
384  model.alpha_k[refno] = fraction;
385  //model.WsumMref[refno]() *= weight;
386  // Default start is half-half mirrored images
387  model.mirror_fraction[refno] = (do_mirror ? 0.5 : 0.);
388  ++refno;
389  }
390 
391  setNumberOfLocalImages();
392  // prepare masks for rotated references
393  mask.resize(dim, dim);
394  mask.setXmippOrigin();
395  BinaryCircularMask(mask, hdim, INNER_MASK);
396  omask.resize(dim, dim);
397  omask.setXmippOrigin();
398  BinaryCircularMask(omask, hdim, OUTSIDE_MASK);
399 
400  // Construct matrices for 0, 90, 180 & 270 degree flipping and mirrors
401  initSamplingStuff();
402 
403  // Set sigdim, i.e. the number of pixels that will be considered in the translations
404  sigdim = 2 * (size_t)ceil(model.sigma_offset * (save_mem2 ? 3 : 6));
405  ++sigdim; // (to get uneven number)
406  sigdim = XMIPP_MIN(dim, sigdim);
407 
408  //Some vectors and matrixes initialization
409  int num_output_refs = model.n_ref * factor_nref;
410  //std::cerr << "DEBUG_JM: num_output_refs: " << num_output_refs << std::endl;
411  refw.resize(num_output_refs);
412  refw2.resize(num_output_refs);
413  refwsc2.resize(num_output_refs);
414  refw_mirror.resize(num_output_refs);
415  sumw_refpsi.resize(num_output_refs * nr_psi);
416  A2.resize(num_output_refs);
417  fref.resize(num_output_refs * nr_psi);
418  mref.resize(num_output_refs * nr_psi);
419  wsum_Mref.resize(num_output_refs);
420  mysumimgs.resize(num_output_refs * nr_psi);
421  Iold.resize(num_output_refs);
422 
423  if (fast_mode)
424  {
425  int mysize = num_output_refs * (do_mirror ? 2 : 1);
426  //if (do_mirror)
427  // mysize *= 2;
428  ioptx_ref.resize(mysize);
429  iopty_ref.resize(mysize);
430  ioptflip_ref.resize(mysize);
431  }
432 
433  randomizeImagesOrder();
434 
435  // Initialize trymindiff for all images
436  imgs_optrefno.clear();
437  imgs_trymindiff.clear();
438  imgs_scale.clear();
439  imgs_bgmean.clear();
440  imgs_offsets.clear();
441  imgs_oldphi.clear();
442  imgs_oldtheta.clear();
443  model.scale.clear();
444 
445  if (model.do_norm)
446  {
447  average_scale = 1.;
448  for (int refno = 0; refno < model.n_ref; refno++)
449  {
450  model.scale.push_back(1.);
451  }
452  }
453 
454  // Initialize imgs_offsets vectors
455  std::vector<double> Vdum;
456  double offx = (zero_offsets ? 0. : -999);
457  int idum = (do_mirror ? 4 : 2) * factor_nref * model.n_ref;
458 
459  FOR_ALL_LOCAL_IMAGES()
460  {
461  imgs_optrefno.push_back(0);
462  imgs_trymindiff.push_back(-1.);
463  imgs_offsets.push_back(Vdum);
464  for (int refno = 0; refno < idum; refno++)
465  {
466  imgs_offsets[IMG_LOCAL_INDEX].push_back(offx);
467  }
468  if (model.do_norm)
469  {
470  // Initialize scale and bgmean for all images
471  // (for now initialize to 1 and 0, below also include doc)
472  imgs_bgmean.push_back(0.);
473  imgs_scale.push_back(1.);
474  }
475  if (limit_rot)
476  {
477  // For limited orientational search: initialize imgs_oldphi & imgs_oldtheta to -999.
478  imgs_oldphi.push_back(-999.);
479  imgs_oldtheta.push_back(-999.);
480  }
481  }
482 
483  //TODO: THINK RESTART LATER, NOW NOT WORKING
484  if (do_restart)
485  {
486  REPORT_ERROR(ERR_NOT_IMPLEMENTED, "Restart option not implement yet in ml2d");
487 
488  // Read optimal image-parameters
489  // FOR_ALL_LOCAL_IMAGES()
490  // {
491  // if (limit_rot)
492  // {
493  // MDimg.getValue(MDL_ANGLE_ROT, imgs_oldphi[IMG_LOCAL_INDEX]);
494  // MDimg.getValue(MDL_ANGLE_TILT, imgs_oldtheta[IMG_LOCAL_INDEX]);
495  // }
496  // if (zero_offsets)
497  // {
498  // idum = (do_mirror ? 2 : 1) * model.n_ref;
499  // double xx, yy;
500  // MDimg.getValue(MDL_SHIFT_X, xx);
501  // MDimg.getValue(MDL_SHIFT_X, yy);
502  // for (int refno = 0; refno < idum; refno++)
503  // {
504  // imgs_offsets[IMG_LOCAL_INDEX][2 * refno] = xx;
505  // imgs_offsets[IMG_LOCAL_INDEX][2 * refno + 1] = yy;
506  // }
507  // }
508  // if (model.do_norm)
509  // {
510  // MDimg.getValue(MDL_BGMEAN, imgs_bgmean[IMG_LOCAL_INDEX]);
511  // MDimg.getValue(MDL_INTSCALE, imgs_scale[IMG_LOCAL_INDEX]);
512  // }
513  // }
514  // // read Model parameters
515  // int refno = 0;
516  // FOR_ALL_OBJECTS_IN_METADATA(MDref)
517  // {
518  // MDref.getValue(MDL_MODELFRAC, model.alpha_k[refno]);
519  // if (do_mirror)
520  // MDref.getValue(MDL_MIRRORFRAC,
521  // model.mirror_fraction[refno]);
522  // if (model.do_norm)
523  // MDref.getValue(MDL_INTSCALE, model.scale[refno]);
524  // refno++;
525  // }
526 
527  }
528  // Call the first iteration 0 if generating initial references from random subsets
529  else if (factor_nref > 1)
530  {
531  istart = SPECIAL_ITER;
532 
533  //Expand MDref because it will be re-used in writeOutputFiles
534  MetaDataDb MDaux(MDref);
535 
536  for (int group = 1; group < factor_nref; ++group)
537  {
538  for (size_t objId : MDaux.ids())
539  MDaux.setValue(MDL_REF3D, group + 1, objId);
540  MDref.unionAll(MDaux);
541  }
542 
543  //std::cerr << "DEBUG_JM: MDref: " << std::endl;
544  //MDref.write(std::cerr);
545  }
546 
547  //--------Setup for Docfile -----------
548  docfiledata.resize(nr_images_local, DATALINELENGTH);
549 
550 }//close function produceSideInfo
551 
552 // Calculate probability density function of all in-plane transformations phi
553 void ProgML2D::calculatePdfInplane()
554 {
555 #ifdef DEBUG
556  std::cerr << "Entering calculatePdfInplane" <<std::endl;
557 #endif
558 
559  double x, y, r2, pdfpix, sum;
560  P_phi.resize(dim, dim);
561  P_phi.setXmippOrigin();
562  Mr2.resize(dim, dim);
563  Mr2.setXmippOrigin();
564 
565  sum = 0.;
566  double _2sigma_offset2 = 2 * model.sigma_offset * model.sigma_offset;
567  double _sigma_offset_denominator = 2 * PI * model.sigma_offset * model.sigma_offset * nr_psi * nr_nomirror_flips;
568 
569  FOR_ALL_ELEMENTS_IN_ARRAY2D(P_phi)
570  {
571  x = (double) j;
572  y = (double) i;
573  r2 = x * x + y * y;
574 
575  if (model.sigma_offset > 0.)
576  {
577  pdfpix = exp(-r2 / _2sigma_offset2);
578  pdfpix /= _sigma_offset_denominator;
579  }
580  else
581  {
582  if (j == 0 && i == 0)
583  pdfpix = 1.;
584  else
585  pdfpix = 0.;
586  }
587 
588  A2D_ELEM(P_phi, i, j) = pdfpix;
589  A2D_ELEM(Mr2, i, j) = (double) r2;
590  sum += pdfpix;
591  }
592 
593  // Normalization
594  P_phi /= sum;
595 #ifdef DEBUG
596 
597  std::cerr << "Leaving calculatePdfInplane" <<std::endl;
598 #endif
599 
600 }
601 
602 // Rotate reference for all models and rotations and fill Fref vectors =============
603 void ProgML2D::rotateReference()
604 {
605 #ifdef DEBUG
606  std::cerr<<"entering rotateReference"<<std::endl;
607 #endif
608 
609  awakeThreads(TH_RR_REFNO, 0, refno_load_param);
610 
611 #ifdef DEBUG
612 
613  std::cerr<<"leaving rotateReference"<<std::endl;
614 #endif
615 }
616 
617 // Collect all rotations and sum to update Iref() for all models ==========
618 void ProgML2D::reverseRotateReference()
619 {
620 
621 #ifdef DEBUG
622  std::cerr<<"entering reverseRotateReference"<<std::endl;
623 #endif
624 
625  awakeThreads(TH_RRR_REFNO, 0, refno_load_param);
626 
627 #ifdef DEBUG
628 
629  std::cerr<<"leaving reverseRotateReference"<<std::endl;
630 #endif
631 
632 }
633 
634 void ProgML2D::preselectLimitedDirections(double &phi, double &theta)
635 {
636 
637  double phi_ref, theta_ref, angle, angle2;
638  Matrix1D<double> u, v;
639 
640  pdf_directions.clear();
641  pdf_directions.resize(model.n_ref);
642 
643  for (int refno = 0; refno < model.n_ref; refno++)
644  {
645  if (!limit_rot || (phi == -999. && theta == -999.))
646  pdf_directions[refno] = 1.;
647  else
648  {
649  phi_ref = model.Iref[refno].rot();
650  theta_ref = model.Iref[refno].tilt();
651  Euler_direction(phi, theta, 0., u);
652  Euler_direction(phi_ref, theta_ref, 0., v);
653  u.selfNormalize();
654  v.selfNormalize();
655  angle = RAD2DEG(acos(dotProduct(u, v)));
656  angle = fabs(realWRAP(angle, -180, 180));
657  // also check mirror
658  angle2 = 180. + angle;
659  angle2 = fabs(realWRAP(angle2, -180, 180));
660  angle = XMIPP_MIN(angle, angle2);
661 
662  if (fabs(angle) > search_rot)
663  pdf_directions[refno] = 0.;
664  else
665  pdf_directions[refno] = 1.;
666  }
667  }
668 
669 }
670 
671 // Pre-selection of significant refno and ipsi, based on current optimal translation =======
672 
673 
674 void ProgML2D::preselectFastSignificant()
675 {
676 
677 #ifdef DEBUG
678  std::cerr<<"entering preselectFastSignificant"<<std::endl;
679 #endif
680 
681  // Initialize Msignificant to all zeros
682  // TODO: check whether this is strictly necessary? Probably not...
683  Msignificant.initZeros();
684  pfs_mindiff = 99.e99;
685  pfs_maxweight.resizeNoCopy((do_mirror ? 2 : 1), model.n_ref);
686  pfs_maxweight.initConstant(-99.e99);
687  pfs_weight.resizeNoCopy(model.n_ref, nr_psi * nr_flip);
688  pfs_weight.initZeros();
689  pfs_count = threads;
690  awakeThreads(TH_PFS_REFNO, 0, refno_load_param);
691 
692 #ifdef DEBUG
693 
694  std::cerr<<"leaving preselectFastSignificant"<<std::endl;
695 #endif
696 }
697 
698 // Maximum Likelihood calculation for one image ============================================
699 // Integration over all translation, given model and in-plane rotation
700 void ProgML2D::expectationSingleImage(Matrix1D<double> &opt_offsets)
701 {
702 #ifdef TIMING
703  timer.tic(ESI_E1);
704 #endif
705 
706  MultidimArray<double> Maux, Mweight;
707  MultidimArray<std::complex<double> > Faux;
708  double my_mindiff;
709  bool is_ok_trymindiff = false;
710  double sigma_noise2 = model.sigma_noise * model.sigma_noise;
711  FourierTransformer local_transformer;
712  ioptx = iopty = 0;
713 
714  // Update sigdim, i.e. the number of pixels that will be considered in the translations
715  sigdim = 2 * CEIL(XMIPP_MAX(1,model.sigma_offset) * (save_mem2 ? 3 : 6));
716  sigdim++; // (to get uneven number)
717  sigdim = XMIPP_MIN(dim, sigdim);
718 
719  // Setup matrices
720  Maux.resize(dim, dim);
721  Maux.setXmippOrigin();
722  Mweight.initZeros(sigdim, sigdim);
723  Mweight.setXmippOrigin();
724 
725  if (!model.do_norm)
726  opt_scale = 1.;
727 
728  // precalculate all flipped versions of the image
729  Fimg_flip.clear();
730 
731  for (size_t iflip = 0; iflip < nr_flip; iflip++)
732  {
733  Maux.setXmippOrigin();
734  applyGeometry(xmipp_transformation::LINEAR, Maux, Mimg, F[iflip], xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
735  local_transformer.FourierTransform(Maux, Faux, false);
736 
737  if (model.do_norm)
738  dAij(Faux,0,0) -= bgmean;
739 
740  Fimg_flip.push_back(Faux);
741 
742  }
743 
744  // The real stuff: loop over all references, rotations and translations
745  int redo_counter = 0;
746 
747 #ifdef TIMING
748 
749  timer.toc(ESI_E1);
750 
751  //timer.tic("WHILE_");
752 #endif
753 
754  while (!is_ok_trymindiff)
755  {
756  // Initialize mindiff, weighted sums and maxweights
757  mindiff = 99.e99;
758  wsum_corr = wsum_offset = wsum_sc = wsum_sc2 = 0.;
759  maxweight = maxweight2 = sum_refw = sum_refw2 = 0.;
760 
761  awakeThreads(TH_ESI_REFNO, opt_refno, refno_load_param);
762 
763  // Now check whether our trymindiff was OK.
764  // The limit of the exp-function lies around
765  // exp(700)=1.01423e+304, exp(800)=inf; exp(-700) = 9.85968e-305; exp(-88) = 0
766  // Use 500 to be on the save side?
767 
768  if (ABS((mindiff - trymindiff) / sigma_noise2) > 500.)
769  //force always redo to use real mindiff for check about LL problem
770  //if (redo_counter==0)
771  {
772  // Re-do whole calculation now with the real mindiff
773  trymindiff = mindiff;
774  redo_counter++;
775  // On iteration 0 images that will go to references other than first
776  // will store optimus references but not yet expanded number of references
777  if (iter == 0)
778  opt_refno = (opt_refno % model.n_ref);
779 
780  // Never re-do more than once!
781  if (redo_counter > 1)
782  REPORT_ERROR(ERR_VALUE_INCORRECT, "ml_align2d BUG% redo_counter > 1");
783  }
784  else
785  {
786  is_ok_trymindiff = true;
787  my_mindiff = trymindiff;
788  trymindiff = mindiff;
789  }
790 
791  }//close while
792 
793  fracweight = maxweight / sum_refw;
794 
795  wsum_sc /= sum_refw;
796 
797  wsum_sc2 /= sum_refw;
798 
799  // Calculate optimal transformation parameters
800  opt_psi = -psi_step * (iopt_flip * nr_psi + iopt_psi) - SMALLANGLE;
801 
802  opt_offsets(0) = -(double) ioptx * MAT_ELEM(F[iopt_flip], 0, 0)
803  - (double) iopty * MAT_ELEM(F[iopt_flip], 0, 1);
804 
805  opt_offsets(1) = -(double) ioptx * MAT_ELEM(F[iopt_flip], 1, 0)
806  - (double) iopty * MAT_ELEM(F[iopt_flip], 1, 1);
807 
808  // Update normalization parameters
809  if (model.do_norm)
810  {
811  // 1. Calculate optimal setting of Mimg
812  MultidimArray<double> Maux2 = Mimg;
813  selfTranslate(xmipp_transformation::LINEAR, Maux2, opt_offsets, true);
814  selfApplyGeometry(xmipp_transformation::LINEAR, Maux2, F[iopt_flip], xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
815  // 2. Calculate optimal setting of Mref
816  int refnoipsi = (opt_refno % model.n_ref) * nr_psi + iopt_psi;
817  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Faux)
818  {
819  dAij(Faux,i,j) = conj(dAij(fref[refnoipsi],i,j));
820  dAij(Faux,i,j) *= opt_scale;
821  }
822 
823  // Still take input from Faux and leave output in Maux
824  local_transformer.inverseFourierTransform();
825  Maux2 = Maux2 - Maux;
826 
827  if (debug == 12)
828  {
829  std::cout << std::endl;
830  std::cout << "scale= " << opt_scale << " changes to " << wsum_sc
831  / wsum_sc2 << std::endl;
832  std::cout << "bgmean= " << bgmean << " changes to "
833  << Maux2.computeAvg() << std::endl;
834  }
835 
836  // non-ML update of bgmean (this is much cheaper than true-ML update...)
837  old_bgmean = bgmean;
838 
839  bgmean = Maux2.computeAvg();
840 
841  // ML-update of opt_scale
842  opt_scale = wsum_sc / wsum_sc2;
843  }
844 
845 #ifdef TIMING
846  timer.toc(ESI_E5);
847 
848  timer.tic(ESI_E6TH);
849 
850 #endif
851  // Update all global weighted sums after division by sum_refw
852  wsum_sigma_noise += (2 * wsum_corr / sum_refw);
853 
854  wsum_sigma_offset += (wsum_offset / sum_refw);
855 
856  sumfracweight += fracweight;
857 
858 
859  // std::cerr << "-------------------- wsum_corr: " << wsum_corr << std::endl;
860  // std::cerr << "-------------------- sum_refw: " << sum_refw << std::endl;
861  // std::cerr << std::endl << "-------------------- wsum_sigma_offset: " << wsum_sigma_offset << std::endl;
862  // std::cerr << "-------------------- wsum_sigma_noise: " << wsum_sigma_noise << std::endl << std::endl;
863  // if (wsum_sigma_noise < 0)
864  // {
865  // std::cerr << "Negative wsum_sigma_noise....exiting" <<std::endl;
866  // exit(1);
867  // }
868 
869  awakeThreads(TH_ESI_UPDATE_REFNO, 0, refno_load_param);
870 
871  if (!model.do_student)
872  // 1st term: log(refw_i)
873  // 2nd term: for subtracting mindiff
874  // 3rd term: for (sqrt(2pi)*sigma_noise)^-1 term in formula (12) Sigworth (1998)
875  dLL = log(sum_refw) - my_mindiff / sigma_noise2 - ddim2 * log(sqrt(2.
876  * PI * sigma_noise2));
877  else
878  // 1st term: log(refw_i)
879  // 2nd term: for dividing by (1 + 2. * mindiff/dfsigma2)^df2
880  // 3rd term: for sigma-dependent normalization term in t-student distribution
881  // 4th&5th terms: gamma functions in t-distribution
882  dLL = log(sum_refw) + df2 * log(1. + (2. * my_mindiff / dfsigma2))
883  - ddim2 * log(sqrt(PI * df * sigma_noise2)) + gammln(-df2)
884  - gammln(df / 2.);
885 
886  //#define DEBUG_JM1
887 #ifdef DEBUG_JM1
888 
889  //if (iter>1)
890  {
891  //static std::ios::openmode mode = std::ios::out|std::ios::trunc;
892  //std::ofstream log;
893  //FileName fn = fn_root + (current_image % 2 == 0 ? "_images.even" : "_images.odd");
894  //log.open(fn.c_str(), mode);
895  //std::cerr << " IMAGE " << current_image << "----------------------------->>>" << std::endl;
896  std::cerr << "----------------------------->>>" << std::endl;
897  std::cerr << " dLL: " << dLL << std::endl;
898  std::cerr << " sum_refw: " << sum_refw << std::endl;
899  std::cerr << " my_mindiff: " << my_mindiff << std::endl;
900  std::cerr << " sigma_noise2: " << sigma_noise2 << std::endl;
901  std::cerr << " ddim2: " << ddim2 << std::endl;
902  //std::cerr << " dfsigma2: " << dfsigma2 << std::endl;
903  std::cerr << " wsum_corr: " << wsum_corr << std::endl;
904  std::cerr << " wsum_offset: " << wsum_offset << std::endl;
905  // std::cerr << " refw: ";
906  // for (int refno = 0; refno < model.n_ref; ++refno)
907  // std::cerr << std::setw(15) << refw[refno];
908  // std::cerr<< std::endl;
909  // std::cerr << " refw_mirror: ";
910  // for (int refno = 0; refno < model.n_ref; ++refno)
911  // std::cerr << std::setw(15) << refw_mirror[refno];
912  // std::cerr<< std::endl;
913  // log.close();
914  // mode = std::ios::out|std::ios::app;
915  }
916 #endif
917 #undef DEBUG_JM1
918 
919  LL += dLL;
920 
921 #ifdef TIMING
922 
923  timer.toc(ESI_E6TH);
924 
925 #endif
926 }//close function expectationSingleImage
927 
929 void ProgML2D::createThreads()
930 {
931 
932  //Initialize some variables for using for threads
933 
934  th_ids = new pthread_t[threads];
935  threads_d = new structThreadTasks[threads];
936 
937  barrier_init(&barrier, threads + 1);
938  barrier_init(&barrier2, threads + 1);
939  barrier_init(&barrier3, threads);//Main thread will not wait here
940 
941  for (int i = 0; i < threads; i++)
942  {
943  threads_d[i].thread_id = i;
944  threads_d[i].prm = this;
945 
946  int result = pthread_create((th_ids + i), nullptr, doThreadsTasks,
947  (void *) (threads_d + i));
948 
949  if (result != 0)
950  {
951  REPORT_ERROR(ERR_THREADS_NOTINIT, "");
952  }
953  }
954 
955 }//close function createThreads
956 
958 void ProgML2D::destroyThreads()
959 {
960  threadTask = TH_EXIT;
961  barrier_wait(&barrier);
962  delete[] th_ids;
963  delete[] threads_d;
964 }
965 
967 void * doThreadsTasks(void * data)
968 {
969  auto * thread_data = (structThreadTasks *) data;
970 
971  ProgML2D * prm = thread_data->prm;
972 
973  barrier_t & barrier = prm->barrier;
974  barrier_t & barrier2 = prm->barrier2;
975 
976  //Loop until the threadTask become TH_EXIT
977 
978  do
979  {
980  //Wait until main threads order to start
981  //debug_print("Waiting on barrier, th ", thread_id);
982  barrier_wait(&barrier);
983 
984  //Check task to do
985 
986  switch (prm->threadTask)
987  {
988 
989  case TH_PFS_REFNO:
990  prm->doThreadPreselectFastSignificantRefno();
991  break;
992 
993  case TH_ESI_REFNO:
994  prm->doThreadExpectationSingleImageRefno();
995  break;
996 
997  case TH_ESI_UPDATE_REFNO:
998  prm->doThreadESIUpdateRefno();
999  break;
1000 
1001  case TH_RR_REFNO:
1002  prm->doThreadRotateReferenceRefno();
1003  break;
1004 
1005  case TH_RRR_REFNO:
1006  prm->doThreadReverseRotateReferenceRefno();
1007  break;
1008 
1009  case TH_EXIT:
1010  pthread_exit(nullptr);
1011  break;
1012 
1013  }
1014 
1015  barrier_wait(&barrier2);
1016 
1017  }
1018  while (1);
1019 
1020 }//close function doThreadsTasks
1021 
1022 
1026 int ProgML2D::getThreadRefnoJob(int &refno)
1027 {
1028  int load = 0;
1029 
1030  pthread_mutex_lock(&refno_mutex);
1031 
1032  if (refno_count < model.n_ref)
1033  {
1034  load = XMIPP_MIN(refno_load, model.n_ref - refno_count);
1035  refno = refno_index;
1036  refno_index = (refno_index + load) % model.n_ref;
1037  refno_count += load;
1038  }
1039 
1040  pthread_mutex_unlock(&refno_mutex);
1041 
1042  return load;
1043 }//close function getThreadRefnoJob
1044 
1046 void ProgML2D::awakeThreads(ThreadTask task, int start_refno, int load)
1047 {
1048  threadTask = task;
1049  refno_index = start_refno;
1050  refno_count = 0;
1051  refno_load = load;
1052  barrier_wait(&barrier);
1053  //Wait until done
1054  barrier_wait(&barrier2);
1055 }//close function awakeThreads
1056 
1057 
1058 void ProgML2D::doThreadRotateReferenceRefno()
1059 {
1060 #ifdef DEBUG
1061  std::cerr << "entering doThreadRotateReference " << std::endl;
1062 #endif
1063 
1064  double AA, stdAA=0., psi, dum, avg;
1065  MultidimArray<double> Maux(dim, dim);
1066  MultidimArray<std::complex<double> > Faux;
1067  FourierTransformer local_transformer;
1068  int refnoipsi;
1069 
1070  Maux.setXmippOrigin();
1071 
1072  FOR_ALL_THREAD_REFNO()
1073  {
1074  computeStats_within_binary_mask(omask, model.Iref[refno](), dum,
1075  dum, avg, dum);
1076  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1077  {
1078  refnoipsi = refno * nr_psi + ipsi;
1079  // Add arbitrary number (small_angle) to avoid 0-degree rotation (lacking interpolation)
1080  psi = (double) (ipsi * psi_max / nr_psi) + SMALLANGLE;
1081  rotate(xmipp_transformation::BSPLINE3, Maux, model.Iref[refno](), -psi, 'Z', xmipp_transformation::WRAP);
1082  apply_binary_mask(mask, Maux, Maux, avg);
1083  // Normalize the magnitude of the rotated references to 1st rot of that ref
1084  // This is necessary because interpolation due to rotation can lead to lower overall Fref
1085  // This would result in lower probabilities for those rotations
1086  AA = Maux.sum2();
1087  if (ipsi == 0)
1088  {
1089  stdAA = AA;
1090  A2[refno] = AA;
1091  }
1092 
1093  if (AA > 0)
1094  Maux *= sqrt(stdAA / AA);
1095 
1096  if (fast_mode)
1097  mref[refnoipsi] = Maux;
1098 
1099  // Do the forward FFT
1100  local_transformer.FourierTransform(Maux, Faux, false);
1101 
1102  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Faux)
1103  {
1104  dAij(Faux, i, j) = conj(dAij(Faux, i, j));
1105  }
1106 
1107  fref[refnoipsi] = Faux;
1108  }
1109 
1110  // If we don't use save_mem1 Iref[refno] is useless from here on
1111  //FIXME: Segmentation fault with blocks
1112  //if (!save_mem1)
1113  // model.Iref[refno]().resize(0, 0);
1114 
1115  }//close while
1116 
1117 }//close function doThreadRotateReferenceRefno
1118 
1119 void ProgML2D::doThreadReverseRotateReferenceRefno()
1120 {
1121  double psi, dum, avg;
1122  MultidimArray<double> Maux(dim, dim), Maux2(dim, dim), Maux3(dim, dim);
1123  MultidimArray<std::complex<double> > Faux;
1124  FourierTransformer local_transformer;
1125 
1126  Maux.setXmippOrigin();
1127  Maux2.setXmippOrigin();
1128  Maux3.setXmippOrigin();
1129 
1130  FOR_ALL_THREAD_REFNO()
1131  {
1132  Maux.initZeros();
1133  wsum_Mref[refno] = Maux;
1134  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1135  {
1136  // Add arbitrary number to avoid 0-degree rotation without interpolation effects
1137  psi = (double) (ipsi * psi_max / nr_psi) + SMALLANGLE;
1138  int refnoipsi = refno * nr_psi + ipsi;
1139  // Do the backward FFT
1140  // The construction with Faux and Maux3 should perhaps not be necessary,
1141  // but I am having irreproducible segmentation faults for the iFFT
1142  Faux = wsumimgs[refnoipsi];
1143 
1144  local_transformer.inverseFourierTransform(Faux, Maux);
1145  Maux3 = Maux;
1146  //CenterFFT(Maux3, true);
1147  centerFFT2(Maux3);
1148  computeStats_within_binary_mask(omask, Maux3, dum, dum, avg, dum);
1149  rotate(xmipp_transformation::BSPLINE3, Maux2, Maux3, psi, 'Z', xmipp_transformation::WRAP);
1150  apply_binary_mask(mask, Maux2, Maux2, avg);
1151  wsum_Mref[refno] += Maux2;
1152  }
1153 
1154  }//close while refno
1155 
1156 }//close function doThreadReverseRotateReference
1157 
1159 #define IIFLIP (imirror * nr_nomirror_flips + iflip)
1160 #define IROT (IIFLIP * nr_psi + ipsi)
1161 #define IREFMIR ()
1162 #define WEIGHT (dAij(pfs_weight, refno, IROT))
1163 #define MAX_WEIGHT (dAij(pfs_maxweight, imirror, refno))
1164 #define MSIGNIFICANT (dAij(Msignificant, refno, IROT))
1165 
1166 void ProgML2D::doThreadPreselectFastSignificantRefno()
1167 {
1168  MultidimArray<double> Mtrans, Mflip;
1169  double ropt, aux, diff, pdf, fracpdf;
1170  double A2_plus_Xi2;
1171  int irefmir;
1172  Matrix1D<double> trans(2);
1173  double local_mindiff;
1174  double sigma_noise2 = model.sigma_noise * model.sigma_noise;
1175  int nr_mirror = (do_mirror) ? 2 : 1;
1176 
1177  Mtrans.resize(dim, dim);
1178  Mtrans.setXmippOrigin();
1179  Mflip.resize(dim, dim);
1180  Mflip.setXmippOrigin();
1181 
1182  local_mindiff = 99.e99;
1183 
1184  // A. Translate image and calculate probabilities for every rotation
1185  FOR_ALL_THREAD_REFNO()
1186  {
1187  if (!limit_rot || pdf_directions[refno] > 0.)
1188  {
1189  A2_plus_Xi2 = 0.5 * (A2[refno] + Xi2);
1190  for (int imirror = 0; imirror < nr_mirror; imirror++)
1191  {
1192  irefmir = imirror * model.n_ref + refno;
1193  // Get optimal offsets
1194  trans(0) = allref_offsets[2 * irefmir];
1195  trans(1) = allref_offsets[2 * irefmir + 1];
1196  ropt = sqrt(trans(0) * trans(0) + trans(1) * trans(1));
1197  // Do not trust optimal offsets if they are larger than 3*sigma_offset:
1198  if (ropt > 3 * model.sigma_offset)
1199  {
1200  for (size_t iflip = 0; iflip < nr_nomirror_flips; iflip++)
1201  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1202  MSIGNIFICANT = 1;
1203  }
1204  else
1205  {
1206  translate(xmipp_transformation::LINEAR, Mtrans, Mimg, trans, true);
1207  for (size_t iflip = 0; iflip < nr_nomirror_flips; iflip++)
1208  {
1209  applyGeometry(xmipp_transformation::LINEAR, Mflip, Mtrans, F[IIFLIP], xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
1210  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1211  {
1212  diff = A2_plus_Xi2;
1213  MultidimArray<double> &mref_ref = mref[refno*nr_psi + ipsi];
1214  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Mflip)
1215  {
1216  diff -= DIRECT_MULTIDIM_ELEM(Mflip, n) * DIRECT_MULTIDIM_ELEM(mref_ref, n);
1217  }
1218  WEIGHT = diff;
1219  if (diff < local_mindiff)
1220  local_mindiff = diff;
1221  }
1222  }
1223  }//close else if ropt > ...
1224  }//close for imirror
1225  }//close if !limit_rot ...
1226  }//close for_all refno
1227 
1229  pthread_mutex_lock(&update_mutex);
1230  pfs_mindiff = XMIPP_MIN(pfs_mindiff, local_mindiff);
1231  refno_index = refno_count = 0;
1232  pthread_mutex_unlock(&update_mutex);
1233 
1235  barrier_wait(&barrier3);
1236  local_mindiff = pfs_mindiff;
1237 
1238  // B. Now that we have local_mindiff, calculate the weights
1239  FOR_ALL_THREAD_REFNO_NODECL()
1240  {
1241 
1242  if (!limit_rot || pdf_directions[refno] > 0.)
1243  {
1244  for (int imirror = 0; imirror < nr_mirror; imirror++)
1245  {
1246  irefmir = imirror * model.n_ref + refno;
1247  // Get optimal offsets
1248  trans(0) = allref_offsets[2 * irefmir];
1249  trans(1) = allref_offsets[2 * irefmir + 1];
1251  for (size_t iflip = 0; iflip < nr_nomirror_flips; iflip++)
1252  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1253  {
1254  if (!MSIGNIFICANT)
1255  {
1256  fracpdf = model.alpha_k[refno] *
1257  (imirror ? model.mirror_fraction[refno] : (1.- model.mirror_fraction[refno]));
1258  pdf = fracpdf * A2D_ELEM(P_phi, (int)trans(1), (int)trans(0));
1259  if (!model.do_student)
1260  {
1261  // normal distribution
1262  aux = (WEIGHT - local_mindiff) / sigma_noise2;
1263  // next line because of numerical precision of exp-function
1264  WEIGHT = (aux > 1000. ? 0. : exp(-aux) * pdf);
1265  }
1266  else
1267  {
1268  // t-student distribution
1269  aux = (dfsigma2 + 2. * WEIGHT) / (dfsigma2 + 2. * local_mindiff);
1270  WEIGHT = pow(aux, df2) * pdf;
1271  }
1272  if (WEIGHT > MAX_WEIGHT)
1273  MAX_WEIGHT = WEIGHT;
1274  }
1275  } // close ipsi
1277  for (size_t iflip = 0; iflip < nr_nomirror_flips; iflip++)
1278  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1279  if (!MSIGNIFICANT)
1280  MSIGNIFICANT = (WEIGHT >= C_fast * MAX_WEIGHT) ? 1 : 0;
1281  }//close for imirror
1282  } //endif limit_rot and pdf_directions
1283  } //end for_all refno
1284 
1285 }//close function doThreadPreselectFastSignificantRefno
1286 
1287 void ProgML2D::doThreadExpectationSingleImageRefno()
1288 {
1289  double diff;
1290  double aux, pdf, fracpdf, A2_plus_Xi2;
1291  double weight, stored_weight, weight2 = 0, my_maxweight;
1292  double my_sumweight, my_sumstoredweight, ref_scale = 1.;
1293  int irot, output_irefmir, refnoipsi, output_refnoipsi;
1294  //Some local variables to store partial sums of global sums variables
1295  double local_mindiff, local_wsum_corr, local_wsum_offset, maxw_ref;
1296  double local_wsum_sc, local_wsum_sc2, local_maxweight, local_maxweight2=0.0;
1297  double sigma_noise2 = model.sigma_noise * model.sigma_noise;
1298  int local_iopty=0, local_ioptx=0, local_iopt_psi=0, local_iopt_flip=0,
1299  local_opt_refno=0;
1300 
1301  MultidimArray<double> Maux, Mweight;
1302  MultidimArray<std::complex<double> > Faux, Fzero(dim, hdim + 1);
1303  FourierTransformer local_transformer;
1304 
1305  // Setup matrices
1306  Maux.resize(dim, dim);
1307  Maux.setXmippOrigin();
1308  Mweight.resize(sigdim, sigdim);
1309  Mweight.setXmippOrigin();
1310  Fzero.initZeros();
1311  local_transformer.setReal(Maux);
1312  local_transformer.getFourierAlias(Faux);
1313 
1314  // Start the loop over all refno at old_optrefno (=opt_refno from the previous iteration).
1315  // This will speed-up things because we will find Pmax probably right away,
1316  // and this will make the if-statement that checks SIGNIFICANT_WEIGHT_LOW
1317  // effective right from the start
1318  //std::cerr << "DEBUG_JM: doThreadExpectationSingleImageRefno: " << std::endl;
1319 
1320  if (nr_nomirror_flips == 0)
1321  {
1322  std::ostringstream msg;
1323  msg << "Division by zero: nr_nomirror_flips == 0";
1324  throw std::runtime_error(msg.str());
1325  }
1326 
1327  FOR_ALL_THREAD_REFNO()
1328  {
1329 
1330  int output_refno = mygroup * model.n_ref + refno;
1331  // std::cerr << "DEBUG_JM: refno: " << refno << std::endl;
1332  // std::cerr << "DEBUG_JM: model.n_ref: " << model.n_ref << std::endl;
1333  // std::cerr << "DEBUG_JM: output_refno: " << output_refno << std::endl;
1334  refw[output_refno] = refw2[output_refno] = refw_mirror[output_refno] = 0.;
1335  local_maxweight = -99.e99;
1336  local_mindiff = 99.e99;
1337  local_wsum_sc = local_wsum_sc2 = local_wsum_corr = local_wsum_offset = 0;
1338  // Initialize my weighted sums
1339  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1340  {
1341  output_refnoipsi = output_refno * nr_psi + ipsi;
1342  mysumimgs[output_refnoipsi] = Fzero;
1343  sumw_refpsi[output_refnoipsi] = 0.;
1344  }
1345 
1346  // This if is for limited rotation options
1347  if (!limit_rot || pdf_directions[refno] > 0.)
1348  {
1349  if (model.do_norm)
1350  ref_scale = opt_scale / model.scale[refno];
1351 
1352  A2_plus_Xi2 = 0.5 * (ref_scale * ref_scale * A2[refno] + Xi2);
1353 
1354  maxw_ref = -99.e99;
1355  for (size_t iflip = 0; iflip < nr_flip; iflip++)
1356  {
1357  if (iflip == nr_nomirror_flips)
1358  maxw_ref = -99.e99;
1359  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1360  {
1361  refnoipsi = refno * nr_psi + ipsi;
1362  output_refnoipsi = output_refno * nr_psi + ipsi;
1363  irot = iflip * nr_psi + ipsi;
1364  output_irefmir = (int)floor(iflip / nr_nomirror_flips)
1365  * factor_nref * model.n_ref + refno;
1366 
1367  //#define DEBUG_JM2
1368 #ifdef DEBUG_JM2
1369 
1370  if (iter > 1)
1371  {
1372  std::cerr << iter << " iflip, ipsi, refno, irot: " << iflip << " " << ipsi << " " << refno << " " << irot << std::endl;
1373  std::cerr << iter << " A2_plus_Xi2: " << A2_plus_Xi2 << std::endl;
1374  std::cerr << "dAij(Msignificant): " << dAij(Msignificant, refno, irot) << std::endl;
1375  }
1376 #endif
1377  // This if is the speed-up caused by the -fast options
1378  if (dAij(Msignificant, refno, irot))
1379  {
1380  if (iflip < nr_nomirror_flips)
1381  fracpdf = model.alpha_k[refno] * (1. - model.mirror_fraction[refno]);
1382  else
1383  fracpdf = model.alpha_k[refno] * model.mirror_fraction[refno];
1384 
1385  // A. Backward FFT to calculate weights in real-space
1386  //Set this references to avoid indexing inside the heavy loop
1387  //MultidimArray<std::complex<double> > & Fimg_flip_aux = Fimg_flip[iflip];
1388  Faux = Fimg_flip[iflip];
1389  MultidimArray<std::complex<double> > & fref_aux = fref[refnoipsi];
1390  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Faux)
1391  {
1392  DIRECT_MULTIDIM_ELEM(Faux, n) *= DIRECT_MULTIDIM_ELEM(fref_aux, n);
1393  }
1394  // Takes the input from Faux, and leaves the output in Maux
1395  local_transformer.inverseFourierTransform();
1396  //CenterFFT(Maux, true);
1397  centerFFT2(Maux);
1398 
1399 
1400 #ifdef DEBUG_JM2
1401 
1402  if (iter > 1)
1403  std::cerr
1404  << "Maux: " << std::endl << Maux
1405  << "Fimg_flip[iflip]" << std::endl<< Fimg_flip[iflip]
1406  << "fref[refnoipsi]" << std::endl<< fref[refnoipsi] << std::endl;
1407 #endif
1408 
1409  // B. Calculate weights for each pixel within sigdim (Mweight)
1410  my_sumweight = my_sumstoredweight = my_maxweight = 0.;
1411 
1412  FOR_ALL_ELEMENTS_IN_ARRAY2D(Mweight)
1413  {
1414  diff = A2_plus_Xi2 - ref_scale * A2D_ELEM(Maux, i, j) * ddim2;
1415  pdf = fracpdf * A2D_ELEM(P_phi, i, j);
1416 
1417 #ifdef DEBUG_JM2
1418 
1419  if (iter >= 2 && current_image == myFirstImg)
1420  std::cerr << "---------------------------------------" << std::endl
1421  << " pdf " << pdf << std::endl
1422  << " A2_plus_Xi2 " << A2_plus_Xi2 << std::endl
1423  << " A2D_ELEM(Maux, i, j) " << A2D_ELEM(Maux, i, j) << std::endl
1424  << " ref_scale " << ref_scale << std::endl
1425  << " ddim2 " << ddim2 << std::endl
1426  << "diff " << diff << std::endl
1427  << "trymindiff " << trymindiff << std::endl
1428  << "sigma_noise2 " << sigma_noise2 << std::endl;
1429 #endif
1430 
1431  if (!model.do_student)
1432  {
1433  // Normal distribution
1434  aux = (diff - trymindiff) / sigma_noise2;
1435  // next line because of numerical precision of exp-function
1436  weight = (aux > 1000.) ? 0. : exp(-aux) * pdf;
1437  //#define DEBUG_JM2
1438 #ifdef DEBUG_JM2
1439 
1440  if (iter >=2 && current_image == myFirstImg && pdf > 0)
1441  std::cerr << "aux = (diff - trymindiff) / sigma_noise2: " << aux << std::endl
1442  << "weight: " << weight << std::endl;
1443 #endif
1444 #undef DEBUG_JM2
1445  // store weight
1446  A2D_ELEM(Mweight, i, j) = stored_weight = weight;
1447  // calculate weighted sum of (X-A)^2 for sigma_noise update
1448  local_wsum_corr += weight * diff;
1449  }
1450  else
1451  {
1452  // t-student distribution
1453  // pdf = (1 + diff2/sigma2*df)^df2
1454  // Correcting for mindiff:
1455  // pdfc = (1 + diff2/sigma2*df)^df2 / (1 + mindiff/sigma2*df)^df2
1456  // = ( (1 + diff2/sigma2*df)/(1 + mindiff/sigma2*df) )^df2
1457  // = ( (sigma2*df + diff2) / (sigma2*df + mindiff) )^df2
1458  // Extra factor two because we saved 0.5*diff2!!
1459  aux = (dfsigma2 + 2. * diff)
1460  / (dfsigma2 + 2. * trymindiff);
1461  weight = pow(aux, df2) * pdf;
1462  // Calculate extra weight acc. to Eq (10) Wang et al.
1463  // Patt. Recognition Lett. 25, 701-710 (2004)
1464  weight2 = (df + ddim2) / (df + (2. * diff / sigma_noise2));
1465  // Store probability weights
1466  stored_weight = weight * weight2;
1467  A2D_ELEM(Mweight, i, j) = stored_weight;
1468  // calculate weighted sum of (X-A)^2 for sigma_noise update
1469  local_wsum_corr += stored_weight * diff;
1470  refw2[output_refno] += stored_weight;
1471  }
1472 
1473  local_mindiff = XMIPP_MIN(local_mindiff, diff);
1474 
1475  // Accumulate sum weights for this (my) matrix
1476  my_sumweight += weight;
1477  my_sumstoredweight += stored_weight;
1478  // calculated weighted sum of offsets as well
1479  local_wsum_offset += weight * A2D_ELEM(Mr2, i, j);
1480 
1481  if (model.do_norm)
1482  {
1483  // weighted sum of Sum_j ( X_ij*A_kj )
1484  local_wsum_sc += stored_weight * (A2_plus_Xi2 - diff) / ref_scale;
1485  // weighted sum of Sum_j ( A_kj*A_kj )
1486  local_wsum_sc2 += stored_weight * A2[refno];
1487  }
1488 
1489  // keep track of optimal parameters
1490  my_maxweight = XMIPP_MAX(my_maxweight, weight);
1491 
1492  if (weight > local_maxweight)
1493  {
1494  if (model.do_student)
1495  local_maxweight2 = weight2;
1496  local_maxweight = weight;
1497  local_iopty = i;
1498  local_ioptx = j;
1499  local_iopt_psi = ipsi;
1500  local_iopt_flip = iflip;
1501  local_opt_refno = output_refno;
1502  }
1503 
1504  if (fast_mode && weight > maxw_ref)
1505  {
1506  maxw_ref = weight;
1507  iopty_ref[output_irefmir] = i;
1508  ioptx_ref[output_irefmir] = j;
1509  ioptflip_ref[output_irefmir] = iflip;
1510  }
1511 
1512  } // close for over all elements in Mweight
1513 
1514 #ifdef DEBUG_JM3
1515  if (iter == 3 && refno==0)
1516  {
1517  std::cerr << Mweight << std::endl;
1518  std::cerr << "maxweight: " << my_maxweight << " " << my_sumstoredweight << " " << my_sumweight << std::endl;
1519  exit(1);
1520  }
1521 #endif
1522  // C. only for significant settings, store weighted sums
1523  if (my_maxweight > SIGNIFICANT_WEIGHT_LOW * maxweight)
1524  {
1525  sumw_refpsi[output_refno * nr_psi + ipsi] += my_sumstoredweight;
1526 
1527  if (iflip < nr_nomirror_flips)
1528  refw[output_refno] += my_sumweight;
1529  else
1530  refw_mirror[output_refno] += my_sumweight;
1531 
1532  // Back from smaller Mweight to original size of Maux
1533  Maux.initZeros();
1534 
1535  FOR_ALL_ELEMENTS_IN_ARRAY2D(Mweight)
1536  {
1537  A2D_ELEM(Maux, i, j) = A2D_ELEM(Mweight, i, j);
1538  }
1539 
1540  // Use forward FFT in convolution theorem again
1541  // Takes the input from Maux and leaves it in Faux
1542  local_transformer.FourierTransform();
1543 
1544  MultidimArray< std::complex<double> > &mysumimgs_ref = mysumimgs[output_refnoipsi];
1545  MultidimArray< std::complex<double> > &Fimg_flip_ref = Fimg_flip[iflip];
1546 
1547  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Faux)
1548  {
1549  DIRECT_MULTIDIM_ELEM(mysumimgs_ref, n) +=
1550  conj(DIRECT_MULTIDIM_ELEM(Faux,n)) * DIRECT_MULTIDIM_ELEM(Fimg_flip_ref,n);
1551  }
1552  }
1553  } // close if Msignificant
1554  } // close for ipsi
1555  } // close for iflip
1556  } // close if pdf_directions
1557 
1558  pthread_mutex_lock(&update_mutex);
1559  //Update maxweight
1560  if (local_maxweight > maxweight)
1561  {
1562  maxweight = local_maxweight;
1563 
1564  if (model.do_student)
1565  maxweight2 = local_maxweight2;
1566  else
1567  maxweight2=0.;
1568  iopty = local_iopty;
1569  ioptx = local_ioptx;
1570  iopt_psi = local_iopt_psi;
1571  iopt_flip = local_iopt_flip;
1572  opt_refno = local_opt_refno;
1573  }
1574 
1575  //Update sums
1576  sum_refw += refw[output_refno] + refw_mirror[output_refno];
1577  wsum_offset += local_wsum_offset;
1578  wsum_corr += local_wsum_corr;
1579 
1580 
1581 
1582  mindiff = XMIPP_MIN(mindiff, local_mindiff);
1583 
1584  if (model.do_norm)
1585  {
1586  wsum_sc += local_wsum_sc;
1587  wsum_sc2 += local_wsum_sc2;
1588  }
1589  pthread_mutex_unlock(&update_mutex);
1590 
1591  //Ask for next job
1592  } // close while refno
1593 
1594 #define DP(x) << std::setw(15) << x
1595  //#define DEBUG_JM
1596 #ifdef DEBUG_JM
1597  if (iter > 1)
1598  {
1599  std::cerr << "DEBUG_JM: ====== iter: " << iter << "======= block: " << current_block << std::endl;
1600 
1601  std::cerr << formatString("====> img: %lu\n", current_image);
1602  std::cerr << "Xi2: " << Xi2 << std::endl;
1603  std::cerr << "sum_refw: " << sum_refw << std::endl;
1604 
1605  std::cerr DP("A2") DP("refw") DP("refw_mirror") << std::endl;
1606  for (int refno = 0; refno < model.n_ref; ++refno)
1607  std::cerr DP(A2[refno]) DP(refw[refno]) DP(refw_mirror[refno]) << std::endl;
1608  // std::cerr << "local_wsum_corr: " << local_wsum_corr << std::endl;
1609  // std::cerr << "wsum_corr: " << wsum_corr << std::endl;
1610  if (iter > 2)
1611  {
1612  std::cerr << "DEBUG_JM: ..............EXITING....................: " << std::endl;
1613 
1614  exit(1);
1615  }
1616  }
1617 #endif
1618 #undef DEBUG_JM
1619 
1620 
1621 
1622 }//close function doThreadExpectationSingleImage
1623 
1624 void ProgML2D::doThreadESIUpdateRefno()
1625 {
1626 
1627  double scale_dim2_sumw = (opt_scale * ddim2) / sum_refw;
1628  int num_refs = model.n_ref * factor_nref;
1629 
1630  FOR_ALL_THREAD_REFNO()
1631  {
1632  int output_refno = mygroup * model.n_ref + refno;
1633 
1634  if (fast_mode)
1635  {
1636  for (int group = 0; group < factor_nref; group++)
1637  {
1638  int group_refno = group * model.n_ref + refno;
1639 
1640  // Update optimal offsets for refno (and its mirror)
1641  allref_offsets[2 * group_refno] = -(double) ioptx_ref[output_refno]
1642  * MAT_ELEM(F[ioptflip_ref[output_refno]], 0, 0)
1643  - (double) iopty_ref[output_refno]
1644  * MAT_ELEM(F[ioptflip_ref[output_refno]], 0, 1);
1645  allref_offsets[2 * group_refno + 1] = -(double) ioptx_ref[output_refno]
1646  * MAT_ELEM(F[ioptflip_ref[output_refno]], 1, 0)
1647  - (double) iopty_ref[output_refno]
1648  * MAT_ELEM(F[ioptflip_ref[output_refno]], 1, 1);
1649  if (do_mirror)
1650  {
1651  allref_offsets[2 * (num_refs + group_refno)]
1652  = -(double) ioptx_ref[num_refs + output_refno]
1653  * MAT_ELEM(F[ioptflip_ref[num_refs + output_refno]], 0, 0)
1654  - (double) iopty_ref[num_refs + output_refno]
1655  * MAT_ELEM(F[ioptflip_ref[num_refs + output_refno]], 0, 1);
1656  allref_offsets[2 * (num_refs + group_refno) + 1]
1657  = -(double) ioptx_ref[num_refs + output_refno]
1658  * MAT_ELEM(F[ioptflip_ref[num_refs + output_refno]], 1, 0)
1659  - (double) iopty_ref[num_refs + output_refno]
1660  * MAT_ELEM(F[ioptflip_ref[num_refs + output_refno]], 1, 1);
1661  }
1662  }
1663  }
1664 
1665  if (!limit_rot || pdf_directions[refno] > 0.)
1666  {
1667  sumw[output_refno] += (refw[output_refno] + refw_mirror[output_refno]) / sum_refw;
1668  sumw2[output_refno] += refw2[output_refno] / sum_refw;
1669  sumw_mirror[output_refno] += refw_mirror[output_refno] / sum_refw;
1670 
1671  if (model.do_student)
1672  {
1673  sumwsc[output_refno] += refw2[output_refno] * opt_scale / sum_refw;
1674  sumwsc2[output_refno] += refw2[output_refno] * (opt_scale * opt_scale)
1675  / sum_refw;
1676  }
1677  else
1678  {
1679  sumwsc[output_refno] += (refw[output_refno] + refw_mirror[output_refno])
1680  * opt_scale / sum_refw;
1681  sumwsc2[output_refno] += (refw[output_refno] + refw_mirror[output_refno])
1682  * (opt_scale * opt_scale) / sum_refw;
1683  }
1684 
1685  std::complex<double> cscale_dim2_sumw=scale_dim2_sumw;
1686  for (size_t ipsi = 0; ipsi < nr_psi; ipsi++)
1687  {
1688  int refnoipsi = output_refno * nr_psi + ipsi;
1689  // Correct weighted sum of images for new bgmean (only first element=origin in Fimg)
1690  if (model.do_norm)
1691  dAij(mysumimgs[refnoipsi],0,0) -= sumw_refpsi[refnoipsi] * (bgmean - old_bgmean) / ddim2;
1692  // Sum mysumimgs to the global weighted sum
1693  wsumimgs[refnoipsi] += (cscale_dim2_sumw * mysumimgs[refnoipsi]);
1694  }
1695  }
1696 
1697  }//close while refno
1698 
1699 }//close function doThreadESIUpdateRefno
1700 
1701 void ProgML2D::iteration()
1702 {
1703  for (current_block = 0; current_block < blocks; current_block++)
1704  {
1705  // Integrate over all images
1706  expectation();
1707  //Maximize the model
1708  maximization();
1709  }//close for blocks
1710 }
1711 
1712 
1713 void ProgML2D::expectation()
1714 {
1715  MultidimArray<std::complex<double> > Fdzero(dim, hdim + 1);
1716  int num_output_refs = factor_nref * model.n_ref;
1717 
1718  LOG(" ProgML2D::expectation BEGIN");
1719 #ifdef DEBUG
1720 
1721  std::cerr<<"entering expectation"<<std::endl;
1722 #endif
1723 
1724 #ifdef TIMING
1725 
1726  timer.tic(E_RR);
1727 #endif
1728 
1729  rotateReference();
1730 #ifdef TIMING
1731 
1732  timer.toc(E_RR);
1733  timer.tic(E_PRE);
1734 #endif
1735  // Pre-calculate pdf of all in-plane transformations
1736  calculatePdfInplane();
1737 
1738  // Initialize weighted sums
1739  LL = 0.;
1740  wsumimgs.clear();
1741 
1742  sumw.assign(num_output_refs, 0.);
1743  sumw2.assign(num_output_refs, 0.);
1744  sumwsc.assign(num_output_refs, 0.);
1745  sumwsc2.assign(num_output_refs, 0.);
1746  sumw_mirror.assign(num_output_refs, 0.);
1747 
1748  wsum_sigma_noise = 0.;
1749  wsum_sigma_offset = 0.;
1750  sumfracweight = 0.;
1751 
1752  Fdzero.initZeros();
1753  wsumimgs.assign(num_output_refs * nr_psi, Fdzero);
1754 
1755  // Local variables of old threadExpectationSingleImage
1756  Image<double> img;
1757  FileName fn_img, fn_tkrans;
1758  Matrix1D<double> opt_offsets(2);
1759  double old_phi = -999., old_theta = -999.;
1760  double opt_flip;
1761 
1762  //Some initializations
1763  opt_scale = 1., bgmean = 0.;
1764 
1765  Msignificant.resizeNoCopy(model.n_ref, nr_psi * nr_flip);
1766 
1767  static size_t img_done;
1768  if (current_block == 0) //when not iem current block is always 0
1769  {
1770  initProgress(nr_images_local);
1771  img_done = 0;
1772  }
1773 
1774  String _msg = formatString("Images: %lu, first: %lu, last: %lu", nr_images_local, myFirstImg, myLastImg);
1775  LOG(_msg.c_str());
1776  //std::cerr << "-----xmipp_current: Expectation, iter " << iter << "-------" << std::endl;
1777  //for (int imgno = 0, img_done = 0; imgno < nn; imgno++)
1778  // Loop over all images
1779  FOR_ALL_LOCAL_IMAGES()
1780  {
1781 
1782  if (IMG_BLOCK(imgno) == current_block)
1783  {
1784  current_image = imgno;
1785  //std::cerr << "\n ======>>> imgno: " << imgno << std::endl;
1786  mygroup = (factor_nref > 1) ? divide_equally_group(nr_images_global, factor_nref, imgno) : 0;
1787 
1788  MDimg.getValue(MDL_IMAGE, fn_img, img_id[imgno]);
1789  img.read(fn_img);
1790  img().setXmippOrigin();
1791  Xi2 = img().sum2();
1792  Mimg = img();
1793 
1794 
1795  //#define DEBUG_JM1
1796 #ifdef DEBUG_JM1
1797 
1798  //if (iter >= 2 && current_image == myFirstImg)
1799  printf(" ====================>>> Iter: %02d Image: %06lu: \n", iter, imgno);
1800  printf(" fn_img: %s", fn_img.c_str());
1801  //printf(" mygroup: %d", mygroup);
1802 #endif
1803 #undef DEBUG_JM1
1804 
1805  // These two parameters speed up expectationSingleImage
1806  opt_refno = imgs_optrefno[IMG_LOCAL_INDEX];
1807  trymindiff = imgs_trymindiff[IMG_LOCAL_INDEX];
1808 
1809  if (trymindiff < 0.)
1810  // 90% of Xi2 may be a good idea (factor half because 0.5*diff is calculated)
1811  trymindiff = trymindiff_factor * 0.5 * Xi2;
1812 
1813  if (model.do_norm)
1814  {
1815  bgmean = imgs_bgmean[IMG_LOCAL_INDEX];
1816  opt_scale = imgs_scale[IMG_LOCAL_INDEX];
1817  }
1818 
1819  // Get optimal offsets for all references
1820  if (fast_mode)
1821  {
1822  allref_offsets = imgs_offsets[IMG_LOCAL_INDEX];
1823  }
1824 
1825  // Read optimal orientations from memory
1826  if (limit_rot)
1827  {
1828  old_phi = imgs_oldphi[IMG_LOCAL_INDEX];
1829  old_theta = imgs_oldtheta[IMG_LOCAL_INDEX];
1830  }
1831  // For limited orientational search: preselect relevant directions
1832  preselectLimitedDirections(old_phi, old_theta);
1833 
1834  // Use a maximum-likelihood target function in real space
1835  // with complete or reduced-space translational searches (-fast)
1836  if (fast_mode)
1837  preselectFastSignificant();
1838  else
1839  Msignificant.initConstant(1);
1840 
1841  expectationSingleImage(opt_offsets);
1842 
1843  // Write optimal offsets for all references to disc
1844  if (fast_mode)
1845  {
1846  imgs_offsets[IMG_LOCAL_INDEX] = allref_offsets;
1847  }
1848 
1849  // Store mindiff for next iteration
1850  imgs_trymindiff[IMG_LOCAL_INDEX] = trymindiff;
1851 
1852  // Store opt_refno for next iteration
1853  imgs_optrefno[IMG_LOCAL_INDEX] = opt_refno;
1854 
1855  // Store optimal phi and theta in memory
1856  if (limit_rot)
1857  {
1858  imgs_oldphi[IMG_LOCAL_INDEX] = model.Iref[opt_refno % model.n_ref].rot();
1859  imgs_oldtheta[IMG_LOCAL_INDEX] = model.Iref[opt_refno % model.n_ref].tilt();
1860  }
1861 
1862  // Store optimal normalization parameters in memory
1863  if (model.do_norm)
1864  {
1865  imgs_scale[IMG_LOCAL_INDEX] = opt_scale;
1866  imgs_bgmean[IMG_LOCAL_INDEX] = bgmean;
1867  }
1868 
1869  // Output docfile
1870  opt_flip = 0.;
1871  if (-opt_psi > 360.)
1872  {
1873  opt_psi += 360.;
1874  opt_flip = 1.;
1875  }
1876 
1877  dAij(docfiledata,IMG_LOCAL_INDEX,0)
1878  = model.Iref[opt_refno % model.n_ref].rot(); // rot
1879  dAij(docfiledata,IMG_LOCAL_INDEX,1)
1880  = model.Iref[opt_refno % model.n_ref].tilt(); // tilt
1881  dAij(docfiledata,IMG_LOCAL_INDEX,2) = opt_psi + 360.; // psi
1882  dAij(docfiledata,IMG_LOCAL_INDEX,3) = opt_offsets(0); // Xoff
1883  dAij(docfiledata,IMG_LOCAL_INDEX,4) = opt_offsets(1); // Yoff
1884  dAij(docfiledata,IMG_LOCAL_INDEX,5) = (double) (opt_refno + 1); // Ref
1885  dAij(docfiledata,IMG_LOCAL_INDEX,6) = opt_flip; // Mirror
1886  dAij(docfiledata,IMG_LOCAL_INDEX,7) = fracweight; // P_max/P_tot
1887  dAij(docfiledata,IMG_LOCAL_INDEX,8) = dLL; // log-likelihood
1888  if (model.do_norm)
1889  {
1890  dAij(docfiledata,IMG_LOCAL_INDEX,9) = bgmean; // background mean
1891  dAij(docfiledata,IMG_LOCAL_INDEX,10) = opt_scale; // image scale
1892  }
1893  if (model.do_student)
1894  {
1895  dAij(docfiledata,IMG_LOCAL_INDEX,11) = maxweight2; // Robustness weight
1896  }
1897 
1898  //Report progress and increment the images done
1899  setProgress(++img_done);
1900 
1901  //#define DEBUG_JM1
1902 #ifdef DEBUG_JM1
1903  // {
1904  // //std::cerr << "---------------------- DEBUG_JM: current_image: " << current_image << std::endl;
1905  std::cerr << " LL: " << LL << std::endl;
1906  std::cerr << " wsum_sigma_noise: " << wsum_sigma_noise << std::endl;
1907  std::cerr << " wsum_sigma_offset: " << wsum_sigma_offset << std::endl;
1908  std::cerr << " sumfracweight: " << sumfracweight << std::endl;
1909  // }
1910 #endif
1911 #undef DEBUG_JM1
1912 
1913  }//close if current_block, also close of for all images
1914  }//close for all images
1915 
1916  if (current_block == (blocks - 1))
1917  endProgress();
1918 
1919  //Changes temporally the model n_ref for the
1920  //refno loop, but not yet n_ref because in iem
1921  //isn't yet the end of iteration
1922  model.n_ref *= factor_nref;
1923  // Rotate back and calculate weighted sums
1924  reverseRotateReference();
1925  //Restore back the model.n_ref
1926  model.n_ref /= factor_nref;
1927  LOG(" ProgML2D::expectation END");
1928 
1929 }//function expectation
1930 
1931 
1932 // Update all model parameters
1933 void ProgML2D::maximizeModel(ModelML2D &local_model)
1934 {
1935 
1936 #ifdef DEBUG
1937  std::cerr<<"entering maximization"<<std::endl;
1938 #endif
1939 
1940  Matrix1D<double> rmean_sigma2, rmean_signal2;
1941  Matrix1D<int> center(2), radial_count;
1942  MultidimArray<std::complex<double> > Faux, Faux2;
1943  MultidimArray<double> Maux;
1944  FileName fn_tmp;
1945 
1946  // After iteration 0, factor_nref will ALWAYS be one
1947  if (factor_nref > 1)
1948  {
1949  int old_refno = local_model.n_ref;
1950  local_model.setNRef(local_model.n_ref * factor_nref);
1951  // Now also expand the Iref vector to contains factor_nref times more images
1952  // Make sure that headers are the same by using = assignments
1953  for (int group = 1; group < factor_nref; group++)
1954  {
1955  for (int refno = 0; refno < old_refno; refno++)
1956  {
1957 
1958  local_model.Iref[group * old_refno + refno] = local_model.Iref[refno];
1959  }
1960  }
1961  }
1962 
1963  // Update the reference images
1964  local_model.sumw_allrefs = 0.;
1965  local_model.sumw_allrefs2 = 0.;
1966 
1967 
1968  //#define ASSIGN(var) if (var > SIGNIFICANT_WEIGHT_LOW) local_model.var = var
1969 #define ASSIGN(var) local_model.var = var
1970 
1971  ASSIGN(dim);
1972  ASSIGN(sumfracweight);
1973  local_model.LL = LL;
1974  // Update sigma of the origin offsets
1975  if (!fix_sigma_offset)
1976  ASSIGN(wsum_sigma_offset);
1977  // Update the noise parameters
1978  if (!fix_sigma_noise)
1979  ASSIGN(wsum_sigma_noise);
1980 
1981 
1982  for (int refno = 0; refno < local_model.n_ref; refno++)
1983  {
1984  double weight = sumw[refno];
1985  if (weight > 0.)
1986  {
1987  if (local_model.do_student)
1988  {
1989  weight = sumw2[refno];
1990  local_model.sumw_allrefs2 += sumw2[refno];
1991  }
1992 
1993  double avg, stddev, min, max;
1994  wsum_Mref[refno].computeStats(avg, stddev, min, max);
1995  //std::cerr << formatString("DEBUG_JM: avg %f, stddev %f, min %f, max %f", avg, stddev, min, max) <<std::endl;
1996  local_model.WsumMref[refno]() = wsum_Mref[refno];
1997  //local_model.Iref[refno]() *= 1/weight;//fixme: sumwsc2[refno];
1998  local_model.sumw_allrefs += sumw[refno];
1999  local_model.WsumMref[refno].setWeight(weight);
2000  }
2001  else
2002  {
2003  local_model.WsumMref[refno].setWeight(0.);
2004  local_model.WsumMref[refno]().initZeros(dim, dim);
2005  local_model.WsumMref[refno]().setXmippOrigin();
2006  }
2007 
2008  // Adjust average scale (nr_classes will be smaller than model.n_ref for the 3D case!)
2009  if (local_model.do_norm)
2010  ASSIGN(sumwsc[refno]);
2011 
2012  if (!fix_fractions)
2013  ASSIGN(sumw_mirror[refno]);
2014  }
2015 
2016 #ifdef DEBUG
2017 
2018  std::cerr<<"leaving maximization"<<std::endl;
2019 
2020 #endif
2021 }//close function maximization
2022 
2023 void ProgML2D::maximization()
2024 {
2025  LOG(" ProgML2D::maximization BEGIN");
2026  if (blocks == 1) //ie not IEM, normal maximization
2027  {
2028  maximizeModel(model);
2029  model.update();
2030  // After iteration 0, factor_nref will ALWAYS be one
2031  factor_nref = 1;
2032  }
2033  else //do IEM
2034  {
2035  bool special_first = (!do_restart && iter == istart);
2036  ModelML2D block_model(model.n_ref);
2037 
2038  LOG(" ProgML2D::maximization: readModel and substractModel");
2039  if (!special_first)
2040  {
2041  readModel(block_model, current_block);
2042  model.substractModel(block_model);
2043  }
2044  LOG(" ProgML2D::maximization: maximizeModel");
2045  maximizeModel(block_model);
2046  LOG(" ProgML2D::maximization: maximizeModel");
2047  writeOutputFiles(block_model, OUT_BLOCK);
2048  LOG(" ProgML2D::maximization: addModels");
2049  if (!special_first)
2050  {
2051  model.addModel(block_model);
2052  model.update();
2053  }
2054  else if (current_block == blocks - 1) //last block
2055  {
2056  for (current_block = 0; current_block < blocks; current_block++)
2057  {
2058  readModel(block_model, current_block);
2059 
2060  if (current_block == 0)
2061  model = block_model;
2062  else
2063  model.addModel(block_model);
2064  }
2065  model.update();
2066  //printModel("GLOBAL model: after addition and update: ", model);
2067  // After iteration 0, factor_nref will ALWAYS be one
2068  factor_nref = 1;
2069  //restore the value of current block
2070  --current_block;
2071  }
2072  }
2073 
2074  if (model.do_norm)
2075  correctScaleAverage();
2076 
2077  LOG(" ProgML2D::maximization END");
2078 
2079  // static int times = 0;
2080  // if (times > 1)
2081  // exit(1);
2082  // ++times;
2083 }//close function maximizationBlocks
2084 
2085 void ProgML2D::correctScaleAverage()
2086 {
2087  //FIXME: This function needs re-implementation for the do_norm parameter
2088  //now the ref3d comes in metadata, refs_per_class can be avoided
2089  // int iclass, nr_classes = ROUND(model.n_ref / refs_per_class);
2090  // std::vector<double> wsum_scale(nr_classes), sumw_scale(nr_classes);
2091  // ldiv_t temp;
2092  // average_scale = 0.;
2093  //
2094  // for (int refno = 0; refno < model.n_ref; refno++)
2095  // {
2096  // average_scale += model.getSumwsc(refno);
2097  // temp = ldiv(refno, refs_per_class);
2098  // iclass = ROUND(temp.quot);
2099  // wsum_scale[iclass] += model.getSumwsc(refno);
2100  // sumw_scale[iclass] += model.getSumw(refno);
2101  // }
2102  // for (int refno = 0; refno < model.n_ref; refno++)
2103  // {
2104  // temp = ldiv(refno, refs_per_class);
2105  // iclass = ROUND(temp.quot);
2106  // if (sumw_scale[iclass] > 0.)
2107  // {
2108  // model.scale[refno] = wsum_scale[iclass] / sumw_scale[iclass];
2109  // model.WsumMref[refno]() *= model.scale[refno];
2110  // }
2111  // else
2112  // {
2113  // model.scale[refno] = 1.;
2114  // }
2115  // }
2116  // average_scale /= model.sumw_allrefs;
2117 }//close function correctScaleAverage
2118 
2120 void ProgML2D::addPartialDocfileData(const MultidimArray<double> &data,
2121  size_t first, size_t last)
2122 {
2123 #ifdef DEBUG
2124  std::cerr << "Entering addPartialDocfileData" <<std::endl;
2125 #endif
2126 
2127  for (size_t imgno = first; imgno <= last; imgno++)
2128  {
2129  size_t index = imgno - first;
2130  size_t id = img_id[imgno];
2131  if (do_ML3D)
2132  {
2133  MDimg.setValue(MDL_ANGLE_ROT, dAij(data, index, 0), id);
2134  MDimg.setValue(MDL_ANGLE_TILT, dAij(data, index, 1), id);
2135  }
2136  //Here we change the sign of the angle becase in the code
2137  //is represent the rotation of the reference and we want
2138  //to store the rotation of the image
2139  double psi = -dAij(data, index, 2);
2140  MDimg.setValue(MDL_ANGLE_PSI, psi, id);
2141  MDimg.setValue(MDL_SHIFT_X, dAij(data, index, 3), id);
2142  MDimg.setValue(MDL_SHIFT_Y, dAij(data, index, 4), id);
2143  MDimg.setValue(MDL_REF, (int)round(dAij(data, index, 5)), id);
2144  if (do_mirror)
2145  {
2146  MDimg.setValue(MDL_FLIP, dAij(data, index, 6) != 0., id);
2147  }
2148  MDimg.setValue(MDL_PMAX, dAij(data, index, 7), id);
2149  MDimg.setValue(MDL_LL, dAij(data, index, 8), id);
2150  if (model.do_norm)
2151  {
2152  MDimg.setValue(MDL_BGMEAN, dAij(data, index, 9), id);
2153  MDimg.setValue(MDL_INTSCALE, dAij(data, index, 10), id);
2154  }
2155  if (model.do_student)
2156  {
2157  MDimg.setValue(MDL_WROBUST, dAij(data, index, 11), id);
2158  }
2159  }
2160 
2161 #ifdef DEBUG
2162  std::cerr << "Leaving addPartialDocfileData" <<std::endl;
2163 #endif
2164 }//close function addDocfileData
2165 
2166 //Some macros
2167 #define ITER_PREFIX "iter"//formatString("iter%06d", iter)
2168 #define FINAL_PREFIX "final"
2169 #define IS_FINAL outputType == OUT_FINAL
2170 
2171 void ProgML2D::writeOutputFiles(const ModelML2D &model, OutputType outputType)
2172 {
2173  FileName fn_tmp, fn_prefix, fn_base;
2174  Image<double> Itmp;
2175  MetaDataVec MDo;
2176  bool write_img_xmd = true, write_refs_log = true, write_conv = !do_ML3D;
2177  bool write_norm = model.do_norm;
2178 
2179  if (iter == 0)
2180  write_conv = false;
2181 
2182  //By default write down the parameters
2183  std::vector<Image<double> > const * ptrImages = &(model.Iref);
2184  std::vector<double> const * ptrMirror = &(model.mirror_fraction);
2185  double avePmax = model.avePmax;
2186  double sigma_noise = model.sigma_noise;
2187  double sigma_offset = model.sigma_offset;
2188 
2189  switch (outputType)
2190  {
2191  case OUT_BLOCK:
2192  write_img_xmd = false;
2193  write_conv = false;
2194  // Sjors 18may2010: why not write scales in blocks? We need this now, don't we?
2195  //write_norm = false;
2196  fn_prefix = formatString("block%03d", current_block + 1);
2197  // For blocks we will write out the summations instead of parameters
2198  ptrImages = &(model.WsumMref);
2199  ptrMirror = &(model.sumw_mirror);
2200  avePmax = model.sumfracweight;
2201  sigma_noise = model.wsum_sigma_noise;
2202  sigma_offset = model.wsum_sigma_offset;
2203  break;
2204  case OUT_ITER:
2205  //fn_base = getBaseName("_it", iter);
2206  fn_prefix = ITER_PREFIX;
2207  break;
2208  case OUT_FINAL:
2209  //fn_base = fn_root;
2210  fn_prefix = FINAL_PREFIX;
2211  break;
2212  case OUT_IMGS:
2213  //std::cerr << "OUT_IMGS" <<std::endl;
2214  write_refs_log = false;
2215  //fn_base = getBaseName("_it", iter);
2216  fn_prefix = ITER_PREFIX;
2217  break;
2218  case OUT_REFS:
2219  //std::cerr << "OUT_REFS" <<std::endl;
2220  LOG("OUT_REFS");
2221  write_img_xmd = false;
2222  write_conv = false;
2223  //fn_base = getBaseName("_it", iter);
2224  fn_prefix = ITER_PREFIX;
2225  break;
2226  }
2227 
2228  fn_base = (fn_prefix == ITER_PREFIX) ? //All intermediate iteration files should go to "extra" folder
2229  getIterExtraPath(fn_root, iter) : fn_root;
2230 
2231  if (write_img_xmd)
2232  {
2233  //static WriteModeMetaData mode = MD_OVERWRITE;
2234  //Write image metadata, for each iteration a new block will be written
2235  fn_tmp = FN_IMAGES_MD(fn_base);
2236  MDimg.write(fn_tmp);
2237  // if (fn_prefix == ITER_PREFIX)
2238  // fn_tmp = formatString("iter%06d@%s_iter_images.xmd", iter, rootStr);
2239  // else
2240  // fn_tmp = formatString("%s_%s_images.xmd", rootStr, prefixStr);
2241  //MDimg.write(fn_tmp, mode);
2242  //mode = MD_APPEND;
2243  }
2244 
2245 
2246  if (write_refs_log)
2247  {
2248  // Write out current reference images and fill sel & log-file
2249  // Re-use the MDref metadata that was read in produceSideInfo2
2250  // This way. MDL_ANGLE_ROT, MDL_ANGLE_TILT, MDL_REF etc are treated ok for do_ML3D
2251  int refno = 0;
2252  int select_img;
2253  double weight;
2254  size_t count;
2255 
2256  fn_tmp = FN_CLASSES_STK(fn_base);
2257  //
2258  // if (fn_prefix == ITER_PREFIX)
2259  // fn_tmp = formatString("%s_iter%06d_refs.stk", rootStr, iter);
2260  // else
2261  // fn_tmp = formatString("%s_%s_refs.stk", rootStr, prefixStr);
2262 
2263  for (size_t objId : MDref.ids())
2264  {
2265  select_img = refno + 1;
2266  //Itmp = model.Iref[refno];
2267  //std::cerr << "DEBUG_JM: refno: " << refno << std::endl;
2268  Itmp = (*ptrImages)[refno];
2269  //std::cerr << "DEBUG_JM: fn_tmp: " << fn_tmp << std::endl;
2270  //std::cerr << "DEBUG_JM: select_img: " << select_img << std::endl;
2271  Itmp.write(fn_tmp, select_img, true, WRITE_REPLACE);
2272  //MDref.setValue(MDL_ITER, iter, objId);//write out iteration number
2273  MDref.setValue(MDL_REF, select_img, objId); //Also write reference number
2274  MDref.setValue(MDL_IMAGE, formatString("%06d@%s", select_img, fn_tmp.c_str()), objId);
2275  MDref.setValue(MDL_ENABLED, 1, objId);
2276  weight = model.WsumMref[refno].weight();
2277  MDref.setValue(MDL_WEIGHT, weight, objId);
2278  count = ROUND(weight);
2279  MDref.setValue(MDL_CLASS_COUNT, count, objId);
2280  if (do_mirror)
2281  MDref.setValue(MDL_MIRRORFRAC, (*ptrMirror)[refno], objId);
2282  if (write_conv)
2283  MDref.setValue(MDL_SIGNALCHANGE, conv[refno]*1000, objId);
2284  if (write_norm)
2285  MDref.setValue(MDL_INTSCALE, model.scale[refno], objId);
2286  ++refno;
2287  }
2288  //fn_tmp.copyFile(formatString("%s_output_block%d.stk", fn_tmp.c_str(), current_block));
2289 
2290  // fn_tmp = formatString("%s_%s", rootStr, prefixStr);
2291  // FileName fn_ref = fn_tmp + "_refs.xmd";
2292  //
2293  // if (!fn_prefix.contains("block"))
2294  // fn_ref = formatString("iter%06d@%s", iter, fn_ref.c_str());
2295  // MDref.write(fn_ref, mode);
2296  fn_tmp = FN_CLASSES_MD(fn_base);
2297  MDref.write(fn_tmp);
2298 
2299  if (outputType == OUT_REFS)
2300  outRefsMd = fn_tmp;
2301 
2302  // Write out log-file
2303  MetaDataVec mdLog;
2304  size_t objId = mdLog.addObject();
2305  mdLog.setValue(MDL_ITER, iter, objId);
2306  mdLog.setValue(MDL_LL, model.LL, objId);
2307  mdLog.setValue(MDL_PMAX, avePmax, objId);
2308  mdLog.setValue(MDL_SIGMANOISE, sigma_noise, objId);
2309  mdLog.setValue(MDL_SIGMAOFFSET, sigma_offset, objId);
2310  mdLog.setValue(MDL_RANDOMSEED, seed, objId);
2311  if (write_norm)
2312  mdLog.setValue(MDL_INTSCALE, average_scale, objId);
2313 
2314  mdLog.write(FN_LOGS_MD(fn_base), MD_APPEND);
2315 
2316  if (write_img_xmd)
2317  {
2318  MetaDataVec mdImgs;
2319  size_t n = MDref.size();
2320  for (size_t ref = 1; ref <= n; ++ref)
2321  {
2322  mdImgs.importObjects(MDimg, MDValueEQ(MDL_REF, (int)ref));
2323  mdImgs.write(FN_CLASS_IMAGES_MD(fn_base, ref), MD_APPEND);
2324  }
2325  }
2326 
2327  // if (fn_prefix.contains("block") || mode == MD_OVERWRITE)
2328  // mdLog.write(fn_tmp);
2329  // else
2330  // mdLog.append(fn_tmp);
2331  // mode = MD_APPEND;
2332  }
2333 
2334 }//close function writeModel
2335 
2336 void ProgML2D::readModel(ModelML2D &model, int block)
2337 {
2338 
2339  // First read general model parameters from _log.xmd
2340  FileName fn_base = formatString("%s_block%03d", fn_root.c_str(), (block + 1));
2341  MetaDataVec MDi;
2342  MDi.read(fn_base + "_logs.xmd");
2343  model.dim = dim;
2344 
2345  {
2346  size_t id = MDi.firstRowId();
2347  MDi.getValue(MDL_LL, model.LL, id);
2348  MDi.getValue(MDL_PMAX, model.sumfracweight, id);
2349  MDi.getValue(MDL_SIGMANOISE, model.wsum_sigma_noise, id);
2350  MDi.getValue(MDL_SIGMAOFFSET, model.wsum_sigma_offset, id);
2351  }
2352 
2353  //Then read reference model parameters from _ref.xmd
2354  FileName fn_img;
2355  Image<double> img;
2356  MDi.clear();
2357  MDi.read(fn_base + "_refs.xmd");
2358  int refno = 0;
2359  double weight = 0;
2360 
2361  model.sumw_allrefs = 0.;
2362  for (size_t objId : MDi.ids())
2363  {
2364  MDi.getValue(MDL_IMAGE, fn_img, objId);
2365  img.read(fn_img);
2366  img().setXmippOrigin();
2367  model.WsumMref[refno] = img;
2368  MDi.getValue(MDL_WEIGHT, weight, objId);
2369  model.WsumMref[refno].setWeight(weight);
2370  model.sumw_allrefs += weight;
2371  model.sumw_mirror[refno] = 0.;
2372  if (do_mirror)
2373  MDi.getValue(MDL_MIRRORFRAC, model.sumw_mirror[refno], objId);
2374  refno++;
2375  }
2376 }//close function readModel
2377 
2378 FileName ProgML2D::getBaseName(String suffix, int number)
2379 {
2380  FileName fn_base = fn_root + suffix;
2381  if (number >= 0)
2382  fn_base.compose(fn_base, number, "");
2383  return fn_base;
2384 }
2385 
ProgML2D::getThreadRefnoJob
int getThreadRefnoJob(int &refno)
Assign refno jobs to threads.
Definition: ml_align2d.cpp:1026
ML2DBaseProgram::nr_psi
size_t nr_psi
Definition: ml2d.h:154
selfApplyGeometry
void selfApplyGeometry(int Splinedegree, MultidimArray< std::complex< double > > &V1, const Matrix2D< double > &A, bool inv, bool wrap, std::complex< double > outside)
Definition: transformations.cpp:704
barrier_init
int barrier_init(barrier_t *barrier, int needed)
Definition: xmipp_threads.cpp:366
ProgML2D::threadTask
int threadTask
Definition: ml_align2d.h:67
ProgML2D::printModel
virtual void printModel(const String &msg, const ModelML2D &model)
Definition: ml_align2d.cpp:281
MDL_LL
contribution of an image to log-likelihood value
Definition: metadata_label.h:269
ML2DBaseProgram::img_id
std::vector< size_t > img_id
Definition: ml2d.h:232
ML2DBaseProgram::ddim2
double ddim2
Definition: ml2d.h:152
ML2DBaseProgram::MDimg
MetaDataDb MDimg
Definition: ml2d.h:172
ML2DBaseProgram::refs_per_class
int refs_per_class
Definition: ml2d.h:217
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
OUT_REFS
Definition: ml2d.h:76
min
void min(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:118
ProgML2D::preselectLimitedDirections
void preselectLimitedDirections(double &phi, double &theta)
Definition: ml_align2d.cpp:634
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
OUT_ITER
Definition: ml2d.h:76
ProgML2D::addPartialDocfileData
virtual void addPartialDocfileData(const MultidimArray< double > &data, size_t first, size_t last)
Add docfiledata to docfile.
Definition: ml_align2d.cpp:2120
MultidimArray::resize
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
Definition: multidim_array.cpp:983
XMIPP_MAX
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
FINAL_PREFIX
#define FINAL_PREFIX
Definition: ml_align2d.cpp:2168
ProgML2D::barrier2
barrier_t barrier2
Definition: ml_align2d.h:68
ProgML2D::opt_scale
double opt_scale
Definition: ml_align2d.h:85
ERR_NOT_IMPLEMENTED
Case or algorithm not implemented yet.
Definition: xmipp_error.h:177
ML2DBaseProgram::nr_flip
size_t nr_flip
Definition: ml2d.h:156
metadata_sql.h
alglib::conj
alglib::complex conj(const alglib::complex &z)
Definition: ap.cpp:4988
ProgML2D::dLL
double dLL
Definition: ml_align2d.h:86
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
ProgML2D::old_bgmean
double old_bgmean
Definition: ml_align2d.h:93
ML2DBaseProgram::zero_offsets
bool zero_offsets
Definition: ml2d.h:186
FourierTransformer::inverseFourierTransform
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
MetaDataVec::read
void read(const FileName &inFile, const std::vector< MDLabel > *desiredLabels=nullptr, bool decomposeStack=true) override
Definition: metadata_vec.cpp:120
ModelML2D::LL
double LL
Definition: ml2d.h:99
ProgML2D::wsum_sc
double wsum_sc
Definition: ml_align2d.h:93
SIGNIFICANT_WEIGHT_LOW
#define SIGNIFICANT_WEIGHT_LOW
Definition: ml_tomo.h:46
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
FOR_ALL_LOCAL_IMAGES
#define FOR_ALL_LOCAL_IMAGES()
Definition: ml2d.h:64
apply_binary_mask
void apply_binary_mask(const MultidimArray< int > &mask, const MultidimArray< T > &m_in, MultidimArray< T > &m_out, T subs_val=(T) 0)
Definition: mask.h:857
FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
Definition: multidim_array_base.h:509
ProgML2D::Mimg
MultidimArray< double > Mimg
Definition: ml_align2d.h:78
ProgML2D::sumw_mirror
std::vector< double > sumw_mirror
Definition: ml_align2d.h:75
Euler_direction
void Euler_direction(double alpha, double beta, double gamma, Matrix1D< double > &v)
Definition: geometry.cpp:721
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
dAij
#define dAij(M, i, j)
Definition: multidim_array_base.h:345
MD_APPEND
Definition: metadata_writemode.h:34
MetaDataDb::findObjects
void findObjects(std::vector< size_t > &objectsOut, const MDQuery &query) const override
Definition: metadata_db.cpp:709
ProgML2D::sigdim
size_t sigdim
Definition: ml_align2d.h:95
ProgML2D::iopty_ref
std::vector< int > iopty_ref
Definition: ml_align2d.h:97
ProgML2D::pfs_mindiff
double pfs_mindiff
Definition: ml_align2d.h:99
ProgML2D::wsum_sigma_offset
double wsum_sigma_offset
Definition: ml_align2d.h:74
ML2DBaseProgram::Mr2
MultidimArray< double > Mr2
Definition: ml2d.h:178
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
BinaryCircularMask
void BinaryCircularMask(MultidimArray< int > &mask, double radius, int mode, double x0, double y0, double z0)
Definition: mask.cpp:193
MetaDataDb::getValue
bool getValue(MDObject &mdValueOut, size_t id) const override
Definition: metadata_db.cpp:109
ML2DBaseProgram::current_block
size_t current_block
Definition: ml2d.h:229
ML2DBaseProgram::dfsigma2
double dfsigma2
Definition: ml2d.h:204
MDL_SIGNALCHANGE
Signal change for an image.
Definition: metadata_label.h:447
ML2DBaseProgram::model
ModelML2D model
Definition: ml2d.h:224
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
ProgML2D::createThreads
void createThreads()
Create working threads.
Definition: ml_align2d.cpp:929
y
static double * y
Definition: numerical_recipes.cpp:8487
ProgML2D::ioptx_ref
std::vector< int > ioptx_ref
Definition: ml_align2d.h:97
ML2DBaseProgram::fix_fractions
bool fix_fractions
Definition: ml2d.h:139
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
ProgML2D::sum_refw
double sum_refw
Definition: ml_align2d.h:92
IMG_BLOCK
#define IMG_BLOCK(imgno)
Definition: ml2d.h:69
ProgML2D::wsum_corr
double wsum_corr
Definition: ml_align2d.h:92
ModelML2D::dim
int dim
Definition: ml2d.h:102
ProgML2D::iopt_flip
int iopt_flip
Definition: ml_align2d.h:84
MDL_SHIFT_X
Shift for the image in the X axis (double)
Definition: metadata_label.h:431
ML2DBaseProgram::average_scale
double average_scale
Definition: ml2d.h:213
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
XmippProgram::initProgress
void initProgress(size_t total, size_t stepBin=60)
Definition: xmipp_program.cpp:251
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
ProgML2D::readParams
void readParams()
Read arguments from command line.
Definition: ml_align2d.cpp:75
ProgML2D::fref
std::vector< MultidimArray< std::complex< double > > > fref
Definition: ml_align2d.h:82
getImageSize
void getImageSize(const MetaData &md, size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim, MDLabel image_label)
Definition: metadata_extension.cpp:273
ModelML2D::do_student_sigma_trick
bool do_student_sigma_trick
Definition: ml2d.h:104
centerFFT2
void centerFFT2(MultidimArray< double > &v)
Definition: xmipp_fft.cpp:276
ML2DBaseProgram::myLastImg
size_t myLastImg
Definition: ml2d.h:168
FN_LOGS_MD
#define FN_LOGS_MD(base)
Definition: ml2d.h:56
ProgML2D::maximization
virtual void maximization()
Update all model parameters, adapted for IEM blocks use.
Definition: ml_align2d.cpp:2023
ModelML2D::mirror_fraction
std::vector< double > mirror_fraction
Definition: ml2d.h:93
FileName::compose
void compose(const String &str, const size_t no, const String &ext="")
Definition: xmipp_filename.cpp:52
ProgML2D::doThreadExpectationSingleImageRefno
void doThreadExpectationSingleImageRefno()
Thread code to parallelize refno loop in expectationSingleImage.
Definition: ml_align2d.cpp:1287
MDL_ANGLE_PSI
Special label to be used when gathering MDs in MpiMetadataPrograms.
Definition: metadata_label.h:51
MultidimArray::initConstant
void initConstant(T val)
Definition: multidim_array.h:2705
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
ModelML2D::substractModel
void substractModel(const ModelML2D &model)
Definition: ml2d.cpp:409
ML2DBaseProgram::fn_img
FileName fn_img
Definition: ml2d.h:132
ModelML2D::n_ref
int n_ref
Definition: ml2d.h:83
ProgML2D::omask
MultidimArray< int > omask
Definition: ml_align2d.h:65
WRITE_REPLACE
Definition: xmipp_image_base.h:78
ML2DBaseProgram::psi_step
double psi_step
Definition: ml2d.h:158
Matrix1D::selfNormalize
void selfNormalize()
Definition: matrix1d.cpp:723
OutputType
OutputType
Definition: ml2d.h:76
MDL_WROBUST
Weight of t-student distribution in robust Maximum likelihood.
Definition: metadata_label.h:494
ProgML2D::destroyThreads
void destroyThreads()
Exit threads and free memory.
Definition: ml_align2d.cpp:958
MDL_BGMEAN
<Bfactor of a map, or even a local bfactor
Definition: metadata_label.h:76
ThreadTask
ThreadTask
Definition: ml2d.h:74
ModelML2D::setNRef
void setNRef(int n_ref)
Definition: ml2d.cpp:342
ProgML2D::produceSideInfo2
virtual void produceSideInfo2()
Try to merge produceSideInfo1 and 2.
Definition: ml_align2d.cpp:357
ProgML2D::threads_d
structThreadTasks * threads_d
Definition: ml_align2d.h:70
ProgML2D::readModel
virtual void readModel(ModelML2D &model, int block)
Read model from file.
Definition: ml_align2d.cpp:2336
ProgML2D::pfs_count
int pfs_count
Definition: ml_align2d.h:102
ML2DBaseProgram::hdim
size_t hdim
Definition: ml2d.h:151
MetaData::ids
virtual IdIteratorProxy< false > ids()
Definition: metadata_base.h:754
ML2DBaseProgram::fn_root
FileName fn_root
Definition: ml2d.h:132
MAX_WEIGHT
#define MAX_WEIGHT
Definition: ml_align2d.cpp:1163
ML2DBaseProgram::psi_max
double psi_max
Definition: ml2d.h:160
WEIGHT
#define WEIGHT
Definition: ml_align2d.cpp:1162
ProgML2D::Fimg_flip
std::vector< MultidimArray< std::complex< double > > > Fimg_flip
Definition: ml_align2d.h:90
ProgML2D::barrier
barrier_t barrier
Definition: ml_align2d.h:68
ProgML2D::expectation
virtual void expectation()
Integrate over all experimental images.
Definition: ml_align2d.cpp:1713
ML2DBaseProgram::defineAdditionalParams
virtual void defineAdditionalParams(XmippProgram *prog, const char *sectionLine)
Definition: ml2d.cpp:266
ProgML2D::sumfracweight
double sumfracweight
Definition: ml_align2d.h:73
FN_CLASSES_MD
#define FN_CLASSES_MD(base)
Definition: ml2d.h:55
x
doublereal * x
Definition: numerical_recipes.cpp:2230
ModelML2D::update
void update()
Definition: ml2d.cpp:414
i
#define i
Definition: numerical_recipes.cpp:2493
ProgML2D::refno_index
size_t refno_index
Definition: ml_align2d.h:110
MDL_ENABLED
Is this image enabled? (int [-1 or 1])
Definition: metadata_label.h:231
rotate
#define rotate(a, i, j, k, l)
computeStats_within_binary_mask
void computeStats_within_binary_mask(const MultidimArray< T1 > &mask, const MultidimArray< T > &m, double &min_val, double &max_val, double &avg, double &stddev)
Definition: mask.h:799
ProgML2D::mygroup
int mygroup
Definition: ml_align2d.h:113
FOR_ALL_ELEMENTS_IN_ARRAY2D
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
Definition: multidim_array_base.h:390
IMG_LOCAL_INDEX
#define IMG_LOCAL_INDEX
Definition: ml2d.h:67
XmippProgram::addSeeAlsoLine
void addSeeAlsoLine(const char *seeAlso)
Definition: xmipp_program.cpp:290
divide_equally_group
size_t divide_equally_group(size_t N, size_t size, size_t myself)
Definition: xmipp_funcs.cpp:1169
theta
double theta
Definition: numerical_recipes.cpp:6214
ProgML2D::ioptflip_ref
std::vector< int > ioptflip_ref
Definition: ml_align2d.h:97
ProgML2D::barrier3
barrier_t barrier3
Definition: ml_align2d.h:68
ModelML2D::print
void print(int tabs=0) const
Just for debugging now.
Definition: ml2d.cpp:475
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
ML2DBaseProgram::save_mem2
bool save_mem2
Definition: ml2d.h:192
ML2DBaseProgram::do_ML3D
bool do_ML3D
Definition: ml2d.h:196
MetaDataVec::clear
void clear() override
Definition: metadata_vec.cpp:170
ProgML2D::LL
double LL
Definition: ml_align2d.h:73
ML2DBaseProgram::nr_images_local
size_t nr_images_local
Definition: ml2d.h:166
first
glob_log first
Definition: numerical_recipes.cpp:5201
XmippProgram::argc
int argc
Original command line arguments.
Definition: xmipp_program.h:86
FOR_ALL_THREAD_REFNO
#define FOR_ALL_THREAD_REFNO()
Definition: ml_align2d.h:33
ML2DBaseProgram::search_rot
double search_rot
Definition: ml2d.h:190
ProgML2D::ioptx
int ioptx
Definition: ml_align2d.h:96
ProgML2D::wsumimgs
std::vector< MultidimArray< std::complex< double > > > wsumimgs
Definition: ml_align2d.h:83
ProgML2D::doThreadRotateReferenceRefno
void doThreadRotateReferenceRefno()
Thread code to parallelize refno loop in rotateReference.
Definition: ml_align2d.cpp:1058
ML2DBaseProgram::limit_rot
bool limit_rot
Definition: ml2d.h:188
ML2DBaseProgram::myFirstImg
size_t myFirstImg
Definition: ml2d.h:168
log
void log(Image< double > &op)
Definition: image_operate.cpp:228
getParameter
const char * getParameter(int argc, const char **argv, const char *param, const char *option)
Definition: args.cpp:30
ProgML2D::mindiff
double mindiff
Definition: ml_align2d.h:94
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
ProgML2D::refno_load
int refno_load
Definition: ml_align2d.h:111
FourierTransformer::setReal
void setReal(MultidimArray< double > &img)
Definition: xmipp_fftw.cpp:129
ProgML2D::bgmean
double bgmean
Definition: ml_align2d.h:85
ProgML2D::sumw
std::vector< double > sumw
Definition: ml_align2d.h:75
ProgML2D::doThreadReverseRotateReferenceRefno
void doThreadReverseRotateReferenceRefno()
Thread code to parallelize refno loop in reverseRotateReference.
Definition: ml_align2d.cpp:1119
index
viol index
Definition: numerical_recipes.cpp:5197
ML2DBaseProgram::imgs_oldtheta
std::vector< double > imgs_oldtheta
Definition: ml2d.h:194
selfTranslate
void selfTranslate(int SplineDegree, MultidimArray< T > &V1, const Matrix1D< double > &v, bool wrap=xmipp_transformation::WRAP, T outside=0)
Definition: transformations.h:1394
ProgML2D::opt_refno
int opt_refno
Definition: ml_align2d.h:84
MetaDataVec::setValue
bool setValue(const MDObject &mdValueIn, size_t id)
Definition: metadata_vec.cpp:282
CEIL
#define CEIL(x)
Definition: xmipp_macros.h:225
MetaDataVec::addObject
size_t addObject() override
Definition: metadata_vec.cpp:435
ProgML2D::maximizeModel
void maximizeModel(ModelML2D &model)
Update all model parameters.
Definition: ml_align2d.cpp:1933
ML2DBaseProgram::initSamplingStuff
void initSamplingStuff()
Definition: ml2d.cpp:44
MetaDataVec::firstRowId
size_t firstRowId() const override
Definition: metadata_vec.cpp:508
ProgML2D::trymindiff
double trymindiff
Definition: ml_align2d.h:85
ProgML2D::wsum_sigma_noise
double wsum_sigma_noise
Definition: ml_align2d.h:74
ML2DBaseProgram::defineBasicParams
virtual void defineBasicParams(XmippProgram *prog)
Definition: ml2d.cpp:226
MetaDataVec::importObjects
void importObjects(const MetaData &md, const std::vector< size_t > &objectsToAdd, bool doClear=true) override
Definition: metadata_vec.cpp:458
ML2DBaseProgram::F
std::vector< Matrix2D< double > > F
Definition: ml2d.h:174
ProgML2D::preselectFastSignificant
void preselectFastSignificant()
Definition: ml_align2d.cpp:674
SPECIAL_ITER
#define SPECIAL_ITER
Definition: ml_align2d.h:51
SMALLANGLE
#define SMALLANGLE
Definition: ml_tomo.h:48
MDL_FLIP
Flip the image? (bool)
Definition: metadata_label.h:236
ModelML2D::sumw_allrefs
double sumw_allrefs
Definition: ml2d.h:95
getIterExtraPath
FileName getIterExtraPath(const FileName &fn_root, int iter, bool makePath)
Definition: ml2d.cpp:305
ProgML2D::iopt_psi
int iopt_psi
Definition: ml_align2d.h:84
ML2DBaseProgram::do_restart
bool do_restart
Definition: ml2d.h:236
ProgML2D::refw_mirror
std::vector< double > refw_mirror
Definition: ml_align2d.h:91
ModelML2D::wsum_sigma_offset
double wsum_sigma_offset
Definition: ml2d.h:89
MDL_PMAX
Maximum value of normalized probability function (now called "Pmax/sumP") (double) ...
Definition: metadata_label.h:352
ModelML2D::sumw_allrefs2
double sumw_allrefs2
Definition: ml2d.h:95
max
void max(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:128
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ML2DBaseProgram::docfiledata
MultidimArray< double > docfiledata
Definition: ml2d.h:234
FN_CLASS_IMAGES_MD
#define FN_CLASS_IMAGES_MD(base, ref)
Definition: ml2d.h:57
XmippProgram::argv
const char ** argv
Definition: xmipp_program.h:87
ML2DBaseProgram::trymindiff_factor
double trymindiff_factor
Definition: ml2d.h:202
MDL_SIGMANOISE
Standard deviation of the noise in ML model.
Definition: metadata_label.h:445
ML2DBaseProgram::imgs_offsets
std::vector< std::vector< double > > imgs_offsets
Definition: ml2d.h:200
ITER_PREFIX
#define ITER_PREFIX
Definition: ml_align2d.cpp:2167
ML2DBaseProgram::istart
int istart
Definition: ml2d.h:145
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
ML2DBaseProgram::df
double df
Definition: ml2d.h:204
ProgML2D::defineParams
void defineParams()
Params definition.
Definition: ml_align2d.cpp:45
MSIGNIFICANT
#define MSIGNIFICANT
Definition: ml_align2d.cpp:1164
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
XmippProgram::addExampleLine
void addExampleLine(const char *example, bool verbatim=true)
Definition: xmipp_program.cpp:286
ROUND
#define ROUND(x)
Definition: xmipp_macros.h:210
ModelML2D::alpha_k
std::vector< double > alpha_k
Definition: ml2d.h:93
MDL_INTSCALE
Intensity scale for an image.
Definition: metadata_label.h:260
ProgML2D::wsum_sc2
double wsum_sc2
Definition: ml_align2d.h:93
ML2DBaseProgram::dim
size_t dim
Definition: ml2d.h:151
ModelML2D::avePmax
double avePmax
Definition: ml2d.h:97
ProgML2D::Xi2
double Xi2
Definition: ml_align2d.h:107
ProgML2D::maxweight2
double maxweight2
Definition: ml_align2d.h:86
ML2DBaseProgram::factor_nref
int factor_nref
Definition: ml2d.h:215
XmippProgram::verbose
int verbose
Verbosity level.
Definition: xmipp_program.h:161
FOR_ALL_THREAD_REFNO_NODECL
#define FOR_ALL_THREAD_REFNO_NODECL()
Same macro as before, but without declaring refno and load.
Definition: ml_align2d.h:38
ProgML2D::refno_count
int refno_count
Definition: ml_align2d.h:111
MetaDataDb::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const override
Definition: metadata_db.cpp:748
ML2DBaseProgram::imgs_oldphi
std::vector< double > imgs_oldphi
Definition: ml2d.h:194
IIFLIP
#define IIFLIP
Some macro definitions for the following function.
Definition: ml_align2d.cpp:1159
ML2DBaseProgram::MDref
MetaDataDb MDref
Definition: ml2d.h:172
ProgML2D::A2
std::vector< double > A2
Definition: ml_align2d.h:105
MetaDataDb::size
size_t size() const override
Definition: metadata_db.cpp:549
ProgML2D::doThreadPreselectFastSignificantRefno
void doThreadPreselectFastSignificantRefno()
Thread code to parallelize refno loop in preselectFastSignificant.
Definition: ml_align2d.cpp:1166
FourierTransformer::FourierTransform
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
ProgML2D::writeOutputFiles
virtual void writeOutputFiles(const ModelML2D &model, OutputType outputType=OUT_FINAL)
Write model parameters.
Definition: ml_align2d.cpp:2171
ML2DBaseProgram::imgs_optrefno
std::vector< int > imgs_optrefno
Definition: ml2d.h:211
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
ML2DBaseProgram::defineHiddenParams
virtual void defineHiddenParams(XmippProgram *prog)
Definition: ml2d.cpp:290
MDL_SIGMAOFFSET
Standard deviation of the offsets in ML model.
Definition: metadata_label.h:446
ML2DBaseProgram::fix_sigma_noise
bool fix_sigma_noise
Definition: ml2d.h:143
ProgML2D::pdf_directions
std::vector< double > pdf_directions
Definition: ml_align2d.h:80
ProgML2D::sumwsc2
std::vector< double > sumwsc2
Definition: ml_align2d.h:75
j
#define j
Definition: numerical_recipes.cpp:2493
MDL_REF3D
3D Class to which the image belongs (int)
Definition: metadata_label.h:376
ProgML2D::mask
MultidimArray< int > mask
Definition: ml_align2d.h:65
ProgML2D::refwsc2
std::vector< double > refwsc2
Definition: ml_align2d.h:91
ProgML2D::sum_refw2
double sum_refw2
Definition: ml_align2d.h:92
ProgML2D::wsum_offset
double wsum_offset
Definition: ml_align2d.h:93
ml_align2d.h
ML2DBaseProgram::df2
double df2
Definition: ml2d.h:204
ProgML2D::mref
std::vector< MultidimArray< double > > mref
Definition: ml_align2d.h:81
ERR_THREADS_NOTINIT
Threads cannot be initiated.
Definition: xmipp_error.h:188
MetaDataVec::getValue
bool getValue(MDObject &mdValueOut, size_t id) const override
Definition: metadata_vec.cpp:293
ImageBase::setEulerAngles
void setEulerAngles(double rot, double tilt, double psi, const size_t n=0)
Definition: xmipp_image_base.cpp:509
ProgML2D::allref_offsets
std::vector< double > allref_offsets
Definition: ml_align2d.h:79
ML2DBaseProgram::fn_scratch
FileName fn_scratch
Definition: ml2d.h:132
ProgML2D::refw2
std::vector< double > refw2
Definition: ml_align2d.h:91
MDL_REF
Class to which the image belongs (int)
Definition: metadata_label.h:377
ASSIGN
#define ASSIGN(var)
LOG
#define LOG(msg)
Definition: mpi_write_test.cpp:49
OUT_BLOCK
Definition: ml2d.h:76
XmippProgram::setProgress
void setProgress(size_t value=0)
Definition: xmipp_program.cpp:263
MultidimArray::sum2
double sum2() const
Definition: multidim_array.h:3424
ML2DBaseProgram::setNumberOfLocalImages
virtual void setNumberOfLocalImages()
Set the number of images, this function is useful only for MPI.
Definition: ml2d.cpp:110
ML2DBaseProgram::fn_frac
FileName fn_frac
Definition: ml2d.h:132
mybarrier_t
Definition: xmipp_threads.h:472
ProgML2D::opt_psi
double opt_psi
Definition: ml_align2d.h:85
ProgML2D::maxweight
double maxweight
Definition: ml_align2d.h:92
alglib::round
int round(double x)
Definition: ap.cpp:7245
MetaDataDb::setValue
bool setValue(const MDObject &mdValueIn, size_t id) override
Definition: metadata_db.cpp:90
translate
void translate(PDB *pdb_original, PDB *pdb_move, unsigned num_atoms, double x0, double y0, double z0)
Definition: lib_pwk.cpp:125
ProgML2D::calculatePdfInplane
void calculatePdfInplane()
Calculate probability density distribution for in-plane transformations.
Definition: ml_align2d.cpp:553
MDL_CLASS_COUNT
Number of images assigned to the same class as this image.
Definition: metadata_label.h:81
DP
#define DP(x)
FN_IMAGES_MD
#define FN_IMAGES_MD(base)
Definition: ml2d.h:54
ProgML2D::iopty
int iopty
Definition: ml_align2d.h:96
ModelML2D::sigma_noise
double sigma_noise
Definition: ml2d.h:87
ModelML2D::Iref
std::vector< Image< double > > Iref
Definition: ml2d.h:85
dotProduct
T dotProduct(const Matrix1D< T > &v1, const Matrix1D< T > &v2)
Definition: matrix1d.h:1140
MetaData::removeDisabled
virtual void removeDisabled()
Definition: metadata_base.cpp:123
TH_EXIT
Definition: ml2d.h:74
RAD2DEG
#define RAD2DEG(r)
Definition: xmipp_macros.h:320
ProgML2D::rotateReference
void rotateReference()
Fill vector of matrices with all rotations of reference.
Definition: ml_align2d.cpp:603
ML2DBaseProgram::C_fast
double C_fast
Definition: ml2d.h:182
ModelML2D::sumfracweight
double sumfracweight
Definition: ml2d.h:97
ProgML2D::no_iem
bool no_iem
Definition: ml_align2d.h:63
ML2DBaseProgram::fn_ref
FileName fn_ref
Definition: ml2d.h:132
FourierTransformer::getFourierAlias
void getFourierAlias(T &V)
Definition: xmipp_fftw.h:207
ProgML2D::reverseRotateReference
void reverseRotateReference()
Apply reverse rotations to all matrices in vector and fill new matrix with their sum.
Definition: ml_align2d.cpp:618
ProgML2D::expectationSingleImage
void expectationSingleImage(Matrix1D< double > &opt_offsets)
ML-integration over all (or -fast) translations.
Definition: ml_align2d.cpp:700
ModelML2D::sigma_offset
double sigma_offset
Definition: ml2d.h:89
update_mutex
pthread_mutex_t update_mutex
Definition: ml_align2d.cpp:31
ProgML2D::sumwsc
std::vector< double > sumwsc
Definition: ml_align2d.h:75
ProgML2D::wsum_Mref
std::vector< MultidimArray< double > > wsum_Mref
Definition: ml_align2d.h:76
ProgML2D::doThreadESIUpdateRefno
void doThreadESIUpdateRefno()
Thread code to parallelize update loop in ESI.
Definition: ml_align2d.cpp:1624
String
std::string String
Definition: xmipp_strings.h:34
alglib::psi
double psi(const double x)
Definition: specialfunctions.cpp:2492
ModelML2D::do_student
bool do_student
Definition: ml2d.h:104
ProgML2D::refno_load_param
size_t refno_load_param
Definition: ml_align2d.h:110
ProgML2D::awakeThreads
void awakeThreads(ThreadTask task, int start_refno, int load=1)
Awake threads for different tasks.
Definition: ml_align2d.cpp:1046
FN_CLASSES_STK
#define FN_CLASSES_STK(base)
Definition: ml2d.h:58
doThreadsTasks
void * doThreadsTasks(void *data)
Function for threads do different tasks.
Definition: ml_align2d.cpp:967
ProgML2D::refw
std::vector< double > refw
Definition: ml_align2d.h:91
formatString
String formatString(const char *format,...)
Definition: xmipp_strings.cpp:602
ModelML2D
Definition: ml2d.h:79
MetaDataDb::read
void read(const FileName &inFile, const std::vector< MDLabel > *desiredLabels=NULL, bool decomposeStack=true) override
Definition: metadata_db.cpp:912
ProgML2D::sumw_refpsi
std::vector< double > sumw_refpsi
Definition: ml_align2d.h:91
realWRAP
#define realWRAP(x, x0, xF)
Definition: xmipp_macros.h:304
ML2DBaseProgram::P_phi
MultidimArray< double > P_phi
Definition: ml2d.h:178
u
doublereal * u
Definition: numerical_recipes.cpp:2230
ML2DBaseProgram::outRefsMd
String outRefsMd
Definition: ml2d.h:250
ML2DBaseProgram::randomizeImagesOrder
virtual void randomizeImagesOrder()
Randomize initial images order, only once.
Definition: ml2d.cpp:96
ML2DBaseProgram::do_mirror
bool do_mirror
Definition: ml2d.h:137
ProgML2D::th_ids
pthread_t * th_ids
Definition: ml_align2d.h:69
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
MetaDataDb::unionAll
void unionAll(const MetaDataDb &mdIn)
Definition: metadata_db.cpp:1366
OUTSIDE_MASK
constexpr int OUTSIDE_MASK
Definition: mask.h:48
ML2DBaseProgram::fast_mode
bool fast_mode
Definition: ml2d.h:180
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
ML2DBaseProgram::fix_sigma_offset
bool fix_sigma_offset
Definition: ml2d.h:141
r2
float r2
Definition: image_find_center.cpp:54
ML2DBaseProgram::nr_images_global
size_t nr_images_global
Definition: ml2d.h:164
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
ML2DBaseProgram::imgs_trymindiff
std::vector< double > imgs_trymindiff
Definition: ml2d.h:209
ProgML2D::correctScaleAverage
void correctScaleAverage()
Correct references scale.
Definition: ml_align2d.cpp:2085
ModelML2D::sumw_mirror
std::vector< double > sumw_mirror
Definition: ml2d.h:91
MDL_MIRRORFRAC
Mirror fraction for a Maximum Likelihood model.
Definition: metadata_label.h:293
structThreadTasks::prm
ProgML2D * prm
Definition: ml_align2d.h:47
idum
int idum
Definition: xmipp_funcs.cpp:391
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
MultidimArray::computeAvg
double computeAvg() const
Definition: multidim_array.h:1739
DATALINELENGTH
#define DATALINELENGTH
Definition: ml2d.h:50
ModelML2D::WsumMref
std::vector< Image< double > > WsumMref
Definition: ml2d.h:85
ProgML2D::getBaseName
FileName getBaseName(String suffix="", int number=-1)
Get base name based on fn_root and some number.
Definition: ml_align2d.cpp:2378
ImageBase::setWeight
void setWeight(double weight, const size_t n=0)
Definition: xmipp_image_base.h:770
OUT_FINAL
Definition: ml2d.h:76
ProgML2D::show
virtual void show()
Show info at starting program.
Definition: ml_align2d.cpp:196
ML2DBaseProgram::save_mem1
bool save_mem1
Definition: ml2d.h:192
ProgML2D::pfs_weight
MultidimArray< double > pfs_weight
Definition: ml_align2d.h:101
ML2DBaseProgram::dim2
size_t dim2
Definition: ml2d.h:151
barrier_wait
int barrier_wait(barrier_t *barrier)
Definition: xmipp_threads.cpp:382
PI
#define PI
Definition: tools.h:43
MDL_ITER
Current iteration number (int)
Definition: metadata_label.h:262
ProgML2D::iteration
virtual void iteration()
Perform an iteration.
Definition: ml_align2d.cpp:1701
ML2DBaseProgram::nr_nomirror_flips
size_t nr_nomirror_flips
Definition: ml2d.h:162
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331
ProgML2D::fracweight
double fracweight
Definition: ml_align2d.h:86
ERR_VALUE_INCORRECT
Incorrect value received.
Definition: xmipp_error.h:195
ML2DBaseProgram::Iold
std::vector< Image< double > > Iold
Definition: ml2d.h:176
refno_mutex
pthread_mutex_t refno_mutex
Definition: ml_align2d.cpp:34
ProgML2D::mysumimgs
std::vector< MultidimArray< std::complex< double > > > mysumimgs
Definition: ml_align2d.h:90
ML2DBaseProgram::conv
std::vector< double > conv
Definition: ml2d.h:219
ML2DBaseProgram::blocks
size_t blocks
Definition: ml2d.h:227
ModelML2D::scale
std::vector< double > scale
Definition: ml2d.h:93
n
int * n
Definition: numerical_recipes.cpp:2229
MDL_IMAGE
Name of an image (std::string)
Definition: metadata_label.h:246
OUT_IMGS
Definition: ml2d.h:76
ML2DBaseProgram::threads
int threads
Definition: ml2d.h:244
ML2DBaseProgram::imgs_bgmean
std::vector< double > imgs_bgmean
Definition: ml2d.h:209
ModelML2D::do_norm
bool do_norm
Definition: ml2d.h:104
gammln
double gammln(double xx)
Definition: numerical_recipes.cpp:589
ML2DBaseProgram::current_image
size_t current_image
Definition: ml2d.h:221
ModelML2D::addModel
void addModel(const ModelML2D &model)
Definition: ml2d.cpp:404
ML2DBaseProgram::eps
double eps
Definition: ml2d.h:170
ProgML2D::Msignificant
MultidimArray< int > Msignificant
Definition: ml_align2d.h:77
ProgML2D::sumw2
std::vector< double > sumw2
Definition: ml_align2d.h:75
ProgML2D::produceSideInfo
virtual void produceSideInfo()
Try to merge produceSideInfo1 and 2.
Definition: ml_align2d.cpp:290
MDL_WEIGHT
< Score 4 for volumes
Definition: metadata_label.h:480
ML2DBaseProgram::imgs_scale
std::vector< double > imgs_scale
Definition: ml2d.h:209
xmipp_transformation::BSPLINE3
Definition: transformations_defines.h:32
ProgML2D::pfs_maxweight
MultidimArray< double > pfs_maxweight
Definition: ml_align2d.h:100
ML2DBaseProgram::cline
String cline
Definition: ml2d.h:134
ModelML2D::wsum_sigma_noise
double wsum_sigma_noise
Definition: ml2d.h:87
INNER_MASK
constexpr int INNER_MASK
Definition: mask.h:47
MDL_RANDOMSEED
Seed for random number generator.
Definition: metadata_label.h:375
Generated by   doxygen 1.8.13