Xmipp  v3.23.11-Nereus
tomo_align_tilt_series.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Carlos Oscar coss@cnb.csic.es (2008)
4  *
5  * Unidad de Bioinformatica of 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.csic.es'
24  ***************************************************************************/
25 
26 #include <fstream>
27 #include <queue>
28 #include "tomo_align_tilt_series.h"
29 #include "core/metadata_sql.h"
30 #include "core/xmipp_image.h"
31 #include "core/transformations.h"
32 #include "core/utils/memory_utils.h"
33 #include "data/fourier_filter.h"
34 #include "data/mask.h"
35 #include "data/morphology.h"
36 #include "data/numerical_tools.h"
37 
38 /* Generate mask ----------------------------------------------------------- */
39 //#define DEBUG
40 void generateMask(const MultidimArray<double> &I, MultidimArray<unsigned char> &mask,
41  int patchSize)
42 {
43  MultidimArray<double> dmask;
44  dmask.initZeros(I);
45  FOR_ALL_ELEMENTS_IN_ARRAY2D(dmask)
46  if (I(i,j)!=0)
47  dmask(i,j)=1;
48 
49  MultidimArray<double> maskEroded;
50  maskEroded.initZeros(dmask);
51  erode2D(dmask,maskEroded,8,0,patchSize);
52  typeCast(maskEroded,mask);
53  mask.setXmippOrigin();
54 
55 #ifdef DEBUG
56 
57  ImageXmipp save;
58  save()=I;
59  save.write("PPP.xmp");
60  save()=maskEroded;
61  save.write("PPPmask.xmp");
62  std::cout << "Masks saved. Press any key\n";
63  char c;
64  std::cin >> c;
65 #endif
66 }
67 #undef DEBUG
68 
69 /* Compute affine matrix --------------------------------------------------- */
70 class AffineFitness
71 {
72 public:
73  Matrix1D<double> minAllowed;
74  Matrix1D<double> maxAllowed;
75  std::vector< MultidimArray<double> *> I1;
76  std::vector< MultidimArray<double> *> I2;
77  std::vector< MultidimArray<double> *> Mask1;
78  std::vector< MultidimArray<double> *> Mask2;
79  bool showMode;
80  bool checkRotation;
81 
82  AffineFitness()
83  {
84  showMode=false;
85  checkRotation=true;
86  }
87 
88  ~AffineFitness()
89  {
90  for (size_t i=0; i<I1.size(); i++)
91  delete I1[i];
92  for (size_t i=0; i<I2.size(); i++)
93  delete I2[i];
94  for (size_t i=0; i<Mask1.size(); i++)
95  delete Mask1[i];
96  for (size_t i=0; i<Mask2.size(); i++)
97  delete Mask2[i];
98  }
99 
100  double affine_fitness_individual(double *p)
101  {
102  // Check limits
103  if (!showMode)
104  FOR_ALL_ELEMENTS_IN_MATRIX1D(minAllowed)
105  if (p[i]<minAllowed(i) || p[i]>maxAllowed(i))
106  return 1e20;
107 
108  // Separate solution
109  Matrix2D<double> A12, A21;
110  A12.initIdentity(3);
111  A12(0,0)=p[0];
112  A12(0,1)=p[1];
113  A12(0,2)=p[4];
114  A12(1,0)=p[2];
115  A12(1,1)=p[3];
116  A12(1,2)=p[5];
117 
118  A21=A12.inv();
119 
120  // Check it is approximately a rotation
121  if (!showMode && checkRotation)
122  {
123  // Check that the eigenvalues are close to 1
124  std::complex<double> a=A12(0,0);
125  std::complex<double> b=A12(0,1);
126  std::complex<double> c=A12(1,0);
127  std::complex<double> d=A12(1,1);
128  std::complex<double> eig1=(a+d)/2.0+sqrt(4.0*b*c+(a-d)*(a-d))/2.0;
129  std::complex<double> eig2=(a+d)/2.0-sqrt(4.0*b*c+(a-d)*(a-d))/2.0;
130  double abs_eig1=abs(eig1);
131  double abs_eig2=abs(eig2);
132 
133  if (abs_eig1<0.95 || abs_eig1>1.05)
134  return 1e20*ABS(abs_eig1-1);
135  if (abs_eig2<0.95 || abs_eig2>1.05)
136  return 1e20*ABS(abs_eig2-1);
137 
138  // Compute the eigenvalues of A*A^t disregarding the shifts
139  // A*A^t=[a b; c d]
140 
141  a=A12(0,0)*A12(0,0)+A12(0,1)*A12(0,1);
142  b=A12(0,0)*A12(1,0)+A12(0,1)*A12(1,1);
143  c=b;
144  d=A12(1,0)*A12(1,0)+A12(1,1)*A12(1,1);
145  eig1=(a+d)/2.0+sqrt(4.0*b*c+(a-d)*(a-d))/2.0;
146  eig2=(a+d)/2.0-sqrt(4.0*b*c+(a-d)*(a-d))/2.0;
147  abs_eig1=abs(eig1);
148  abs_eig2=abs(eig2);
149 
150  if (abs_eig1<0.95 || abs_eig1>1.05)
151  return 1e20*ABS(abs_eig1-1);
152  if (abs_eig2<0.95 || abs_eig2>1.05)
153  return 1e20*ABS(abs_eig2-1);
154  if (abs(c)>0.05)
155  return 1e20*ABS(abs(c));
156  }
157 
158  // For each pyramid level
159  double dist=0;
160  Matrix2D<double> A12level=A12;
161  Matrix2D<double> A21level=A21;
162  for (size_t level=0; level<I1.size(); level++)
163  {
164  const MultidimArray<double> &Mask1_level=*(Mask1[level]);
165  const MultidimArray<double> &Mask2_level=*(Mask2[level]);
166 
167  // Produce the transformed images
168  MultidimArray<double> transformedI1, transformedI2;
169  applyGeometry(xmipp_transformation::LINEAR,transformedI1,*(I1[level]),A12level,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
170  applyGeometry(xmipp_transformation::LINEAR,transformedI2,*(I2[level]),A21level,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
171 
172  // Produce masks for the comparison
173  MultidimArray<int> maskInTheSpaceOf1, maskInTheSpaceOf2;
174  MultidimArray<double> maskAux;
175  applyGeometry(xmipp_transformation::LINEAR,maskAux,Mask1_level,A12level,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
176  maskInTheSpaceOf2.initZeros(YSIZE(maskAux),XSIZE(maskAux));
177  maskInTheSpaceOf2.setXmippOrigin();
178  FOR_ALL_ELEMENTS_IN_ARRAY2D(maskAux)
179  {
180  if (A2D_ELEM(maskAux,i,j)>0.5)
181  A2D_ELEM(maskInTheSpaceOf2,i,j)=(int)A2D_ELEM(Mask2_level,i,j);
182  }
183 
184  maskAux.initZeros();
185  applyGeometry(xmipp_transformation::LINEAR,maskAux,Mask2_level,A21level,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
186  maskInTheSpaceOf1.initZeros(YSIZE(maskAux),XSIZE(maskAux));
187  maskInTheSpaceOf1.setXmippOrigin();
188  FOR_ALL_ELEMENTS_IN_ARRAY2D(maskAux)
189  {
190  if (A2D_ELEM(maskAux,i,j)>0.5)
191  A2D_ELEM(maskInTheSpaceOf1,i,j)=(int)A2D_ELEM(Mask1_level,i,j);
192  }
193 
194  // Compare the two images
195  double distLevel=0.5*(
196  1-correlationIndex(transformedI1,*(I2[level]),&maskInTheSpaceOf2)+
197  1-correlationIndex(transformedI2,*(I1[level]),&maskInTheSpaceOf1));
198 
199  if (showMode)
200  {
201  Image<double> save;
202  save()=*(I1[level]);
203  save.write("PPPimg1.xmp");
204  save()=*(I2[level]);
205  save.write("PPPimg2.xmp");
206  save()=*(Mask1[level]);
207  save.write("PPPmask1.xmp");
208  save()=*(Mask2[level]);
209  save.write("PPPmask2.xmp");
210  save()=transformedI1;
211  save.write("PPPTransformedImg1.xmp");
212  save()=transformedI2;
213  save.write("PPPTransformedImg2.xmp");
214  typeCast(maskInTheSpaceOf1,save());
215  save.write("PPPmaskInTheSpaceOf1.xmp");
216  typeCast(maskInTheSpaceOf2,save());
217  save.write("PPPmaskInTheSpaceOf2.xmp");
218  save()=*(I1[level])-transformedI2;
219  save.write("PPPDiffImg1.xmp");
220  save()=*(I2[level])-transformedI1;
221  save.write("PPPDiffImg2.xmp");
222  std::cout << "Level=" << level << "\nA12level=\n" << A12level << "A21level=\n" << A21level << std::endl;
223  std::cout << "distLevel=" << distLevel << std::endl;
224  std::cout << "Do you like the alignment? (y/n)\n";
225  char c;
226  std::cin >> c;
227  if (c=='n')
228  {
229  dist=-1;
230  break;
231  }
232  }
233  dist+=distLevel;
234 
235  A12level(0,2)/=2;
236  A12level(1,2)/=2;
237  A21level(0,2)/=2;
238  A21level(1,2)/=2;
239  }
240  dist/=I1.size();
241 
242  return dist;
243  }
244 
245  static double Powell_affine_fitness_individual(double *p, void *prm)
246  {
247  return ((AffineFitness*)prm)->affine_fitness_individual(p+1);
248  }
249 };
250 
251 void setupAffineFitness(AffineFitness &fitness, const MultidimArray<double> &I1,
252  const MultidimArray<double> &I2, int maxShift, bool isMirror,
253  bool checkRotation, int pyramidLevel)
254 {
255  fitness.checkRotation=checkRotation;
256 
257  // Set images
258  int level=0;
259  MultidimArray<double> I1aux=I1;
260  MultidimArray<double> I2aux=I2;
261 
262  // Remove the borders
263  int borderY=CEIL(0.1*YSIZE(I1));
264  int borderX=CEIL(0.1*XSIZE(I1));
265  I1aux.selfWindow(STARTINGY(I1)+borderY,STARTINGX(I1)+borderX,
266  FINISHINGY(I1)-borderY,FINISHINGX(I1)-borderX);
267  I2aux.selfWindow(STARTINGY(I1)+borderY,STARTINGX(I1)+borderX,
268  FINISHINGY(I1)-borderY,FINISHINGX(I1)-borderX);
269 
270  MultidimArray<double> Mask1;
271  Mask1.initZeros(I1aux);
272  FOR_ALL_ELEMENTS_IN_ARRAY2D(I1aux)
273  if (A2D_ELEM(I1,i,j)!=0)
274  A2D_ELEM(Mask1,i,j)=1;
275 
276  MultidimArray<double> Mask2;
277  Mask2.initZeros(I2aux);
278  FOR_ALL_ELEMENTS_IN_ARRAY2D(I2aux)
279  if (A2D_ELEM(I2,i,j)!=0)
280  A2D_ELEM(Mask2,i,j)=1;
281 
282  do
283  {
284  // Push back the current images
285  I1aux.setXmippOrigin();
286  I2aux.setXmippOrigin();
287  Mask1.setXmippOrigin();
288  Mask2.setXmippOrigin();
289 
290  MultidimArray<double> *dummy=nullptr;
291  dummy=new MultidimArray<double>;
292  *dummy=I1aux;
293  fitness.I1.push_back(dummy);
294  dummy=new MultidimArray<double>;
295  *dummy=I2aux;
296  fitness.I2.push_back(dummy);
297  dummy=new MultidimArray<double>;
298  *dummy=Mask1;
299  fitness.Mask1.push_back(dummy);
300  dummy=new MultidimArray<double>;
301  *dummy=Mask2;
302  fitness.Mask2.push_back(dummy);
303 
304  // Prepare for next level
305  level++;
306  if (pyramidLevel>=level)
307  {
308  selfScaleToSize(xmipp_transformation::LINEAR,I1aux,YSIZE(I1aux)/2,XSIZE(I1aux)/2);
309  selfScaleToSize(xmipp_transformation::LINEAR,I2aux,YSIZE(I2aux)/2,XSIZE(I2aux)/2);
310  selfScaleToSize(xmipp_transformation::LINEAR,Mask1,YSIZE(Mask1)/2,XSIZE(Mask1)/2);
311  selfScaleToSize(xmipp_transformation::LINEAR,Mask2,YSIZE(Mask2)/2,XSIZE(Mask2)/2);
312  FOR_ALL_ELEMENTS_IN_ARRAY2D(Mask1)
313  Mask1(i,j)=(Mask1(i,j)>0.5)? 1:0;
314  FOR_ALL_ELEMENTS_IN_ARRAY2D(Mask2)
315  Mask2(i,j)=(Mask2(i,j)>0.5)? 1:0;
316  }
317  }
318  while (level<=pyramidLevel);
319 
320  // Set limits for the affine matrices
321  // Order: 1->2: 4 affine params+2 translations
322  // Order: 2->1: 4 affine params+2 translations
323  fitness.minAllowed.resize(6);
324  fitness.maxAllowed.resize(6);
325 
326  // Scale factors
327  fitness.minAllowed(0)=fitness.minAllowed(3)=0.9;
328  fitness.maxAllowed(0)=fitness.maxAllowed(3)=1.1;
329  if (isMirror)
330  {
331  fitness.minAllowed(3)=-1.1;
332  fitness.maxAllowed(3)=-0.9;
333  }
334 
335  // Rotation factors
336  fitness.minAllowed(1)=fitness.minAllowed(2)=-0.2;
337  fitness.maxAllowed(1)=fitness.maxAllowed(2)= 0.2;
338 
339  // Shifts
340  fitness.minAllowed(4)=fitness.minAllowed(5)=-maxShift;
341  fitness.maxAllowed(4)=fitness.maxAllowed(5)= maxShift;
342 }
343 
344 static pthread_mutex_t globalAffineMutex = PTHREAD_MUTEX_INITIALIZER;
345 double computeAffineTransformation(const MultidimArray<unsigned char> &I1,
346  const MultidimArray<unsigned char> &I2, int maxShift, int maxIterDE,
347  Matrix2D<double> &A12, Matrix2D<double> &A21, bool show,
348  double thresholdAffine, bool globalAffine, bool isMirror,
349  bool checkRotation, int pyramidLevel)
350 {
351  try
352  {
353  AffineFitness affy;
354  MultidimArray<double> I1d, I2d;
355  typeCast(I1, I1d);
356  typeCast(I2, I2d);
357  setupAffineFitness(affy, I1d, I2d, maxShift, isMirror, checkRotation,
358  pyramidLevel);
359 
360  // Return result
361  double cost = 0.0;
362 
363  // Optimize with differential evolution
364  Matrix1D<double> A(6);
365  A(0)=A(3)=1;
366  if (isMirror)
367  A(3)*=-1;
368  if (MAT_XSIZE(A12)==0)
369  {
370  if (globalAffine)
371  {
372  // Exhaustive search
373  double bestCost=2;
374  double stepX=(affy.maxAllowed(4)-affy.minAllowed(4))/40.0;
375  double stepY=(affy.maxAllowed(5)-affy.minAllowed(5))/40.0;
376  for (double shiftY=affy.minAllowed(5); shiftY<=affy.maxAllowed(5); shiftY+=stepY)
377  for (double shiftX=affy.minAllowed(4); shiftX<=affy.maxAllowed(4); shiftX+=stepX)
378  {
379  A(4)=shiftX;
380  A(5)=shiftY;
381  double cost=affy.affine_fitness_individual(MATRIX1D_ARRAY(A));
382  if (cost<bestCost)
383  bestCost=cost;
384  }
385  }
386  else
387  {
388  // Initialize with cross correlation
389  double tx, ty;
390  pthread_mutex_lock( &globalAffineMutex );
391  CorrelationAux aux;
392  if (!isMirror)
393  bestShift(I1d,I2d,tx,ty,aux);
394  else
395  {
396  MultidimArray<double> auxI2d=I2d;
397  auxI2d.selfReverseY();
398  STARTINGX(auxI2d)=STARTINGX(I2d);
399  STARTINGY(auxI2d)=STARTINGY(I2d);
400  bestShift(I1d,auxI2d,tx,ty,aux);
401  ty=-ty;
402  }
403  pthread_mutex_unlock( &globalAffineMutex );
404  A(4)=-tx;
405  A(5)=-ty;
406  }
407 
408  // Optimize with Powell
409  Matrix1D<double> steps(A);
410  steps.initConstant(1);
411  int iter;
412  powellOptimizer(A, 1, VEC_XSIZE(A),
413  AffineFitness::Powell_affine_fitness_individual, &affy, 0.005,
414  cost, iter, steps, false);
415 
416  // Separate solution
417  A12.initIdentity(3);
418  A12(0,0)=A(0);
419  A12(0,1)=A(1);
420  A12(0,2)=A(4);
421  A12(1,0)=A(2);
422  A12(1,1)=A(3);
423  A12(1,2)=A(5);
424 
425  A21=A12.inv();
426  }
427 
428  if (show)
429  {
430  affy.showMode=true;
431  Matrix1D<double> p(6);
432  p(0)=A12(0,0);
433  p(1)=A12(0,1);
434  p(4)=A12(0,2);
435  p(2)=A12(1,0);
436  p(3)=A12(1,1);
437  p(5)=A12(1,2);
438  cost = affy.affine_fitness_individual(MATRIX1D_ARRAY(p));
439  affy.showMode=false;
440  }
441 
442  return cost;
443  }
444  catch (XmippError &XE)
445  {
446  std::cout << XE;
447  exit(1);
448  }
449 }
450 
451 /* Parameters -------------------------------------------------------------- */
452 void ProgTomographAlignment::readParams()
453 {
454  fnSel=getParam("-i");
455  fnSelOrig=getParam("--iorig");
456  fnRoot=getParam("--oroot");
457  if (fnRoot=="")
458  fnRoot=fnSel.withoutExtension();
459  globalAffine=checkParam("--globalAffine");
460  useCriticalPoints=checkParam("--useCriticalPoints");
461  if (useCriticalPoints)
462  Ncritical=getIntParam("--useCriticalPoints");
463  else
464  Ncritical=0;
465  seqLength=getIntParam("--seqLength");
466  blindSeqLength=getIntParam("--blindSeqLength");
467  maxStep=getIntParam("--maxStep");
468  gridSamples=getIntParam("--gridSamples");
469  psiMax=getDoubleParam("--psiMax");
470  deltaRot=getDoubleParam("--deltaRot");
471  localSize=getDoubleParam("--localSize");
472  optimizeTiltAngle=checkParam("--optimizeTiltAngle");
473  isCapillar=checkParam("--isCapillar");
474  dontNormalize=checkParam("--dontNormalize");
475  difficult=checkParam("--difficult");
476  corrThreshold=getDoubleParam("--threshold");
477  maxShiftPercentage=getDoubleParam("--maxShiftPercentage");
478  maxIterDE=getIntParam("--maxIterDE");
479  showAffine=checkParam("--showAffine");
480  thresholdAffine=getDoubleParam("--thresholdAffine");
481  identifyOutliersZ = getDoubleParam("--identifyOutliers");
482  doNotIdentifyOutliers = checkParam("--noOutliers");
483  pyramidLevel = getIntParam("--pyramid");
484  numThreads = getIntParam("--thr");
485  if (numThreads<1)
486  numThreads = 1;
487  lastStep=getIntParam("--lastStep");
488 }
489 
490 void ProgTomographAlignment::show()
491 {
492  std::cout << "Input images: " << fnSel << std::endl
493  << "Original images: " << fnSelOrig << std::endl
494  << "Output rootname: " << fnRoot << std::endl
495  << "Global affine: " << globalAffine << std::endl
496  << "Use critical points:" << useCriticalPoints << std::endl
497  << "Num critical points:" << Ncritical << std::endl
498  << "SeqLength: " << seqLength << std::endl
499  << "BlindSeqLength: " << blindSeqLength << std::endl
500  << "MaxStep: " << maxStep << std::endl
501  << "Grid samples: " << gridSamples << std::endl
502  << "Maximum psi: " << psiMax << std::endl
503  << "Delta rot: " << deltaRot << std::endl
504  << "Local size: " << localSize << std::endl
505  << "Optimize tilt angle:" << optimizeTiltAngle << std::endl
506  << "isCapillar: " << isCapillar << std::endl
507  << "DontNormalize: " << dontNormalize << std::endl
508  << "Difficult: " << difficult << std::endl
509  << "Threshold: " << corrThreshold << std::endl
510  << "MaxShift Percentage:" << maxShiftPercentage << std::endl
511  << "MaxIterDE: " << maxIterDE << std::endl
512  << "Show Affine: " << showAffine << std::endl
513  << "Threshold Affine: " << thresholdAffine << std::endl
514  << "Identify outliers Z:" << identifyOutliersZ << std::endl
515  << "No outliers: " << doNotIdentifyOutliers << std::endl
516  << "Pyramid level: " << pyramidLevel << std::endl
517  << "Threads to use: " << numThreads << std::endl
518  << "Last step: " << lastStep << std::endl
519  ;
520 }
521 
522 void ProgTomographAlignment::defineParams()
523 {
524  addUsageLine("Align a single-axis tilt series without any marker.");
525  addSeeAlsoLine("tomo_align_dual_tilt_series");
526  addParamsLine(" == General Options == ");
527  addParamsLine(" -i <metadatafile> : Input images");
528  addParamsLine(" : The selfile must contain the list of micrographs");
529  addParamsLine(" : and its tilt angles");
530  addParamsLine(" [--iorig <metadatafile=\"\">] : Metadata with images at original scale");
531  addParamsLine(" [--oroot <fn_out=\"\">] : Output alignment");
532  addParamsLine(" : If not given, the input selfile without extension");
533  addParamsLine(" [--thr <num=1>] : Parallel processing using \"num\" threads");
534  addParamsLine(" [--lastStep+ <step=-1>] : Last step to perform");
535  addParamsLine(" : Step -1 -> Perform all steps");
536  addParamsLine(" : Step 0 -> Determination of affine transformations");
537  addParamsLine(" : Step 1 -> Determination of landmark chains");
538  addParamsLine(" : Step 2 -> Determination of alignment parameters");
539  addParamsLine(" : Step 3 -> Writing aligned images");
540  addParamsLine(" == Step 0 (Affine alignment) Options == ");
541  addParamsLine(" [--maxShiftPercentage <p=0.2>] : Maximum shift as percentage of image size");
542  addParamsLine(" [--thresholdAffine <th=0.85>] : Threshold affine");
543  addParamsLine(" [--globalAffine] : Look globally for affine transformations");
544  addParamsLine(" [--difficult+] : Apply some filters before affine alignment");
545  addParamsLine(" [--maxIterDE+ <n=30>] : Maximum number of iteration in Differential Evolution");
546  addParamsLine(" [--showAffine+] : Show affine transformations as PPP*");
547  addParamsLine(" [--identifyOutliers+ <z=5>] : Z-score to be an outlier");
548  addParamsLine(" [--noOutliers+] : Do not identify outliers");
549  addParamsLine(" [--pyramid+ <level=1>] : Multiresolution for affine transformations");
550  addParamsLine(" == Step 1 (Landmark chain) Options == ");
551  addParamsLine(" [--seqLength <n=5>] : Sequence length");
552  addParamsLine(" [--localSize <size=0.04>] : In percentage");
553  addParamsLine(" [--useCriticalPoints <n=0>] : Use critical points instead of a grid");
554  addParamsLine(" : n is the number of critical points to choose");
555  addParamsLine(" : in each image");
556  addParamsLine(" [--threshold+ <th=-1>] : Correlation threshold");
557  addParamsLine(" [--blindSeqLength+ <n=-1>] : Blind sequence length, -1=No blind landmarks");
558  addParamsLine(" [--maxStep+ <step=4>] : Maximum step for chain refinement");
559  addParamsLine(" [--gridSamples+ <n=40>] : Total number of samples=n*n");
560  addParamsLine(" [--isCapillar+] : Set this flag if the tilt series is of a capillar");
561  addParamsLine(" == Step 2 (Determination of alignment parameters) Options == ");
562  addParamsLine(" [--psiMax+ <psi=-1>] : Maximum psi in absolute value (degrees)");
563  addParamsLine(" : -1 -> do not optimize for psi");
564  addParamsLine(" [--deltaRot+ <rot=5>] : In degrees. For the first optimization stage");
565  addParamsLine(" [--optimizeTiltAngle+] : Optimize tilt angle");
566  addParamsLine(" == Step 3 (Produce aligned images) Options == ");
567  addParamsLine(" [--dontNormalize+] : Don't normalize the output images");
568  addExampleLine("Typical run",false);
569  addExampleLine("xmipp_tomo_align_tilt_series -i tiltseries.sel --thr 8");
570  addExampleLine("If there are image with large shifts",false);
571  addExampleLine("xmipp_tomo_align_tilt_series -i tiltseries.sel --thr 8 --globalAffine");
572  addExampleLine("If there are clear landmarks that can be tracked",false);
573  addExampleLine("xmipp_tomo_align_tilt_series -i tiltseries.sel --thr 8 --criticalPoints");
574 }
575 
576 /* Produce side info ------------------------------------------------------- */
577 static pthread_mutex_t printingMutex = PTHREAD_MUTEX_INITIALIZER;
578 struct ThreadComputeTransformParams
579 {
580  int myThreadID;
581  ProgTomographAlignment * parent;
582 };
583 
584 void * threadComputeTransform( void * args )
585 {
586  auto * master =
587  (ThreadComputeTransformParams *) args;
588 
589  ProgTomographAlignment * parent = master->parent;
590  int thread_id = master->myThreadID;
591  int localnumThreads = parent->numThreads;
592  bool isCapillar = parent->isCapillar;
593  int Nimg = parent->Nimg;
594  double maxShiftPercentage = parent->maxShiftPercentage;
595  int maxIterDE = parent->maxIterDE;
596  bool showAffine = parent->showAffine;
597  double thresholdAffine = parent->thresholdAffine;
598  std::vector < MultidimArray<unsigned char> *> & img = parent->img;
599  std::vector< std::vector< Matrix2D<double> > > & affineTransformations = parent->affineTransformations;
600  double globalAffine = parent->globalAffine;
601 
602  int maxShift=FLOOR(XSIZE(*img[0])*maxShiftPercentage);
603  int initjj=1;
604  if (isCapillar)
605  initjj=0;
606 
607  initjj += thread_id;
608 
609  double cost;
610  for (int jj=initjj; jj<Nimg; jj+= localnumThreads)
611  {
612  int jj_1;
613  if (isCapillar)
614  jj_1=intWRAP(jj-1,0,Nimg-1);
615  else
616  jj_1=jj-1;
617  MultidimArray<unsigned char>& img_i=*img[jj_1];
618  MultidimArray<unsigned char>& img_j=*img[jj];
619  bool isMirror=(jj==0) && (jj_1==Nimg-1);
620 
621  if (MAT_XSIZE(affineTransformations[jj_1][jj])==0)
622  {
623  Matrix2D<double> Aij, Aji;
624  cost = computeAffineTransformation(img_i, img_j, maxShift,
625  maxIterDE, Aij, Aji,
626  showAffine, thresholdAffine, globalAffine,
627  isMirror,true, parent->pyramidLevel);
628  parent->correlationList[jj]=1-cost;
629 
630  pthread_mutex_lock( &printingMutex );
631  affineTransformations[jj_1][jj]=Aij;
632  affineTransformations[jj][jj_1]=Aji;
633  if (cost<1)
634  parent->writeTransformations(
635  parent->fnRoot+"_transformations.txt");
636  std::cout << "Cost for [" << jj_1 << "] - ["
637  << jj << "] = " << cost << std::endl;
638  parent->iteration++;
639  pthread_mutex_unlock( &printingMutex );
640  }
641  else
642  {
643  Matrix2D<double> Aij;
644  Aij=affineTransformations[jj_1][jj];
645 
646  MultidimArray<double> img_id, img_jd;
647  typeCast(img_i, img_id);
648  typeCast(img_j, img_jd);
649 
650  AffineFitness fitness;
651  setupAffineFitness(fitness, img_id, img_jd, maxShift, isMirror,
652  false, parent->pyramidLevel);
653  Matrix1D<double> p(6);
654  p(0)=Aij(0,0);
655  p(1)=Aij(0,1);
656  p(4)=Aij(0,2);
657  p(2)=Aij(1,0);
658  p(3)=Aij(1,1);
659  p(5)=Aij(1,2);
660  cost = fitness.affine_fitness_individual(MATRIX1D_ARRAY(p));
661  parent->correlationList[jj]=1-cost;
662 
663  pthread_mutex_lock( &printingMutex );
664  std::cout << "Cost for [" << jj_1 << "] - ["
665  << jj << "] = " << cost << std::endl;
666  pthread_mutex_unlock( &printingMutex );
667  parent->iteration++;
668  }
669  }
670 
671  return nullptr;
672 }
673 
674 void ProgTomographAlignment::computeAffineTransformations(
675  bool globalAffineToUse)
676 {
677  bool oldglobalAffine=globalAffine;
678  globalAffine=globalAffineToUse;
679 
680  auto * th_ids = new pthread_t[numThreads];
681  auto * th_args = new ThreadComputeTransformParams[numThreads];
682 
683  for( int nt = 0 ; nt < numThreads ; nt ++ )
684  {
685  // Passing parameters to each thread
686  th_args[nt].parent = this;
687  th_args[nt].myThreadID = nt;
688  pthread_create( (th_ids+nt) , nullptr, threadComputeTransform, (void *)(th_args+nt) );
689  }
690 
691  // Waiting for threads to finish
692  for( int nt = 0 ; nt < numThreads ; nt ++ )
693  pthread_join(*(th_ids+nt), nullptr);
694 
695  // Threads structures are not needed any more
696  delete[] th_ids;
697  delete[] th_args;
698  globalAffine=oldglobalAffine;
699 }
700 
701 void ProgTomographAlignment::identifyOutliers(bool mark)
702 {
703  isOutlier.initZeros(Nimg);
704  MultidimArray<double> correlationListAux(Nimg);
705  for (int i=0; i<Nimg; i++)
706  correlationListAux(i)=correlationList[i];
707 
708  MultidimArray<double> diff;
709  double medianCorr=correlationListAux.computeMedian();
710  diff=correlationListAux-medianCorr;
711  diff.selfABS();
712  double madCorr=diff.computeMedian();
713 
714  std::cout << "Cost distribution= " << 1-medianCorr << " +- "
715  << madCorr << std::endl;
716 
717  double thresholdCorr=medianCorr-identifyOutliersZ*madCorr;
718  for (size_t i=1; i<VEC_XSIZE(isOutlier); i++)
719  {
720  bool potentialOutlier=(correlationListAux(i)<thresholdCorr);
721  if (potentialOutlier)
722  {
723  affineTransformations[i-1][i].clear();
724  if (mark)
725  {
726  isOutlier(i)=true;
727  std::cout << name_list[i-1] << " [" << i
728  << "] is considered as an outlier. "
729  << "Its cost is " << 1-correlationListAux(i)
730  << std::endl;
731  }
732  else
733  std::cout << name_list[i-1] << " [" << i
734  << "] might be an outlier. "
735  << "Its cost is " << 1-correlationListAux(i)
736  << std::endl;
737  }
738  iteration++;
739  }
740 }
741 
742 void ProgTomographAlignment::produceSideInfo()
743 {
744  // Difficult images?
745  if (difficult)
746  {
747  if (pyramidLevel==0)
748  pyramidLevel=2;
749  if (localSize<0.05)
750  localSize=0.08;
751  if (seqLength==5)
752  seqLength=11;
753  }
754 
755  bestPreviousAlignment=new Alignment(this);
756  // Read input data
757  SF.read(fnSel,nullptr);
758  if (SF.containsLabel(MDL_ENABLED))
759  SF.removeObjects(MDValueEQ(MDL_ENABLED, -1));
760  Nimg=SF.size();
761  if (Nimg!=0)
762  {
763  // Clear the list of images if not empty
764  if (!img.empty())
765  {
766  for (size_t i=0; i<img.size(); i++)
767  delete img[i];
768  img.clear();
769  }
770 
771  // Clear the list of masks if not empty
772  if (!maskImg.empty())
773  {
774  for (size_t i=0; i<maskImg.size(); i++)
775  delete maskImg[i];
776  maskImg.clear();
777  }
778 
779  std::cerr << "Reading input data\n";
780  init_progress_bar(Nimg);
781  iteration=0;
782  int n=0;
783 
784  iMinTilt=-1;
785  double minTilt=1000;
786  bool nonZeroTilt=false;
787  Image<double> imgaux;
788  FileName fn;
789  for (size_t objId : SF.ids())
790  {
791  SF.getValue( MDL_IMAGE, fn, objId);
792  imgaux.read(fn);
793  if (difficult)
794  {
795  //ImageXmipp save;
796  //save()=imgaux(); save.write("PPPoriginal.xmp");
797  MultidimArray<double> Ifiltered;
798  Ifiltered=imgaux();
799  Ifiltered.setXmippOrigin();
800 
801  // Reject outliers
802  reject_outliers(Ifiltered,0.5);
803 
804  // Subtract background plane
805  substractBackgroundPlane(Ifiltered);
806 
807  // Subtract the background (rolling ball)
808  substractBackgroundRollingBall(Ifiltered,
809  XSIZE(Ifiltered)/10);
810 
811  // Bandpass the image
812  FourierFilter FilterBP;
813  FilterBP.FilterBand=BANDPASS;
814  FilterBP.w1=1.0/XSIZE(Ifiltered);
815  FilterBP.w2=100.0/XSIZE(Ifiltered);
816  FilterBP.raised_w=1.0/XSIZE(Ifiltered);
817  FilterBP.generateMask(Ifiltered);
818  FilterBP.applyMaskSpace(Ifiltered);
819 
820  // Equalize histogram
821  //histogram_equalization(Ifiltered,8);
822  imgaux()=Ifiltered;
823  //save()=Ifiltered; save.write("PPPpreprocessed.xmp");
824  //std::cout << "Press\n";
825  //char c; std::cin >> c;
826  }
827 
828  if (!useCriticalPoints)
829  {
830  auto* mask_i=new MultidimArray<unsigned char>;
831  generateMask(imgaux(),*mask_i,
832  XMIPP_MAX(ROUND(localSize*XSIZE(imgaux()))/2,5));
833  maskImg.push_back(mask_i);
834  }
835 
836  auto* img_i=new MultidimArray<unsigned char>;
837  imgaux().rangeAdjust(0,255);
838  typeCast(imgaux(),*img_i);
839  img_i->setXmippOrigin();
840  img.push_back(img_i);
841 
842  double tiltAngle;
843  if (SF.containsLabel(MDL_ANGLE_TILT))
844  SF.getValue(MDL_ANGLE_TILT,tiltAngle, objId);
845  else
846  tiltAngle=imgaux.tilt();
847  tiltList.push_back(tiltAngle);
848  if (tiltAngle!=0)
849  nonZeroTilt=true;
850  if (ABS(tiltAngle)<minTilt)
851  {
852  iMinTilt=n;
853  minTilt=ABS(tiltAngle);
854  }
855  name_list.push_back(imgaux.name());
856 
857  progress_bar(n++);
858  iteration++;
859  }
860  progress_bar(Nimg);
861  if (!nonZeroTilt)
862  REPORT_ERROR(ERR_VALUE_NOTSET,"Tilt angles have not been assigned to the input selfile");
863  }
864 
865  // Read images at original scale
866  if (!fnSelOrig.empty())
867  {
868  SForig.read(fnSelOrig,nullptr);
869  SForig.removeDisabled();
870 
871  if (SForig.size()!=SF.size())
872  REPORT_ERROR(ERR_MD_OBJECTNUMBER,"The number of images in both selfiles (-i and -iorig) is different");
873  }
874 
875  // Fill the affine transformations with empty matrices
876  std::vector< Matrix2D<double> > emptyRow;
877  for (int i=0; i<Nimg; i++)
878  {
879  Matrix2D<double> A;
880  emptyRow.push_back(A);
881  }
882  for (int i=0; i<Nimg; i++)
883  {
884  affineTransformations.push_back(emptyRow);
885  correlationList.push_back(-1);
886  }
887 
888  // Check if there are transformations already calculated
889  FileName fn_tmp = fnRoot+"_transformations.txt";
890  if (fn_tmp.exists())
891  {
892  std::ifstream fhIn;
893  fhIn.open(fn_tmp.c_str());
894  if (!fhIn)
895  REPORT_ERROR(ERR_IO_NOTEXIST,(String)"Cannot open "+ fn_tmp);
896  size_t linesRead=0;
897  while (!fhIn.eof())
898  {
899  std::string line;
900  getline(fhIn, line);
901  if (linesRead<affineTransformations.size()-1 && line!="")
902  {
903  std::vector< double > data;
904  readFloatList(line.c_str(), 8, data);
905  double tilt=data[0];
906  int i=0;
907  double bestDistance=ABS(tilt-tiltList[0]);
908  for (int j=0; j<Nimg; j++)
909  {
910  double distance=ABS(tilt-tiltList[j]);
911  if (bestDistance>distance)
912  {
913  bestDistance=distance;
914  i=j;
915  }
916  }
917  affineTransformations[i][i+1].resize(3,3);
918  affineTransformations[i][i+1].initIdentity(3);
919  affineTransformations[i][i+1](0,0)=data[2];
920  affineTransformations[i][i+1](1,0)=data[3];
921  affineTransformations[i][i+1](0,1)=data[4];
922  affineTransformations[i][i+1](1,1)=data[5];
923  affineTransformations[i][i+1](0,2)=data[6];
924  affineTransformations[i][i+1](1,2)=data[7];
925 
926  affineTransformations[i+1][i]=
927  affineTransformations[i][i+1].inv();
928 
929  if (showAffine)
930  {
931  MultidimArray<unsigned char>& img_i=*img[i];
932  MultidimArray<unsigned char>& img_j=*img[i+1];
933  int maxShift=FLOOR(XSIZE(*img[0])*maxShiftPercentage);
934  Matrix2D<double> Aij, Aji;
935  Aij=affineTransformations[i][i+1];
936  Aji=affineTransformations[i+1][i];
937  bool isMirror=false;
938  std::cout << "Transformation between " << i << " ("
939  << tiltList[i] << ") and " << i+1 << " ("
940  << tiltList[i+1] << std::endl;
941  double cost = computeAffineTransformation(img_i, img_j,
942  maxShift, maxIterDE, Aij, Aji,
943  showAffine, thresholdAffine, globalAffine,
944  isMirror,true, pyramidLevel);
945  if (cost<0)
946  {
947  affineTransformations[i][i+1].clear();
948  affineTransformations[i+1][i].clear();
949  writeTransformations(fn_tmp);
950  }
951  }
952  }
953  linesRead++;
954  }
955  fhIn.close();
956  }
957 
958  // Compute the rest of transformations
959  computeAffineTransformations(globalAffine);
960 
961  // Do not show refinement
962  showRefinement = false;
963 
964  // Check which is the distribution of correlation
965  if (!useCriticalPoints)
966  {
967  auto X0=(int)(STARTINGX(*(img[0]))+2*localSize*XSIZE(*(img[0])));
968  auto XF=(int)(FINISHINGX(*(img[0]))-2*localSize*XSIZE(*(img[0])));
969  auto Y0=(int)(STARTINGY(*(img[0]))+2*localSize*YSIZE(*(img[0])));
970  auto YF=(int)(FINISHINGY(*(img[0]))-2*localSize*YSIZE(*(img[0])));
971  avgForwardPatchCorr.initZeros(Nimg);
972  avgBackwardPatchCorr.initZeros(Nimg);
973  avgForwardPatchCorr.initConstant(1);
974  avgBackwardPatchCorr.initConstant(1);
975 
976  for (int ii=0; ii<Nimg; ++ii)
977  {
978  // Compute average forward patch corr
979  MultidimArray<double> corrList;
980  if (ii<=Nimg-2)
981  {
982  corrList.initZeros(200);
983  FOR_ALL_ELEMENTS_IN_ARRAY1D(corrList)
984  {
985  Matrix1D<double> rii(3), rjj(3);
986  do
987  {
988  XX(rii)=rnd_unif(X0,XF);
989  YY(rii)=rnd_unif(Y0,YF);
990  rjj=affineTransformations[ii][ii+1]*rii;
991  refineLandmark(ii,ii+1,rii,rjj,corrList(i),false);
992  }
993  while (corrList(i)<-0.99);
994  }
995  avgForwardPatchCorr(ii)=corrList.computeAvg();
996  }
997 
998  // Compute average backward patch corr
999  if (ii>0)
1000  {
1001  corrList.initZeros(200);
1002  FOR_ALL_ELEMENTS_IN_ARRAY1D(corrList)
1003  {
1004  Matrix1D<double> rii(3), rjj(3);
1005  do
1006  {
1007  XX(rii)=rnd_unif(X0,XF);
1008  YY(rii)=rnd_unif(Y0,YF);
1009  rjj=affineTransformations[ii][ii-1]*rii;
1010  refineLandmark(ii,ii-1,rii,rjj,corrList(i),false);
1011  }
1012  while (corrList(i)<-0.99);
1013  }
1014  avgBackwardPatchCorr(ii)=corrList.computeAvg();
1015  }
1016 
1017  std::cout << "Image " << ii << " Average correlation forward="
1018  << avgForwardPatchCorr(ii)
1019  << " backward=" << avgBackwardPatchCorr(ii) << std::endl;
1020  iteration++;
1021  }
1022  }
1023 
1024  // Identify outliers
1025  if (!doNotIdentifyOutliers)
1026  {
1027  identifyOutliers(false);
1028  computeAffineTransformations(false);
1029  identifyOutliers(true);
1030  }
1031  else
1032  isOutlier.initZeros(Nimg);
1033  if (lastStep==0)
1034  exit(0);
1035 }
1036 
1037 /* Generate landmark set --------------------------------------------------- */
1038 //#define DEBUG
1039 struct ThreadGenerateLandmarkSetParams
1040 {
1041  int myThreadID;
1042  ProgTomographAlignment * parent;
1043 
1044  std::vector<LandmarkChain> *chainList;
1045 };
1046 
1047 void * threadgenerateLandmarkSetGrid( void * args )
1048 {
1049  auto * master =
1050  (ThreadGenerateLandmarkSetParams *) args;
1051  ProgTomographAlignment * parent = master->parent;
1052  int thread_id = master->myThreadID;
1053  int Nimg=parent->Nimg;
1054  int numThreads=parent->numThreads;
1055  const std::vector< std::vector< Matrix2D<double> > > &affineTransformations=
1056  parent->affineTransformations;
1057  int gridSamples=parent->gridSamples;
1058 
1059  auto deltaShift=(int)floor(XSIZE(*(parent->img)[0])/gridSamples);
1060  master->chainList=new std::vector<LandmarkChain>;
1061  Matrix1D<double> rii(3), rjj(3);
1062  ZZ(rii)=1;
1063  ZZ(rjj)=1;
1064  if (thread_id==0)
1065  init_progress_bar(gridSamples);
1066  int includedPoints=0;
1067  Matrix1D<int> visited(Nimg);
1068  for (int nx=thread_id; nx<gridSamples; nx+=numThreads)
1069  {
1070  XX(rii)=STARTINGX(*(parent->img)[0])+ROUND(deltaShift*(0.5+nx));
1071  for (int ny=0; ny<gridSamples; ny+=1)
1072  {
1073  YY(rii)=STARTINGY(*(parent->img)[0])+ROUND(deltaShift*(0.5+ny));
1074  for (int ii=0; ii<=Nimg-1; ++ii)
1075  {
1076  // Check if inside the mask
1077  if (!(*(parent->maskImg[ii]))((int)YY(rii),(int)XX(rii)))
1078  continue;
1079  if (parent->isOutlier(ii))
1080  continue;
1081 
1082  LandmarkChain chain;
1083  chain.clear();
1084  Landmark l;
1085  l.x=XX(rii);
1086  l.y=YY(rii);
1087  l.imgIdx=ii;
1088  chain.push_back(l);
1089  visited.initZeros();
1090  visited(ii)=1;
1091 
1092  // Follow this landmark backwards
1093  bool acceptLandmark=true;
1094  int jj;
1095  if (parent->isCapillar)
1096  jj=intWRAP(ii-1,0,Nimg-1);
1097  else
1098  jj=ii-1;
1099  Matrix2D<double> Aij, Aji;
1100  Matrix1D<double> rcurrent=rii;
1101  while (jj>=0 && acceptLandmark && !visited(jj) &&
1102  !parent->isOutlier(jj))
1103  {
1104  // Compute the affine transformation between ii and jj
1105  int jj_1;
1106  if (parent->isCapillar)
1107  jj_1=intWRAP(jj+1,0,Nimg-1);
1108  else
1109  jj_1=jj+1;
1110  visited(jj)=1;
1111  Aij=affineTransformations[jj][jj_1];
1112  Aji=affineTransformations[jj_1][jj];
1113  rjj=Aji*rcurrent;
1114  double corr;
1115  acceptLandmark=parent->refineLandmark(jj_1,jj,rcurrent,rjj,
1116  corr,true);
1117  if (acceptLandmark)
1118  {
1119  l.x=XX(rjj);
1120  l.y=YY(rjj);
1121  l.imgIdx=jj;
1122  chain.push_back(l);
1123  if (parent->isCapillar)
1124  jj=intWRAP(jj-1,0,Nimg-1);
1125  else
1126  jj=jj-1;
1127  rcurrent=rjj;
1128  }
1129  }
1130 
1131  // Follow this landmark forward
1132  acceptLandmark=true;
1133  if (parent->isCapillar)
1134  jj=intWRAP(ii+1,0,Nimg-1);
1135  else
1136  jj=ii+1;
1137  rcurrent=rii;
1138  while (jj<Nimg && acceptLandmark && !visited(jj) &&
1139  !parent->isOutlier(jj))
1140  {
1141  // Compute the affine transformation between ii and jj
1142  int jj_1;
1143  if (parent->isCapillar)
1144  jj_1=intWRAP(jj-1,0,Nimg-1);
1145  else
1146  jj_1=jj-1;
1147  visited(jj)=1;
1148  Aij=affineTransformations[jj_1][jj];
1149  Aji=affineTransformations[jj][jj_1];
1150  rjj=Aij*rcurrent;
1151  double corr;
1152  acceptLandmark=parent->refineLandmark(jj_1,jj,rcurrent,rjj,
1153  corr,true);
1154  if (acceptLandmark)
1155  {
1156  l.x=XX(rjj);
1157  l.y=YY(rjj);
1158  l.imgIdx=jj;
1159  chain.push_back(l);
1160  if (parent->isCapillar)
1161  jj=intWRAP(jj+1,0,Nimg-1);
1162  else
1163  jj=jj+1;
1164  rcurrent=rjj;
1165  }
1166  }
1167 
1168 #ifdef DEBUG
1169  std::cout << "img=" << ii << " chain length="
1170  << chain.size() << " [" << jjleft
1171  << " - " << jjright << "]";
1172 #endif
1173 
1174  if (chain.size()>parent->seqLength)
1175  {
1176  double corrChain;
1177  bool accepted=parent->refineChain(chain,corrChain);
1178  if (accepted)
1179  {
1180 #ifdef DEBUG
1181  std::cout << " Accepted with length= "
1182  << chain.size()
1183  << " [" << chain[0].imgIdx << " - "
1184  << chain[chain.size()-1].imgIdx << "]= ";
1185  for (int i=0; i<chain.size(); i++)
1186  std::cout << chain[i].imgIdx << " ";
1187 #endif
1188 
1189  master->chainList->push_back(chain);
1190  includedPoints+=chain.size();
1191  }
1192  }
1193 #ifdef DEBUG
1194  std::cout << std::endl;
1195 #endif
1196 
1197  }
1198 #ifdef DEBUG
1199  std::cout << "Point nx=" << nx << " ny=" << ny
1200  << " Number of points="
1201  << includedPoints
1202  << " Number of chains=" << master->chainList->size()
1203  << " ( " << ((double) includedPoints)/
1204  master->chainList->size() << " )\n";
1205 #endif
1206 
1207  }
1208  if (thread_id==0)
1209  progress_bar(nx);
1210  }
1211  if (thread_id==0)
1212  progress_bar(gridSamples);
1213  return nullptr;
1214 }
1215 
1216 void * threadgenerateLandmarkSetBlind( void * args )
1217 {
1218  auto * master =
1219  (ThreadGenerateLandmarkSetParams *) args;
1220  ProgTomographAlignment * parent = master->parent;
1221  int thread_id = master->myThreadID;
1222  int Nimg=parent->Nimg;
1223  int numThreads=parent->numThreads;
1224  const std::vector< std::vector< Matrix2D<double> > > &affineTransformations=
1225  parent->affineTransformations;
1226  int gridSamples=parent->gridSamples;
1227 
1228  auto deltaShift=(int)floor(XSIZE(*(parent->img)[0])/gridSamples);
1229  master->chainList=new std::vector<LandmarkChain>;
1230  Matrix1D<double> rii(3), rjj(3);
1231  ZZ(rii)=1;
1232  ZZ(rjj)=1;
1233  if (thread_id==0)
1234  init_progress_bar(gridSamples);
1235  int includedPoints=0;
1236  int maxSideLength=(parent->blindSeqLength-1)/2;
1237  for (int nx=thread_id; nx<gridSamples; nx+=numThreads)
1238  {
1239  XX(rii)=STARTINGX(*(parent->img)[0])+ROUND(deltaShift*(0.5+nx));
1240  for (int ny=0; ny<gridSamples; ny+=1)
1241  {
1242  YY(rii)=STARTINGY(*(parent->img)[0])+ROUND(deltaShift*(0.5+ny));
1243  for (int ii=0; ii<=Nimg-1; ++ii)
1244  {
1245  // Check if inside the mask
1246  if (!(*(parent->maskImg[ii]))((int)YY(rii),(int)XX(rii)))
1247  continue;
1248  if (parent->isOutlier(ii))
1249  continue;
1250 
1251  LandmarkChain chain;
1252  chain.clear();
1253  Landmark l;
1254  l.x=XX(rii);
1255  l.y=YY(rii);
1256  l.imgIdx=ii;
1257  chain.push_back(l);
1258 
1259  // Follow this landmark backwards
1260  bool acceptLandmark=true;
1261  int jj;
1262  int sideLength=0;
1263  if (parent->isCapillar)
1264  jj=intWRAP(ii-1,0,Nimg-1);
1265  else
1266  jj=ii-1;
1267  Matrix2D<double> Aij, Aji;
1268  Matrix1D<double> rcurrent=rii;
1269  while (jj>=0 && acceptLandmark && sideLength<maxSideLength &&
1270  !parent->isOutlier(jj))
1271  {
1272  // Compute the affine transformation between ii and jj
1273  int jj_1;
1274  if (parent->isCapillar)
1275  jj_1=intWRAP(jj+1,0,Nimg-1);
1276  else
1277  jj_1=jj+1;
1278  Aij=affineTransformations[jj][jj_1];
1279  Aji=affineTransformations[jj_1][jj];
1280  rjj=Aji*rcurrent;
1281  auto iYYrjj=(int)YY(rjj);
1282  auto iXXrjj=(int)XX(rjj);
1283  if (!(*(parent->maskImg[jj])).outside(iYYrjj,iXXrjj))
1284  acceptLandmark=(*(parent->maskImg[jj]))(iYYrjj,iXXrjj);
1285  else
1286  acceptLandmark=false;
1287  if (acceptLandmark)
1288  {
1289  l.x=XX(rjj);
1290  l.y=YY(rjj);
1291  l.imgIdx=jj;
1292  chain.push_back(l);
1293  if (parent->isCapillar)
1294  jj=intWRAP(jj-1,0,Nimg-1);
1295  else
1296  jj=jj-1;
1297  rcurrent=rjj;
1298  sideLength++;
1299  }
1300  }
1301 
1302  // Follow this landmark forward
1303  acceptLandmark=true;
1304  if (parent->isCapillar)
1305  jj=intWRAP(ii+1,0,Nimg-1);
1306  else
1307  jj=ii+1;
1308  rcurrent=rii;
1309  sideLength=0;
1310  while (jj<Nimg && acceptLandmark && sideLength<maxSideLength &&
1311  !parent->isOutlier(jj))
1312  {
1313  // Compute the affine transformation between ii and jj
1314  int jj_1;
1315  if (parent->isCapillar)
1316  jj_1=intWRAP(jj-1,0,Nimg-1);
1317  else
1318  jj_1=jj-1;
1319  Aij=affineTransformations[jj_1][jj];
1320  Aji=affineTransformations[jj][jj_1];
1321  rjj=Aij*rcurrent;
1322  auto iYYrjj=(int)YY(rjj);
1323  auto iXXrjj=(int)XX(rjj);
1324  if (!(*(parent->maskImg[jj])).outside(iYYrjj,iXXrjj))
1325  acceptLandmark=(*(parent->maskImg[jj]))(iYYrjj,iXXrjj);
1326  else
1327  acceptLandmark=false;
1328  if (acceptLandmark)
1329  {
1330  l.x=XX(rjj);
1331  l.y=YY(rjj);
1332  l.imgIdx=jj;
1333  chain.push_back(l);
1334  if (parent->isCapillar)
1335  jj=intWRAP(jj+1,0,Nimg-1);
1336  else
1337  jj=jj+1;
1338  rcurrent=rjj;
1339  sideLength++;
1340  }
1341  }
1342 
1343 #ifdef DEBUG
1344  std::cout << "img=" << ii << " chain length="
1345  << chain.size() << " [" << jjleft
1346  << " - " << jjright << "]\n";
1347 #endif
1348 
1349  master->chainList->push_back(chain);
1350  includedPoints+=chain.size();
1351  }
1352 #ifdef DEBUG
1353  std::cout << "Point nx=" << nx << " ny=" << ny
1354  << " Number of points="
1355  << includedPoints
1356  << " Number of chains=" << master->chainList->size()
1357  << " ( " << ((double) includedPoints)/
1358  master->chainList->size() << " )\n";
1359 #endif
1360 
1361  }
1362  if (thread_id==0)
1363  progress_bar(nx);
1364  }
1365  if (thread_id==0)
1366  progress_bar(gridSamples);
1367  return nullptr;
1368 }
1369 
1370 //#define DEBUG
1371 void * threadgenerateLandmarkSetCriticalPoints( void * args )
1372 {
1373  auto * master =
1374  (ThreadGenerateLandmarkSetParams *) args;
1375  ProgTomographAlignment * parent = master->parent;
1376  int thread_id = master->myThreadID;
1377  int Nimg=parent->Nimg;
1378  int numThreads=parent->numThreads;
1379  const std::vector< std::vector< Matrix2D<double> > > &affineTransformations=
1380  parent->affineTransformations;
1381 
1382  master->chainList=new std::vector<LandmarkChain>;
1383  std::vector<LandmarkChain> candidateChainList;
1384  if (thread_id==0)
1385  init_progress_bar(Nimg);
1386  int halfSeqLength=parent->seqLength/2;
1387 
1388  // Design a mask for the dilation
1389  int radius=4;
1390  MultidimArray<int> mask;
1391  mask.resize(2*radius+1,2*radius+1);
1392  mask.setXmippOrigin();
1393  BinaryCircularMask(mask,4,OUTSIDE_MASK);
1394 
1395  Image<double> I;
1396  for (int ii=thread_id; ii<=Nimg-1; ii+=numThreads)
1397  {
1398  if (parent->isOutlier(ii))
1399  continue;
1400  I.read(parent->name_list[ii]);
1401 
1402  // Generate mask
1403  MultidimArray<unsigned char> largeMask;
1404  generateMask(I(),largeMask,
1405  XMIPP_MAX(ROUND(parent->localSize*XSIZE(I()))/2,5));
1406 
1407  // Filter the image
1408  MultidimArray<double> Ifiltered;
1409  Ifiltered=I();
1410  Ifiltered.setXmippOrigin();
1411  FourierFilter FilterBP;
1412  FilterBP.FilterBand=BANDPASS;
1413  FilterBP.w1=2.0/XSIZE(Ifiltered);
1414  FilterBP.w2=128.0/XSIZE(Ifiltered);
1415  FilterBP.raised_w=1.0/XSIZE(Ifiltered);
1416  ;
1417  FilterBP.generateMask(Ifiltered);
1418  FilterBP.applyMaskSpace(Ifiltered);
1419 
1420  // Identify low valued points and perform dilation
1421  MultidimArray<double> Iaux=Ifiltered;
1422  Iaux.selfWindow(
1423  -ROUND(0.45*YSIZE(Ifiltered)),-ROUND(0.45*XSIZE(Ifiltered)),
1424  ROUND(0.45*YSIZE(Ifiltered)), ROUND(0.45*XSIZE(Ifiltered)));
1425  Histogram1D hist;
1426  compute_hist(Iaux, hist, 400);
1427  double th=hist.percentil(2);
1428  std::vector< Matrix1D<double> > Q;
1429  FOR_ALL_ELEMENTS_IN_ARRAY2D(Iaux)
1430  if (Ifiltered(i,j)<th && largeMask(i,j))
1431  {
1432  // Check if it is a local minimum
1433  bool localMinimum=true;
1434  int x0=j;
1435  int y0=i;
1436  FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
1437  {
1438  int x=x0+j;
1439  int y=y0+i;
1440  if (x>=STARTINGX(Iaux) && x<=FINISHINGX(Iaux) &&
1441  y>=STARTINGY(Iaux) && y<=FINISHINGY(Iaux))
1442  if (A2D_ELEM(Iaux,y,x)<A2D_ELEM(Iaux,y0,x0))
1443  {
1444  localMinimum=false;
1445  break;
1446  }
1447  }
1448  if (localMinimum)
1449  Q.push_back(vectorR2(x0,y0));
1450  }
1451 
1452  // Check that the list is not too long
1453  if (Q.size()>10*parent->Ncritical)
1454  {
1455  size_t qmax=Q.size();
1456  MultidimArray<double> minDistance;
1457  minDistance.resize(qmax);
1458  minDistance.initConstant(1e20);
1459  for (size_t q1=0; q1<qmax; q1++)
1460  for (size_t q2=q1+1; q2< qmax; q2++)
1461  {
1462  double diffX=XX(Q[q1])-XX(Q[q2]);
1463  double diffY=YY(Q[q1])-YY(Q[q2]);
1464  double d12=diffX*diffX+diffY*diffY;
1465  if (d12<minDistance(q1))
1466  minDistance(q1)=d12;
1467  if (d12<minDistance(q2))
1468  minDistance(q2)=d12;
1469  }
1470  MultidimArray<int> idxDistanceSort;
1471  minDistance.indexSort(idxDistanceSort);
1472  std::vector< Matrix1D<double> > Qaux;
1473  int qlimit=XMIPP_MIN(10*(parent->Ncritical),qmax);
1474  for (int q=0; q<qlimit; q++)
1475  Qaux.push_back(Q[idxDistanceSort(qmax-1-q)-1]);
1476  Q.clear();
1477  Q=Qaux;
1478  }
1479 
1480  int qmax=Q.size();
1481  Matrix1D<double> rii(3), rjj(3);
1482  ZZ(rii)=1;
1483  ZZ(rjj)=1;
1484  MultidimArray<double> corrQ;
1485  corrQ.initZeros(qmax);
1486  for (int q=0; q<qmax; q++)
1487  {
1488  XX(rii)=XX(Q[q]);
1489  YY(rii)=YY(Q[q]);
1490 
1491  // Initiate a chain here
1492  LandmarkChain chain;
1493  chain.clear();
1494  Landmark l;
1495  l.x=XX(rii);
1496  l.y=YY(rii);
1497  l.imgIdx=ii;
1498  chain.push_back(l);
1499 
1500  // Follow this landmark backwards
1501  int jjmin=XMIPP_MAX(0,ii-halfSeqLength);
1502  int jjmax=ii-1;
1503  Matrix2D<double> Aij, Aji;
1504  Matrix1D<double> rcurrent=rii;
1505  for (int jj=jjmax; jj>=jjmin; --jj)
1506  {
1507  if (parent->isOutlier(jj))
1508  break;
1509 
1510  // Compute the affine transformation between jj and jj+1
1511  int jj_1=jj+1;
1512  Aij=affineTransformations[jj][jj_1];
1513  Aji=affineTransformations[jj_1][jj];
1514  rjj=Aji*rcurrent;
1515  double corr;
1516  parent->refineLandmark(jj_1,jj,rcurrent,rjj, corr,true);
1517  l.x=XX(rjj);
1518  l.y=YY(rjj);
1519  l.imgIdx=jj;
1520  chain.push_back(l);
1521  rcurrent=rjj;
1522  }
1523 
1524  // Follow this landmark forwards
1525  jjmin=ii+1;
1526  jjmax=XMIPP_MIN(Nimg-1,ii+halfSeqLength);
1527  rcurrent=rii;
1528  for (int jj=jjmin; jj<=jjmax; ++jj)
1529  {
1530  if (parent->isOutlier(jj))
1531  break;
1532 
1533  // Compute the affine transformation between jj-1 and jj
1534  int jj_1=jj-1;
1535  Aij=affineTransformations[jj_1][jj];
1536  Aji=affineTransformations[jj][jj_1];
1537  rjj=Aij*rcurrent;
1538  double corr;
1539  parent->refineLandmark(jj_1,jj,rcurrent,rjj,corr,true);
1540  l.x=XX(rjj);
1541  l.y=YY(rjj);
1542  l.imgIdx=jj;
1543  chain.push_back(l);
1544  rcurrent=rjj;
1545  }
1546 
1547  // Refine chain
1548  parent->refineChain(chain,corrQ(q));
1549  candidateChainList.push_back(chain);
1550  }
1551  if (thread_id==0)
1552  progress_bar(ii);
1553 
1554  // Sort all chains according to its correlation
1555  MultidimArray<int> idx;
1556  corrQ.indexSort(idx);
1557  int imax=XMIPP_MIN(XSIZE(idx),parent->Ncritical);
1558  for (int iq=0; iq<imax; iq++)
1559  {
1560  int q=idx(XSIZE(idx)-1-iq)-1;
1561  if (corrQ(q)>0.5)
1562  {
1563  master->chainList->push_back(candidateChainList[q]);
1564 #ifdef DEBUG
1565 
1566  std::cout << "Corr " << iq << ": " << corrQ(q) << ":";
1567  for (int i=0; i<candidateChainList[q].size(); i++)
1568  std::cout << candidateChainList[q][i].imgIdx << " ";
1569  std::cout << std::endl;
1570 #endif
1571 
1572  }
1573  else
1574  std::cout << "It's recommended to reduce the number of critical points\n";
1575  }
1576  candidateChainList.clear();
1577 
1578 #ifdef DEBUG
1579 
1580  Image<double> save;
1581  typeCast(*(parent->img[ii]),save());
1582  save.write("PPPoriginal.xmp");
1583  save()=Ifiltered;
1584  save.write("PPPfiltered.xmp");
1585  double minval=Ifiltered.computeMin();
1586  FOR_ALL_ELEMENTS_IN_MATRIX2D(Ifiltered)
1587  if (Ifiltered(i,j)>=th || !largeMask(i,j))
1588  save(i,j)=th;
1589  for (int q=0; q<Q.size(); q++)
1590  FOR_ALL_ELEMENTS_IN_MATRIX2D(mask)
1591  if (YY(Q[q])+i>=STARTINGY(Ifiltered) && YY(Q[q])+i<=FINISHINGY(Ifiltered) &&
1592  XX(Q[q])+j>=STARTINGX(Ifiltered) && XX(Q[q])+j<=FINISHINGX(Ifiltered))
1593  save(YY(Q[q])+i,XX(Q[q])+j)=minval;
1594  imax=XMIPP_MIN(XSIZE(idx),parent->Ncritical);
1595  for (int iq=0; iq<imax; iq++)
1596  {
1597  int q=idx(XSIZE(idx)-1-iq)-1;
1598  FOR_ALL_ELEMENTS_IN_MATRIX2D(mask)
1599  if (YY(Q[q])+i>=STARTINGY(Ifiltered) && YY(Q[q])+i<=FINISHINGY(Ifiltered) &&
1600  XX(Q[q])+j>=STARTINGX(Ifiltered) && XX(Q[q])+j<=FINISHINGX(Ifiltered))
1601  save(YY(Q[q])+i,XX(Q[q])+j)=(minval+th)/2;
1602  }
1603 
1604  save.write("PPPcritical.xmp");
1605  std::cout << "Number of critical points=" << Q.size() << std::endl;
1606  std::cout << "CorrQ stats:";
1607  corrQ.printStats();
1608  std::cout << std::endl;
1609  std::cout << "Press any key\n";
1610  char c;
1611  std::cin >> c;
1612 #endif
1613 
1614  }
1615  if (thread_id==0)
1616  progress_bar(Nimg);
1617  return nullptr;
1618 }
1619 #undef DEBUG
1620 
1621 ProgTomographAlignment::~ProgTomographAlignment()
1622 {
1623  // Clear the list of images if not empty
1624  if (!img.empty())
1625  {
1626  for (size_t i=0; i<img.size(); i++)
1627  delete img[i];
1628  img.clear();
1629  }
1630 
1631  if (!maskImg.empty())
1632  {
1633  for (size_t i=0; i<maskImg.size(); i++)
1634  delete maskImg[i];
1635  maskImg.clear();
1636  }
1637 }
1638 
1639 void ProgTomographAlignment::generateLandmarkSet()
1640 {
1641  FileName fn_tmp = fnRoot+"_landmarks.txt";
1642  if (!fn_tmp.exists())
1643  {
1644  auto * th_ids = new pthread_t[numThreads];
1645  auto * th_args=
1646  new ThreadGenerateLandmarkSetParams[numThreads];
1647  for( int nt = 0 ; nt < numThreads ; nt ++ )
1648  {
1649  th_args[nt].parent = this;
1650  th_args[nt].myThreadID = nt;
1651  if (useCriticalPoints)
1652  pthread_create( (th_ids+nt) , nullptr, threadgenerateLandmarkSetCriticalPoints, (void *)(th_args+nt) );
1653  else
1654  pthread_create( (th_ids+nt) , nullptr, threadgenerateLandmarkSetGrid, (void *)(th_args+nt) );
1655  }
1656 
1657  std::vector<LandmarkChain> chainList;
1658  int includedPoints=0;
1659  for( int nt = 0 ; nt < numThreads ; nt ++ )
1660  {
1661  pthread_join(*(th_ids+nt), nullptr);
1662  int imax=th_args[nt].chainList->size();
1663  for (int i=0; i<imax; i++)
1664  {
1665  chainList.push_back((*th_args[nt].chainList)[i]);
1666  includedPoints+=(*th_args[nt].chainList)[i].size();
1667  }
1668  delete th_args[nt].chainList;
1669  }
1670  delete[] th_ids;
1671  delete[] th_args;
1672 
1673  // Add blind landmarks
1674  if (blindSeqLength>0)
1675  {
1676  th_ids = new pthread_t[numThreads];
1677  th_args = new ThreadGenerateLandmarkSetParams[numThreads];
1678  for( int nt = 0 ; nt < numThreads ; nt ++ )
1679  {
1680  th_args[nt].parent = this;
1681  th_args[nt].myThreadID = nt;
1682  pthread_create( (th_ids+nt) , nullptr, threadgenerateLandmarkSetBlind, (void *)(th_args+nt) );
1683  }
1684 
1685  for( int nt = 0 ; nt < numThreads ; nt ++ )
1686  {
1687  pthread_join(*(th_ids+nt), nullptr);
1688  int imax=th_args[nt].chainList->size();
1689  for (int i=0; i<imax; i++)
1690  {
1691  chainList.push_back((*th_args[nt].chainList)[i]);
1692  includedPoints+=(*th_args[nt].chainList)[i].size();
1693  }
1694  delete th_args[nt].chainList;
1695  }
1696  delete[] th_ids;
1697  delete[] th_args;
1698  }
1699 
1700  // Generate the landmark "matrix"
1701  allLandmarksX.resize(chainList.size(),Nimg);
1702  allLandmarksY.resize(chainList.size(),Nimg);
1703  allLandmarksX.initConstant(XSIZE(*img[0]));
1704  allLandmarksY.initConstant(YSIZE(*img[0]));
1705  for (size_t i=0; i<chainList.size(); i++)
1706  {
1707  for (size_t j=0; j<chainList[i].size(); j++)
1708  {
1709  int idx=chainList[i][j].imgIdx;
1710  allLandmarksX(i,idx)=chainList[i][j].x;
1711  allLandmarksY(i,idx)=chainList[i][j].y;
1712  }
1713  }
1714 
1715  // Write landmarks
1716  if (includedPoints>0)
1717  {
1718  writeLandmarkSet(fnRoot+"_landmarks.txt");
1719  std::cout << " Number of points="
1720  << includedPoints
1721  << " Number of chains=" << chainList.size()
1722  << " ( " << ((double) includedPoints)/chainList.size() << " )\n";
1723  }
1724  else
1725  REPORT_ERROR(ERR_VALUE_INCORRECT,"There are no landmarks meeting this threshold. Try to lower it");
1726  }
1727  else
1728  {
1729  readLandmarkSet(fnRoot+"_landmarks.txt");
1730  }
1731 }
1732 #undef DEBUG
1733 
1734 /* Refine landmark --------------------------------------------------------- */
1735 bool ProgTomographAlignment::refineLandmark(int ii, int jj,
1736  const Matrix1D<double> &rii, Matrix1D<double> &rjj, double &maxCorr,
1737  bool tryFourier) const
1738 {
1739  maxCorr=-1;
1740  int halfSize=XMIPP_MAX(ROUND(localSize*XSIZE(*img[ii]))/2,5);
1741  if (XX(rii)-halfSize<STARTINGX(*img[ii]) ||
1742  XX(rii)+halfSize>FINISHINGX(*img[ii]) ||
1743  YY(rii)-halfSize<STARTINGY(*img[ii]) ||
1744  YY(rii)+halfSize>FINISHINGY(*img[ii]) ||
1745  XX(rjj)-halfSize<STARTINGX(*img[jj]) ||
1746  XX(rjj)+halfSize>FINISHINGX(*img[jj]) ||
1747  YY(rjj)-halfSize<STARTINGY(*img[jj]) ||
1748  YY(rjj)+halfSize>FINISHINGY(*img[jj]))
1749  return 0;
1750 
1751  // Check if the two pieces are reversed
1752  bool reversed=isCapillar && ABS(ii-jj)>Nimg/2;
1753 
1754  // Select piece in image ii, compute its statistics and normalize
1755  MultidimArray<double> pieceii(2*halfSize+1,2*halfSize+1);
1756  pieceii.setXmippOrigin();
1757  const MultidimArray<unsigned char> &Iii=(*img[ii]);
1758  FOR_ALL_ELEMENTS_IN_ARRAY2D(pieceii)
1759  A2D_ELEM(pieceii,i,j)=A2D_ELEM(Iii,
1760  (int)(YY(rii)+i),(int)(XX(rii)+j));
1761  if (showRefinement)
1762  {
1763  Image<double> save;
1764  save()=pieceii;
1765  save.write("PPPpieceii.xmp");
1766  std::cout << "ii=" << ii << " jj=" << jj << std::endl;
1767  std::cout << "rii=" << rii.transpose() << std::endl;
1768  }
1769 
1770  // Choose threshold
1771  double actualCorrThreshold=corrThreshold;
1772  if (actualCorrThreshold<0 && VEC_XSIZE(avgForwardPatchCorr)>0)
1773  {
1774  if (ii<jj)
1775  actualCorrThreshold=XMIPP_MIN(avgForwardPatchCorr(ii),
1776  avgBackwardPatchCorr(jj));
1777  else
1778  actualCorrThreshold=XMIPP_MIN(avgBackwardPatchCorr(ii),
1779  avgForwardPatchCorr(jj));
1780  }
1781  if (showRefinement)
1782  std::cout << "actualCorrThreshold=" << actualCorrThreshold << std::endl;
1783 
1784  // Try Fourier
1785  if (tryFourier)
1786  {
1787  const MultidimArray<unsigned char> &Ijj=(*img[jj]);
1788  if (XX(rjj)-halfSize>=STARTINGX(Ijj) &&
1789  XX(rjj)+halfSize<=FINISHINGX(Ijj) &&
1790  YY(rjj)-halfSize>=STARTINGY(Ijj) &&
1791  YY(rjj)+halfSize<=FINISHINGY(Ijj))
1792  {
1793  // Take the piece at jj
1794  MultidimArray<double> piecejj(2*halfSize+1,2*halfSize+1);
1795  piecejj.setXmippOrigin();
1796  FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
1797  A2D_ELEM(piecejj,i,j)=A2D_ELEM(Ijj,
1798  (int)(YY(rjj)+i),(int)(XX(rjj)+j));
1799 
1800  // Compute its correlation with pieceii
1801  double corrOriginal=correlationIndex(pieceii,piecejj);
1802  if (showRefinement)
1803  {
1804  Image<double> save;
1805  save()=piecejj;
1806  save.write("PPPpiecejjOriginal.xmp");
1807  std::cout << "Corr original=" << corrOriginal << std::endl;
1808  }
1809 
1810  // Now try with the best shift
1811  double shiftX,shiftY;
1812  CorrelationAux aux;
1813  bestNonwrappingShift(pieceii,piecejj,shiftX,shiftY,aux);
1814  Matrix1D<double> fftShift(2);
1815  VECTOR_R2(fftShift,shiftX,shiftY);
1816  selfTranslate(xmipp_transformation::LINEAR,piecejj,fftShift,xmipp_transformation::WRAP);
1817  double corrFFT=correlationIndex(pieceii,piecejj);
1818  if (corrFFT>corrOriginal)
1819  {
1820  XX(rjj)-=shiftX;
1821  YY(rjj)-=shiftY;
1822  }
1823 
1824  if (showRefinement)
1825  {
1826  Image<double> save;
1827  save()=piecejj;
1828  save.write("PPPpiecejjFFT.xmp");
1829  std::cout << "FFT shift=" << fftShift.transpose() << std::endl;
1830  std::cout << "Corr FFT=" << corrFFT << std::endl;
1831  if (corrFFT>corrOriginal)
1832  {
1833  if (XX(rjj)-halfSize>=STARTINGX(Ijj) &&
1834  XX(rjj)+halfSize<=FINISHINGX(Ijj) &&
1835  YY(rjj)-halfSize>=STARTINGY(Ijj) &&
1836  YY(rjj)+halfSize<=FINISHINGY(Ijj))
1837  {
1838  FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
1839  A2D_ELEM(piecejj,i,j)=A2D_ELEM(Ijj,
1840  (int)(YY(rjj)+i),(int)(XX(rjj)+j));
1841  save()=piecejj;
1842  save.write("PPPpiecejjNew.xmp");
1843  }
1844  }
1845  }
1846  }
1847  }
1848 
1849  bool retval=refineLandmark(pieceii,jj,rjj,actualCorrThreshold,
1850  reversed,maxCorr);
1851  return retval;
1852 }
1853 
1854 bool ProgTomographAlignment::refineLandmark(const MultidimArray<double> &pieceii,
1855  int jj, Matrix1D<double> &rjj, double actualCorrThreshold,
1856  bool reversed, double &maxCorr) const
1857 {
1858  int halfSize=XSIZE(pieceii)/2;
1859 
1860  double mean_ii=0, stddev_ii=0;
1861  FOR_ALL_ELEMENTS_IN_ARRAY2D(pieceii)
1862  {
1863  mean_ii+=A2D_ELEM(pieceii,i,j);
1864  stddev_ii+=A2D_ELEM(pieceii,i,j)*A2D_ELEM(pieceii,i,j);
1865  }
1866  mean_ii/=MULTIDIM_SIZE(pieceii);
1867  stddev_ii = stddev_ii / MULTIDIM_SIZE(pieceii) - mean_ii * mean_ii;
1868  stddev_ii *= MULTIDIM_SIZE(pieceii) / (MULTIDIM_SIZE(pieceii) - 1);
1869  stddev_ii = sqrt(static_cast<double>((ABS(stddev_ii))));
1870  if (stddev_ii>XMIPP_EQUAL_ACCURACY)
1871  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(pieceii)
1872  DIRECT_MULTIDIM_ELEM(pieceii,n)=
1873  (DIRECT_MULTIDIM_ELEM(pieceii,n)-mean_ii)/stddev_ii;
1874 
1875  // Try all possible shifts
1876  MultidimArray<double> corr((int)(1.5*(2*halfSize+1)),(int)(1.5*(2*halfSize+1)));
1877  corr.setXmippOrigin();
1878  corr.initConstant(-1.1);
1879  bool accept=false;
1880  double maxval=-1;
1881  MultidimArray<double> piecejj(2*halfSize+1,2*halfSize+1);
1882  piecejj.setXmippOrigin();
1883  if (stddev_ii>XMIPP_EQUAL_ACCURACY)
1884  {
1885  int imax=0, jmax=0;
1886  std::queue< Matrix1D<double> > Q;
1887  Q.push(vectorR2(0,0));
1888  while (!Q.empty())
1889  {
1890  // Get the first position to evaluate
1891  auto shifty=(int)YY(Q.front());
1892  auto shiftx=(int)XX(Q.front());
1893  Q.pop();
1894  if (!corr(shifty,shiftx)<-1)
1895  continue;
1896 
1897  // Select piece in image jj, compute its statistics and normalize
1898  double mean_jj=0, stddev_jj=0;
1899  const MultidimArray<unsigned char> &Ijj=(*img[jj]);
1900  if (XX(rjj)+shiftx+STARTINGX(piecejj)<STARTINGX(Ijj) ||
1901  YY(rjj)+shifty+STARTINGY(piecejj)<STARTINGY(Ijj) ||
1902  XX(rjj)+shiftx+FINISHINGX(piecejj)>FINISHINGX(Ijj) ||
1903  YY(rjj)+shifty+FINISHINGY(piecejj)>FINISHINGY(Ijj))
1904  continue;
1905  FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
1906  {
1907  double pixval=A2D_ELEM(Ijj,
1908  (int)(YY(rjj)+shifty+i),
1909  (int)(XX(rjj)+shiftx+j));
1910  A2D_ELEM(piecejj,i,j)=pixval;
1911  mean_jj+=pixval;
1912  stddev_jj+=pixval*pixval;
1913  }
1914  mean_jj/=MULTIDIM_SIZE(piecejj);
1915  stddev_jj = stddev_jj / MULTIDIM_SIZE(piecejj) - mean_jj * mean_jj;
1916  stddev_jj *= MULTIDIM_SIZE(piecejj) / (MULTIDIM_SIZE(piecejj) - 1);
1917  stddev_jj = sqrt(static_cast<double>((ABS(stddev_jj))));
1918  if (reversed)
1919  piecejj.selfReverseY();
1920 
1921  // Compute the correlation
1922  corr(shifty,shiftx)=0;
1923  double &corrRef=corr(shifty,shiftx);
1924  if (stddev_jj>XMIPP_EQUAL_ACCURACY)
1925  {
1926  double istddev_jj=1.0/stddev_jj;
1927  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(piecejj)
1928  corrRef+=
1929  DIRECT_MULTIDIM_ELEM(pieceii,n)*
1930  (DIRECT_MULTIDIM_ELEM(piecejj,n)-mean_jj)*istddev_jj;
1931  corrRef/=MULTIDIM_SIZE(piecejj);
1932  }
1933 
1934  if (corrRef>maxval)
1935  {
1936  maxval=corrRef;
1937  imax=shifty;
1938  jmax=shiftx;
1939  for (int step=1; step<=5; step+=2)
1940  for (int stepy=-1; stepy<=1; stepy++)
1941  for (int stepx=-1; stepx<=1; stepx++)
1942  {
1943  int newshifty=shifty+stepy*step;
1944  int newshiftx=shiftx+stepx*step;
1945  if (newshifty>=STARTINGY(corr) &&
1946  newshifty<=FINISHINGY(corr) &&
1947  newshiftx>=STARTINGX(corr) &&
1948  newshiftx<=FINISHINGX(corr) &&
1949  (XX(rjj)+newshiftx-halfSize)>=STARTINGX(Ijj) &&
1950  (XX(rjj)+newshiftx+halfSize)<=FINISHINGX(Ijj) &&
1951  (YY(rjj)+newshifty-halfSize)>=STARTINGY(Ijj) &&
1952  (YY(rjj)+newshifty+halfSize)<=FINISHINGY(Ijj))
1953  if (corr(newshifty,newshiftx)<-1)
1954  Q.push(vectorR2(newshiftx,newshifty));
1955  }
1956  }
1957  }
1958 
1959  if (maxval>actualCorrThreshold)
1960  {
1961  XX(rjj)+=jmax;
1962  if (reversed)
1963  YY(rjj)-=imax;
1964  else
1965  YY(rjj)+=imax;
1966  accept=true;
1967  }
1968  if (showRefinement)
1969  {
1970  Image<double> save;
1971  FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
1972  piecejj(i,j)=(*img[jj])((int)(YY(rjj)+i),(int)(XX(rjj)+j));
1973  if (reversed)
1974  piecejj.selfReverseY();
1975  save()=piecejj;
1976  save.write("PPPpiecejj.xmp");
1977  save()=corr;
1978  save.write("PPPcorr.xmp");
1979  std::cout << "jj=" << jj << " rjj=" << rjj.transpose() << std::endl;
1980  std::cout << "imax=" << imax << " jmax=" << jmax << std::endl;
1981  std::cout << "maxval=" << maxval << std::endl;
1982  }
1983  }
1984  if (showRefinement)
1985  {
1986  std::cout << "Accepted=" << accept << std::endl;
1987  std::cout << "Press any key\n";
1988  char c;
1989  std::cin >> c;
1990  }
1991  maxCorr=maxval;
1992  return accept;
1993 }
1994 
1995 /* Refine chain ------------------------------------------------------------ */
1996 //#define DEBUG
1997 bool ProgTomographAlignment::refineChain(LandmarkChain &chain,
1998  double &corrChain)
1999 {
2000 #ifdef DEBUG
2001  std::cout << "Chain for refinement: ";
2002  for (int i=0; i<chain.size(); i++)
2003  std::cout << chain[i].imgIdx << " ";
2004  std::cout << std::endl;
2005 #endif
2006 
2007  for (int K=0; K<2 && (chain.size()>seqLength || useCriticalPoints); K++)
2008  {
2009  sort(chain.begin(), chain.end());
2010  Matrix1D<double> rii(2), rjj(2), newrjj(2);
2011 
2012  if (useCriticalPoints)
2013  {
2014  // compute average piece
2015  int halfSize=XMIPP_MAX(ROUND(localSize*XSIZE(*img[0]))/2,5);
2016  int chainLength=chain.size();
2017  MultidimArray<double> avgPiece(2*halfSize+1,2*halfSize+1), pieceAux;
2018  avgPiece.setXmippOrigin();
2019  pieceAux=avgPiece;
2020  for (int n=0; n<chainLength; n++)
2021  {
2022  int ii=chain[n].imgIdx;
2023  VECTOR_R2(rii,chain[n].x,chain[n].y);
2024  if (XX(rii)-halfSize<STARTINGX(*img[ii]) ||
2025  XX(rii)+halfSize>FINISHINGX(*img[ii]) ||
2026  YY(rii)-halfSize<STARTINGY(*img[ii]) ||
2027  YY(rii)+halfSize>FINISHINGY(*img[ii]))
2028  {
2029  corrChain=0;
2030  return false;
2031  }
2032  FOR_ALL_ELEMENTS_IN_ARRAY2D(avgPiece)
2033  pieceAux(i,j)=(*img[ii])((int)(YY(rii)+i),(int)(XX(rii)+j));
2034  if (isCapillar && ABS(chain[n].imgIdx-chain[0].imgIdx)>Nimg/2)
2035  pieceAux.selfReverseY();
2036  avgPiece+=pieceAux;
2037  }
2038  avgPiece/=chainLength;
2039 #ifdef DEBUG
2040 
2041  Image<double> save;
2042  save()=avgPiece;
2043  save.write("PPPavg.xmp");
2044  std::cout << "Average " << K << " ------------------- \n";
2045  showRefinement=true;
2046 #endif
2047 
2048  // Align all images with respect to this average
2049  corrChain=2;
2050  for (int j=0; j<chainLength; j++)
2051  {
2052  int jj=chain[j].imgIdx;
2053  VECTOR_R2(rjj,chain[j].x,chain[j].y);
2054  double corr;
2055  bool accepted=refineLandmark(avgPiece,jj,rjj,0,false,corr);
2056  if (accepted)
2057  {
2058  chain[j].x=XX(rjj);
2059  chain[j].y=YY(rjj);
2060  corrChain=XMIPP_MIN(corrChain,corr);
2061  }
2062  }
2063 
2064 #ifdef DEBUG
2065  showRefinement=false;
2066 #endif
2067 
2068  }
2069  else
2070  {
2071  // Refine every step
2072  for (int step=2; step<=maxStep; step++)
2073  {
2074  int chainLength=chain.size();
2075 
2076  // Refine forwards every step
2077  int ileft=-1;
2078  for (int i=0; i<chainLength-step; i++)
2079  {
2080  int ii=chain[i].imgIdx;
2081  VECTOR_R2(rii,chain[i].x,chain[i].y);
2082  int jj=chain[i+step].imgIdx;
2083  VECTOR_R2(rjj,chain[i+step].x,chain[i+step].y);
2084  newrjj=rjj;
2085  double corr;
2086  bool accepted=refineLandmark(ii,jj,rii,newrjj,corr,false);
2087  if (((newrjj-rjj).module()<4 && accepted) || useCriticalPoints)
2088  {
2089  chain[i+step].x=XX(newrjj);
2090  chain[i+step].y=YY(newrjj);
2091  }
2092  else
2093  ileft=i;
2094  }
2095 
2096 #ifdef DEBUG
2097  // COSS showRefinement=(step==maxStep);
2098 #endif
2099  // Refine backwards all images
2100  int iright=chainLength;
2101  corrChain=2;
2102  for (int i=chainLength-1; i>=1; i--)
2103  {
2104  int ii=chain[i].imgIdx;
2105  VECTOR_R2(rii,chain[i].x,chain[i].y);
2106  int jj=chain[i-1].imgIdx;
2107  VECTOR_R2(rjj,chain[i-1].x,chain[i-1].y);
2108  newrjj=rjj;
2109  double corr;
2110  bool accepted=refineLandmark(ii,jj,rii,newrjj,corr,false);
2111  corrChain=XMIPP_MIN(corrChain,corr);
2112  if (((newrjj-rjj).module()<4 && accepted) || useCriticalPoints)
2113  {
2114  chain[i-1].x=XX(newrjj);
2115  chain[i-1].y=YY(newrjj);
2116  }
2117  else
2118  iright=i;
2119  }
2120 #ifdef DEBUG
2121  // COSS: showRefinement=false;
2122 #endif
2123 
2124  LandmarkChain refinedChain;
2125  for (int ll=ileft+1; ll<=iright-1; ll++)
2126  refinedChain.push_back(chain[ll]);
2127  chain=refinedChain;
2128  double tilt0=tiltList[chain[0].imgIdx];
2129  double tiltF=tiltList[chain[chain.size()-1].imgIdx];
2130  double lengthThreshold=XMIPP_MAX(3.,FLOOR(seqLength*cos(DEG2RAD(0.5*(tilt0+tiltF)))));
2131  if (chain.size()<lengthThreshold && !useCriticalPoints)
2132  return false;
2133  }
2134  }
2135  }
2136 
2137  // Check that the chain is of the desired length
2138  double tilt0=tiltList[chain[0].imgIdx];
2139  double tiltF=tiltList[chain[chain.size()-1].imgIdx];
2140  double lengthThreshold=XMIPP_MAX(3,FLOOR(seqLength*cos(DEG2RAD(0.5*(tilt0+tiltF)))));
2141  return chain.size()>lengthThreshold;
2142 }
2143 #undef DEBUG
2144 
2145 /* Read/Write landmark set ------------------------------------------------- */
2146 void ProgTomographAlignment::writeLandmarkSet(const FileName &fnLandmark) const
2147 {
2148  std::ofstream fhOut;
2149  fhOut.open(fnLandmark.c_str());
2150  if (!fhOut)
2151  REPORT_ERROR(ERR_IO_NOWRITE,(std::string)"Cannot open "+fnLandmark+" for output");
2152  fhOut << "Point x y slice color "
2153  << MAT_YSIZE(allLandmarksX) << " " << MAT_XSIZE(allLandmarksX) << std::endl;
2154  for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
2155  {
2156  int counter=0;
2157  for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
2158  if (allLandmarksX(j,i)!=XSIZE(*img[0]))
2159  {
2160  fhOut << counter << " \t"
2161  << ROUND(allLandmarksX(j,i)-STARTINGX(*img[0])) << " \t"
2162  << ROUND(allLandmarksY(j,i)-STARTINGY(*img[0])) << " \t"
2163  << i+1 << " \t" << j << std::endl;
2164  counter++;
2165  }
2166  }
2167  fhOut.close();
2168 }
2169 
2170 void ProgTomographAlignment::readLandmarkSet(const FileName &fnLandmark)
2171 {
2172  std::ifstream fhIn;
2173  fhIn.open(fnLandmark.c_str());
2174  if (!fhIn)
2175  REPORT_ERROR(ERR_IO_NOTEXIST,(std::string)"Cannot open "+fnLandmark+" for input");
2176  std::string dummyStr;
2177  int Nlandmark;
2178  fhIn >> dummyStr >> dummyStr >> dummyStr >> dummyStr >> dummyStr
2179  >> Nlandmark >> Nimg;
2180  if (Nlandmark<=0)
2181  REPORT_ERROR(ERR_VALUE_INCORRECT,(std::string)"No landmarks are found in "+fnLandmark);
2182  allLandmarksX.resize(Nlandmark,Nimg);
2183  allLandmarksY.resize(Nlandmark,Nimg);
2184  allLandmarksX.initConstant(XSIZE(*img[0]));
2185  allLandmarksY.initConstant(XSIZE(*img[0]));
2186  fhIn.exceptions ( std::ifstream::eofbit |
2187  std::ifstream::failbit |
2188  std::ifstream::badbit );
2189  while (!fhIn.eof())
2190  {
2191  try
2192  {
2193  int dummyInt, x, y, i, j;
2194  fhIn >> dummyInt >> x >> y >> i >> j;
2195  i=i-1;
2196  allLandmarksX(j,i)=x+STARTINGX(*img[0]);
2197  allLandmarksY(j,i)=y+STARTINGY(*img[0]);
2198  }
2199  catch (std::ifstream::failure e)
2200  {
2201  // Do nothing with this line
2202  }
2203  }
2204  fhIn.close();
2205  std::cout << "The file " << fnLandmark << " has been read for the landmarks\n"
2206  << Nlandmark << " landmarks are read\n";
2207 }
2208 
2209 /* Write affine transformations -------------------------------------------- */
2210 void ProgTomographAlignment::writeTransformations(
2211  const FileName &fnTransformations) const
2212 {
2213  std::ofstream fhOut;
2214  fhOut.open(fnTransformations.c_str());
2215  if (!fhOut)
2216  REPORT_ERROR(ERR_IO_NOWRITE,(std::string)"Cannot open "+fnTransformations+" for output");
2217  int imax=affineTransformations.size();
2218  int counter=0;
2219  for (int i=0; i<imax; i++)
2220  {
2221  int jmax=affineTransformations[i].size();
2222  for (int j=i+1; j<jmax; j++)
2223  if (MAT_XSIZE(affineTransformations[i][j])!=0)
2224  {
2225  fhOut << tiltList[i] << "\t0.0\t"
2226  << affineTransformations[i][j](0,0) << "\t"
2227  << affineTransformations[i][j](1,0) << "\t"
2228  << affineTransformations[i][j](0,1) << "\t"
2229  << affineTransformations[i][j](1,1) << "\t"
2230  << affineTransformations[i][j](0,2) << "\t"
2231  << affineTransformations[i][j](1,2) << std::endl;
2232  counter++;
2233  }
2234  }
2235  if (counter==imax-1)
2236  {
2237  fhOut << tiltList[tiltList.size()-1] << "\t0.0\t1.0\t0.0\t0.0\t1.0\t0.0\t0.0\n";
2238  fhOut << "1 0 0 1 0 0\n";
2239  }
2240  fhOut.close();
2241 }
2242 
2243 /* Align images ------------------------------------------------------------ */
2244 //#define DEBUG
2245 void ProgTomographAlignment::alignImages(const Alignment &alignment)
2246 {
2247  // Correct all landmarks
2248  Matrix1D<double> r(2);
2249  for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
2250  {
2251  Matrix2D<double> R;
2252  rotation2DMatrix(90-alignment.rot+alignment.psi(i),R,false);
2253  for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
2254  if (allLandmarksX(j,i)!=XSIZE(*img[0]))
2255  {
2256  VECTOR_R2(r,allLandmarksX(j,i),allLandmarksY(j,i));
2257  r=R*(r-alignment.di[i]+alignment.diaxis[i]);
2258  allLandmarksX(j,i)=XX(r);
2259  allLandmarksY(j,i)=YY(r);
2260  }
2261  }
2262 
2263  // Compute the average height of the 3D landmarks seen at 0 degrees
2264  Matrix1D<double> axis;
2265  Euler_direction(alignment.rot, alignment.tilt, 0, axis);
2266  double z0=0, z0N=0;
2267  for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
2268  if (allLandmarksX(j,iMinTilt)!=XSIZE(*img[0]))
2269  {
2270  Matrix2D<double> Raxismin, Rmin, RtiltYmin;
2271  rotation3DMatrix(tiltList[iMinTilt],axis,Raxismin,false);
2272  rotation2DMatrix(90-alignment.rot+alignment.psi(iMinTilt),Rmin);
2273  rotation3DMatrix(-tiltList[iMinTilt],'Y',RtiltYmin,false);
2274  Matrix1D<double> rjp=RtiltYmin*Rmin*Raxismin*alignment.rj[j];
2275  z0+=ZZ(rjp);
2276  z0N++;
2277  }
2278  if (z0N==0)
2279  REPORT_ERROR(ERR_VALUE_INCORRECT,"There is no landmark at 0 degrees");
2280  z0/=z0N;
2281  std::cout << "Average height of the landmarks at 0 degrees=" << z0 << std::endl;
2282  MetaDataVec DF;
2283 
2284  std::unique_ptr<MetaDataVec::id_iterator> idIter;
2285 
2286  if (!fnSelOrig.empty())
2287  {
2288  idIter = memoryUtils::make_unique<MetaDataVec::id_iterator>(SForig.ids().begin());
2289  }
2290  DF.setComment("First shift by -(shiftX,shiftY), then rotate by psi");
2291 
2292  MultidimArray<double> mask;
2293  MultidimArray<int> iMask;
2294  Image<double> I;
2295  Matrix2D<double> M, M1, M2, M3;
2296  FileName fn_corrected;
2297 
2298  for (int n=0;n<Nimg; n++)
2299  {
2300  // Align the normal image
2301  I.read(name_list[n]);
2302  mask.initZeros(I());
2303  FOR_ALL_ELEMENTS_IN_ARRAY2D(I())
2304  if (I(i,j)!=0)
2305  mask(i,j)=1;
2306  translation2DMatrix(vectorR2(-z0*sin(DEG2RAD(tiltList[n])),0),M1);
2307  rotation2DMatrix(90-alignment.rot+alignment.psi(n),M2);
2308  translation2DMatrix(-(alignment.di[n]+alignment.diaxis[n]),M3);
2309  M=M1*M2*M3;
2310  selfApplyGeometry(xmipp_transformation::BSPLINE3,I(),M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
2311  selfApplyGeometry(xmipp_transformation::LINEAR,mask,M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
2312  mask.binarize(0.5);
2313  typeCast(mask,iMask);
2314  double minval, maxval, avg, stddev;
2315  computeStats_within_binary_mask(iMask,I(),minval, maxval, avg, stddev);
2316  FOR_ALL_ELEMENTS_IN_ARRAY2D(iMask)
2317  if (iMask(i,j)==0)
2318  I(i,j)=0;
2319  else if (!dontNormalize)
2320  I(i,j)-=avg;
2321  double rot=0;
2322  double tilt=tiltList[n];
2323  double psi=0;
2324  I.setEulerAngles(rot, tilt, psi);
2325  fn_corrected.compose(n+1, fnRoot+"_corrected_", "stk");
2326  I.write(fn_corrected);
2327 
2328  // Align the original image
2329  if (fnSelOrig!="")
2330  {
2331  FileName auxFn;
2332  SForig.getValue( MDL_IMAGE, auxFn, **idIter);
2333  Image<double> Iorig;
2334  Iorig.read( auxFn );
2335  //SForig.nextObject();
2336  ++(*idIter);
2337  mask.initZeros(Iorig());
2338  FOR_ALL_ELEMENTS_IN_ARRAY2D(Iorig())
2339  if (Iorig(i,j)!=0)
2340  mask(i,j)=1;
2341  translation2DMatrix(vectorR2(-z0*sin(DEG2RAD(tiltList[n]))*
2342  (((double)XSIZE(Iorig()))/XSIZE(I())),0),M1);
2343  rotation2DMatrix(90-alignment.rot+alignment.psi(n),M2);
2344  translation2DMatrix(-(alignment.di[n]+alignment.diaxis[n])*
2345  (((double)XSIZE(Iorig()))/XSIZE(I())),M3);
2346  M=M1*M2*M3;
2347  selfApplyGeometry(xmipp_transformation::BSPLINE3,Iorig(),M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
2348  selfApplyGeometry(xmipp_transformation::LINEAR,mask,M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
2349  mask.binarize(0.5);
2350  typeCast(mask,iMask);
2351  computeStats_within_binary_mask(iMask,Iorig(),minval, maxval,
2352  avg, stddev);
2353  FOR_ALL_ELEMENTS_IN_ARRAY2D(iMask)
2354  if (iMask(i,j)==0)
2355  Iorig(i,j)=0;
2356  else if (!dontNormalize)
2357  Iorig(i,j)-=avg;
2358  Iorig.setEulerAngles(rot, tilt, psi);
2359  fn_corrected.compose(n+1, fnRoot+"_corrected_originalsize_", "stk");
2360  Iorig.write(fn_corrected);
2361  }
2362 
2363  // Prepare data for the docfile
2364  size_t id = DF.addObject();
2365  DF.setValue(MDL_IMAGE, fn_corrected, id);
2366  DF.setValue(MDL_ANGLE_PSI, 90.-alignment.rot+alignment.psi(n), id);
2367  DF.setValue(MDL_SHIFT_X, XX(alignment.di[n]+alignment.diaxis[n]), id);
2368  DF.setValue(MDL_SHIFT_Y, YY(alignment.di[n]+alignment.diaxis[n]), id);
2369  }
2370  DF.write(fnRoot+"_correction_parameters.txt");
2371 #ifdef DEBUG
2372 
2373  Image<double> save;
2374  save().initZeros(*img[0]);
2375  save().setXmippOrigin();
2376  for (int j=0; j<YSIZE(allLandmarksX); j++)
2377  if (allLandmarksX(j,iMinTilt)!=XSIZE(*img[0]))
2378  {
2379  Matrix2D<double> Raxismin=rotation3DMatrix(tiltList[iMinTilt],axis);
2380  Raxismin.resize(3,3);
2381  Matrix2D<double> Rmin=rotation2DMatrix(90-alignment.rot+alignment.psi(iMinTilt));
2382  Matrix2D<double> RtiltYmin=rotation3DMatrix(-tiltList[iMinTilt],'Y');
2383  RtiltYmin.resize(3,3);
2384  Matrix1D<double> rjp=RtiltYmin*Rmin*Raxismin*alignment.rj[j];
2385  std::cout << rjp.transpose() << std::endl;
2386  for (int i=0; i<XSIZE(allLandmarksX); i++)
2387  if (allLandmarksX(j,i)!=XSIZE(*img[0]))
2388  {
2389  save(allLandmarksY(j,i),allLandmarksX(j,i))=ABS(i-iMinTilt);
2390 
2391  /*
2392  Matrix2D<double> Raxis=rotation3DMatrix(tiltList[i],axis);
2393  Raxis.resize(3,3);
2394  Matrix2D<double> R=rotation2DMatrix(90-alignment.rot+alignment.psi(i));
2395  Matrix1D<double> p=R*Raxis*alignment.rj[j];
2396  if (YY(p)>=STARTINGY(save()) && YY(p)<=FINISHINGY(save()) &&
2397  XX(p)>=STARTINGX(save()) && XX(p)<=FINISHINGX(save()))
2398  save(YY(p),XX(p))=-ABS(i-iMinTilt);
2399  */
2400  Matrix2D<double> RtiltY=rotation3DMatrix(-tiltList[i],'Y');
2401  RtiltY.resize(3,3);
2402  Matrix1D<double> p=RtiltY*rjp;
2403  if (YY(p)>=STARTINGY(save()) && YY(p)<=FINISHINGY(save()) &&
2404  XX(p)>=STARTINGX(save()) && XX(p)<=FINISHINGX(save()))
2405  save(YY(p),XX(p))=-ABS(i-iMinTilt);
2406  }
2407  }
2408  save.write("PPPmovementsLandmarks0.xmp");
2409 #endif
2410 }
2411 #undef DEBUG
2412 
2413 /* Produce information from landmarks -------------------------------------- */
2414 void ProgTomographAlignment::produceInformationFromLandmarks()
2415 {
2416  // Produce V sets
2417  std::vector<int> emptyVector;
2418 
2419  Vseti.clear();
2420  for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
2421  Vseti.push_back(emptyVector);
2422 
2423  Vsetj.clear();
2424  for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
2425  Vsetj.push_back(emptyVector);
2426 
2427  for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
2428  for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
2429  if (allLandmarksX(j,i)!=XSIZE(*img[0]))
2430  {
2431  Vseti[i].push_back(j);
2432  Vsetj[j].push_back(i);
2433  }
2434 
2435  // Count the number of landmarks per image and average projection
2436  // of all landmarks in a given image
2437  ni.initZeros(Nimg);
2438  barpi.clear();
2439  for (int i=0; i<Nimg; i++)
2440  {
2441  ni(i)=static_cast<int>(Vseti[i].size());
2442  Matrix1D<double> pi(2);
2443  for (int jj=0; jj<ni(i); jj++)
2444  {
2445  int j=Vseti[i][jj];
2446  XX(pi)+=allLandmarksX(j,i);
2447  YY(pi)+=allLandmarksY(j,i);
2448  }
2449  pi/=ni(i);
2450  barpi.push_back(pi);
2451  }
2452 }
2453 
2454 /* Remove outliers ---------------------------------------------------------*/
2455 void ProgTomographAlignment::removeOutlierLandmarks(
2456  const Alignment &alignment)
2457 {
2458  std::cout << "Removing outliers ...\n";
2459 
2460  // Compute threshold for outliers
2461  Histogram1D hist;
2462  compute_hist(alignment.errorLandmark, hist, 100);
2463  // double threshold0=hist.percentil(10);
2464  double thresholdF=hist.percentil(90);
2465 
2466  // Identify outliers
2467  std::vector<bool> validLandmark;
2468  int invalidCounter=0;
2469  int Nlandmark=MAT_YSIZE(allLandmarksX);
2470  for (int j=0; j<Nlandmark; j++)
2471  if (//COSS alignment.errorLandmark(j)<threshold0 ||
2472  alignment.errorLandmark(j)>thresholdF)
2473  {
2474  validLandmark.push_back(false);
2475  invalidCounter++;
2476  }
2477  else
2478  validLandmark.push_back(true);
2479 
2480  // Remove outliers
2481  Matrix2D<double> newAllLandmarksX(Nlandmark-invalidCounter,Nimg);
2482  Matrix2D<double> newAllLandmarksY(Nlandmark-invalidCounter,Nimg);
2483  int jj=0;
2484  for (int j=0; j<Nlandmark; j++)
2485  {
2486  if (!validLandmark[j])
2487  continue;
2488  for (int i=0; i<Nimg; i++)
2489  {
2490  newAllLandmarksX(jj,i)=allLandmarksX(j,i);
2491  newAllLandmarksY(jj,i)=allLandmarksY(j,i);
2492  }
2493  jj++;
2494  }
2495  allLandmarksX=newAllLandmarksX;
2496  allLandmarksY=newAllLandmarksY;
2497 
2498  std::cout << invalidCounter << " out of " << Nlandmark
2499  << " landmarks have been removed\n";
2500 }
2501 
2502 /* Run --------------------------------------------------------------------- */
2503 namespace TomographAlignment
2504 {
2505 const ProgTomographAlignment* global_prm;
2506 }
2507 
2508 double wrapperError(double *p, void *prm)
2509 {
2510  Alignment alignment(TomographAlignment::global_prm);
2511  alignment=*(TomographAlignment::global_prm->bestPreviousAlignment);
2512  alignment.rot=p[1];
2513  alignment.tilt=p[2];
2514  return alignment.optimizeGivenAxisDirection();
2515 }
2516 
2517 #define DEBUG
2518 void ProgTomographAlignment::run()
2519 {
2520  produceSideInfo();
2521  std::cerr << "generateLandmarkSet"<< std::endl;
2522  generateLandmarkSet();
2523  std::cerr << "produceInformationFromLandmarks" << std::endl;
2524  produceInformationFromLandmarks();
2525  std::cerr << "alignment" << std::endl;
2526  auto *alignment=new Alignment(this);
2527 
2528  // Exhaustive search for rot
2529  double bestError=0, bestRot=-1;
2530  for (double rot=0; rot<=180-deltaRot; rot+=deltaRot)
2531  {
2532  alignment->clear();
2533  alignment->rot=rot;
2534  double error=alignment->optimizeGivenAxisDirection();
2535 #ifdef DEBUG
2536 
2537  std::cout << "rot= " << rot
2538  << " error= " << error << std::endl;
2539 #endif
2540 
2541  if (bestRot<0 || bestError>error)
2542  {
2543  bestRot=rot;
2544  bestError=error;
2545  *bestPreviousAlignment=*alignment;
2546  }
2547  }
2548  delete alignment;
2549  std::cout << "Best rot=" << bestRot
2550  << " Best error=" << bestError << std::endl;
2551 
2552  // Continuous optimization for the axis direction
2553  Matrix1D<double> axisAngles(2), steps(2);
2554  axisAngles(0)=bestRot;
2555  axisAngles(1)=90;
2556  steps.initConstant(1);
2557  if (!optimizeTiltAngle)
2558  steps(1)=0;
2559  double fitness;
2560  int iter;
2561  TomographAlignment::global_prm=this;
2562  powellOptimizer(axisAngles,1,2,&wrapperError, nullptr,
2563  0.01,fitness,iter,steps,true);
2564 
2565  // Outlier removal
2566  for (int i=0; i<3; i++)
2567  {
2568  // Compute the best alignment
2569  bestPreviousAlignment->rot=axisAngles(0);
2570  bestPreviousAlignment->tilt=axisAngles(1);
2571  fitness=bestPreviousAlignment->optimizeGivenAxisDirection();
2572  bestPreviousAlignment->computeErrorForLandmarks();
2573 
2574  // Remove landmarks that are outliers in the current model
2575  removeOutlierLandmarks(*bestPreviousAlignment);
2576  produceInformationFromLandmarks();
2577  delete bestPreviousAlignment;
2578  bestPreviousAlignment=new Alignment(this);
2579  bestPreviousAlignment->rot=axisAngles(0);
2580  bestPreviousAlignment->tilt=axisAngles(1);
2581  fitness=bestPreviousAlignment->optimizeGivenAxisDirection();
2582 
2583  // Optimize again
2584  powellOptimizer(axisAngles,1,2,&wrapperError,nullptr,
2585  0.01,fitness,iter,steps,true);
2586  }
2587  bestPreviousAlignment->rot=axisAngles(0);
2588  bestPreviousAlignment->tilt=axisAngles(1);
2589  fitness=bestPreviousAlignment->optimizeGivenAxisDirection();
2590 
2591  // Save the alignment
2592  writeLandmarkSet(fnRoot+"_good_landmarks.txt");
2593  std::ofstream fh_out;
2594  fh_out.open((fnRoot+"_alignment.txt").c_str());
2595  if (!fh_out)
2596  REPORT_ERROR(ERR_IO_NOWRITE,
2597  (std::string)"Cannot open "+fnRoot+"_alignment.txt for output");
2598  fh_out << *bestPreviousAlignment;
2599  fh_out.close();
2600 
2601  // Correct the input images
2602  alignImages(*bestPreviousAlignment);
2603  writeLandmarkSet(fnRoot+"_good_landmarks_corrected.txt");
2604  delete bestPreviousAlignment;
2605 }
2606 #undef DEBUG
2607 
2608 /* Optimize for rot -------------------------------------------------------- */
2609 //#define DEBUG
2610 double Alignment::optimizeGivenAxisDirection()
2611 {
2612  double bestError=1e38;
2613  bool firstIteration=true, finish=false;
2614  int Niterations=0;
2615  computeGeometryDependentOfAxis();
2616  do
2617  {
2618  computeGeometryDependentOfRotation();
2619  double error=computeError();
2620 #ifdef DEBUG
2621 
2622  std::cout << "it=" << Niterations << " Error= " << error
2623  << " raxis=" << raxis.transpose() << std::endl;
2624 #endif
2625 
2626  updateModel();
2627  Niterations++;
2628  if (firstIteration)
2629  {
2630  bestError=error;
2631  firstIteration=false;
2632  }
2633  else
2634  {
2635  finish=((error>bestError) || (Niterations>1000) ||
2636  (ABS(error-bestError)/bestError<0.001)) && Niterations>20;
2637  if (error<bestError)
2638  bestError=error;
2639  }
2640  }
2641  while (!finish);
2642  return bestError;
2643 }
2644 #undef DEBUG
2645 
2646 /* Compute the geometry part corresponding to axis ------------------------- */
2647 void Alignment::computeGeometryDependentOfAxis()
2648 {
2649  Matrix1D<double> axis;
2650  Euler_direction(rot, tilt, 0, axis);
2651 
2652  // Compute Aip, Aipt
2653  Matrix2D<double> Raxis, I, Ip, B;
2654  Binvraxis.initZeros(3,3);
2655  I.initIdentity(3);
2656  Ip=I;
2657  Ip(2,2)=0;
2658  for (int i=0; i<Nimg; i++)
2659  {
2660  rotation3DMatrix(prm->tiltList[i],axis,Raxis,false);
2661  Aipt[i](0,0)=Aip[i](0,0)=Raxis(0,0);
2662  Aipt[i](1,0)=Aip[i](0,1)=Raxis(0,1);
2663  Aipt[i](2,0)=Aip[i](0,2)=Raxis(0,2);
2664  Aipt[i](0,1)=Aip[i](1,0)=Raxis(1,0);
2665  Aipt[i](1,1)=Aip[i](1,1)=Raxis(1,1);
2666  Aipt[i](2,1)=Aip[i](1,2)=Raxis(1,2);
2667 
2668  B=I-Raxis;
2669  B=B.transpose()*Ip*B;
2670  Binvraxis+=((double)prm->ni(i))*(B+B.transpose());
2671 
2672  B1i[i]=I-Raxis.transpose();
2673  B2i[i]=B1i[i]*Ip*Raxis;
2674  /*COSS:
2675  std::cout << "Dependent on axis i=" << i << std::endl
2676  << "tilt[i]=" << prm->tiltList[i] << std::endl
2677  << "Raxis=\n" << Raxis
2678  << "Aip[i]=\n" << Aip[i]
2679  << "B1i[i]=\n" << B1i[i] << std::endl
2680  << "B2i[i]=\n" << B2i[i] << std::endl;
2681  */
2682  }
2683  Binvraxis=Binvraxis.inv();
2684 }
2685 
2686 /* Compute the geometry part corresponding to axis ------------------------- */
2687 void Alignment::computeGeometryDependentOfRotation()
2688 {
2689  // Compute Ai, Ait
2690  Matrix2D<double> Rinplane, I(2,3);
2691  I(0,0)=I(1,1)=1;
2692  for (int i=0; i<Nimg; i++)
2693  {
2694  rotation3DMatrix(-psi(i),'Z',Rinplane,false);
2695  B1i[i]=B1i[i]*Rinplane.transpose();
2696 
2697  Rinplane.resize(2,2);
2698  Ai[i]=Rinplane*Aip[i];
2699  Ait[i]=Ai[i].transpose();
2700  diaxis[i]=Rinplane*(I-Aip[i])*raxis;
2701  /* COSS:
2702  std::cout << "Dependent on rotation i=" << i << std::endl
2703  << "B1i[i]=\n" << B1i[i] << std::endl
2704  << "I-Raxis=\n" << (I-Aip[i])
2705  << "raxis=" << raxis.transpose() << std::endl;
2706  std::cout << "diaxis[" << i << "]=" << diaxis[i].transpose() << std::endl;
2707  */
2708  }
2709 }
2710 
2711 /* Compute error ----------------------------------------------------------- */
2712 double Alignment::computeError() const
2713 {
2714  Matrix1D<double> pijp;
2715  double error=0;
2716  double N=0;
2717  for (int i=0; i<Nimg; i++)
2718  {
2719  int jjmax=prm->Vseti[i].size();
2720  for (int jj=0; jj<jjmax; jj++)
2721  {
2722  int j=prm->Vseti[i][jj];
2723  pijp=Ai[i]*rj[j]+di[i]+diaxis[i];
2724  MAT_ELEM(allLandmarksPredictedX,j,i)=XX(pijp);
2725  MAT_ELEM(allLandmarksPredictedY,j,i)=YY(pijp);
2726  double diffx=MAT_ELEM(prm->allLandmarksX,j,i)-XX(pijp);
2727  double diffy=MAT_ELEM(prm->allLandmarksY,j,i)-YY(pijp);
2728  error+=diffx*diffx+diffy*diffy;
2729  N++;
2730  }
2731  }
2732  return sqrt(error/N);
2733 }
2734 
2735 void Alignment::computeErrorForLandmarks()
2736 {
2737  Nlandmark=MAT_YSIZE(prm->allLandmarksX);
2738  errorLandmark.initZeros(Nlandmark);
2739  for (int j=0; j<Nlandmark; j++)
2740  {
2741  int counterj=0;
2742  for (int i=0; i<Nimg; i++)
2743  {
2744  if (prm->allLandmarksX(j,i)!=XSIZE(*(prm->img[0])))
2745  {
2746  double diffx=MAT_ELEM(prm->allLandmarksX,j,i)-
2747  MAT_ELEM(allLandmarksPredictedX,j,i);
2748  double diffy=MAT_ELEM(prm->allLandmarksY,j,i)-
2749  MAT_ELEM(allLandmarksPredictedY,j,i);
2750  errorLandmark(j)+=sqrt(diffx*diffx+diffy*diffy);
2751  counterj++;
2752  }
2753  }
2754  errorLandmark(j)/=counterj;
2755  }
2756 }
2757 
2758 /* Update model ------------------------------------------------------------ */
2759 //#define DEBUG
2760 //UTILS
2761 void Alignment::updateModel()
2762 {
2763  Matrix1D<double> pij(2), piN(2);
2764  Matrix2D<double> A(3,3), AitAi;
2765  Matrix1D<double> b(3);
2766 
2767 #ifdef DEBUG
2768 
2769  std::cout << "Step 0: error=" << computeError() << std::endl;
2770 #endif
2771 
2772  // Update the 3D positions of the landmarks
2773  for (int j=0; j<Nlandmark; j++)
2774  {
2775  A.initZeros();
2776  b.initZeros();
2777  // Compute the part corresponding to Vj
2778  int iimax=prm->Vsetj[j].size();
2779  for (int ii=0; ii<iimax; ii++)
2780  {
2781  int i=prm->Vsetj[j][ii];
2782  XX(pij)=MAT_ELEM(prm->allLandmarksX,j,i);
2783  YY(pij)=MAT_ELEM(prm->allLandmarksY,j,i);
2784  A+=Ait[i]*Ai[i];
2785  b+=Ait[i]*(pij-(di[i]+diaxis[i]));
2786  }
2787 
2788  // Update rj[j]
2789  rj[j]=A.inv()*b;
2790  }
2791 #ifdef DEBUG
2792  std::cout << "Step 1: error=" << computeError() << std::endl;
2793 #endif
2794 
2795  // Compute the average landmarks seen in each image
2796  for (int i=0; i<Nimg; i++)
2797  {
2798  barri[i].initZeros();
2799  for (int jj=0; jj<prm->ni(i); jj++)
2800  {
2801  int j=prm->Vseti[i][jj];
2802  barri[i]+=rj[j];
2803  }
2804  barri[i]/=prm->ni(i);
2805  }
2806 
2807  // Update shifts
2808  if (prm->isCapillar)
2809  {
2810  // Update the raxis
2811  Matrix1D<double> tmp2(3), Pij(3);
2812  for (int i=0; i<Nimg; i++)
2813  {
2814  for (int jj=0; jj<prm->ni(i); jj++)
2815  {
2816  int j=prm->Vseti[i][jj];
2817  XX(Pij)=MAT_ELEM(prm->allLandmarksX,j,i);
2818  YY(Pij)=MAT_ELEM(prm->allLandmarksY,j,i);
2819  tmp2+=B1i[i]*Pij-B2i[i]*rj[j];
2820  }
2821  }
2822  tmp2*=2.0;
2823  tmp2=Binvraxis*tmp2;
2824  XX(raxis)=XX(tmp2);
2825  YY(raxis)=YY(tmp2);
2826  computeGeometryDependentOfRotation();
2827  }
2828  else
2829  {
2830  // Update the individual di
2831  for (int i=0; i<Nimg; i++)
2832  if (i!=prm->iMinTilt)
2833  {
2834  di[i] = prm->barpi[i]-Ai[i]*barri[i]-diaxis[i];
2835 
2836  if (di[i].isAnyNaN())
2837  di[i].initZeros();
2838  }
2839  }
2840 #ifdef DEBUG
2841  std::cout << "Step 2: error=" << computeError() << std::endl;
2842 #endif
2843 
2844  // Update rotations
2845  if (prm->psiMax>0)
2846  {
2847  Matrix2D<double> Ri(2,2);
2848  Matrix2D<double> Aiprj(2,1), Aiprjt(1,2), dim(2,1), Pij(2,1);
2849  for (int i=0; i<Nimg; i++)
2850  {
2851  Ri.initZeros();
2852  dim.fromVector(di[i]+diaxis[i]);
2853  for (int jj=0; jj<prm->ni(i); jj++)
2854  {
2855  int j=prm->Vseti[i][jj];
2856  Aiprj.fromVector(Aip[i]*rj[j]);
2857  Aiprjt=Aiprj.transpose();
2858 
2859  MAT_ELEM(Pij,0,0)=MAT_ELEM(prm->allLandmarksX,j,i);
2860  MAT_ELEM(Pij,1,0)=MAT_ELEM(prm->allLandmarksY,j,i);
2861  Ri+=(dim-Pij)*Aiprjt;
2862  }
2863  psi(i)=CLIP(RAD2DEG(atan(((Ri(0,1)-Ri(1,0))/(Ri(0,0)+Ri(1,1))))),
2864  -(prm->psiMax),prm->psiMax);
2865  Matrix2D<double> Rinplane;
2866  rotation3DMatrix(-psi(i),'Z',Rinplane,false);
2867  }
2868  }
2869 #ifdef DEBUG
2870  std::cout << "Step 3: error=" << computeError() << std::endl;
2871 #endif
2872 
2873  // Update the rotation dependent part
2874  computeGeometryDependentOfRotation();
2875 }
2876 #undef DEBUG
2877 
2878 /* Print ------------------------------------------------------------------- */
2879 std::ostream& operator << (std::ostream &out, Alignment &alignment)
2880 {
2881  out << "Alignment parameters ===========================================\n"
2882  << "rot=" << alignment.rot << " tilt=" << alignment.tilt << std::endl
2883  << "raxis=" << alignment.raxis.transpose() << std::endl;
2884  out << "Images ---------------------------------------------------------\n";
2885  for (int i=0; i<alignment.Nimg; i++)
2886  out << "Image " << i << " psi= " << alignment.psi(i)
2887  << " di= " << alignment.di[i].transpose()
2888  << " diaxis= " << alignment.diaxis[i].transpose()
2889  << " (" << alignment.prm->ni(i) << ")\n";
2890  out << "Landmarks ------------------------------------------------------\n";
2891  alignment.computeErrorForLandmarks();
2892  for (int j=0; j<alignment.Nlandmark; j++)
2893  out << "Landmark " << j << " rj= " << alignment.rj[j].transpose()
2894  << " " << alignment.errorLandmark(j) << std::endl;
2895  return out;
2896 }
FOR_ALL_ELEMENTS_IN_MATRIX2D
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
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
VECTOR_R2
#define VECTOR_R2(v, x, y)
Definition: matrix1d.h:112
ProgTomographAlignment::globalAffine
bool globalAffine
Look for local affine transformation.
Definition: tomo_align_tilt_series.h:92
MAT_YSIZE
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
init_progress_bar
void init_progress_bar(long total)
Definition: xmipp_funcs.cpp:782
MultidimArray::selfWindow
void selfWindow(int n0, int z0, int y0, int x0, int nF, int zF, int yF, int xF, T init_value=0)
Definition: multidim_array.h:538
numerical_tools.h
ProgTomographAlignment::refineLandmark
bool refineLandmark(int ii, int jj, const Matrix1D< double > &rii, Matrix1D< double > &rjj, double &maxCorr, bool tryFourier) const
Definition: tomo_align_tilt_series.cpp:1735
Alignment::di
std::vector< Matrix1D< double > > di
Definition: tomo_align_tilt_series.h:299
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
threadgenerateLandmarkSetGrid
void * threadgenerateLandmarkSetGrid(void *args)
Definition: tomo_align_tilt_series.cpp:1047
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
metadata_sql.h
ProgTomographAlignment::maxShiftPercentage
double maxShiftPercentage
Maxshift percentage.
Definition: tomo_align_tilt_series.h:111
Alignment::prm
const ProgTomographAlignment * prm
Definition: tomo_align_tilt_series.h:298
ThreadGenerateLandmarkSetParams::chainList
std::vector< LandmarkChain > * chainList
Definition: tomo_align_tilt_series.cpp:1044
translation2DMatrix
template void translation2DMatrix(const Matrix1D< float > &, Matrix2D< float > &, bool inverse)
Alignment::rj
std::vector< Matrix1D< double > > rj
Definition: tomo_align_tilt_series.h:300
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
iteration
n The following was calculated during iteration
Definition: numerical_recipes.cpp:8201
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
alignImages
double alignImages(const MultidimArray< double > &Iref, const AlignmentTransforms &IrefTransforms, MultidimArray< double > &I, Matrix2D< double > &M, bool wrap, AlignmentAux &aux, CorrelationAux &aux2, RotationalCorrelationAux &aux3)
Definition: filters.cpp:2047
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
ProgTomographAlignment::seqLength
size_t seqLength
Sequence Length.
Definition: tomo_align_tilt_series.h:101
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
c
doublereal * c
Definition: numerical_recipes.cpp:2230
MULTIDIM_SIZE
#define MULTIDIM_SIZE(v)
Definition: multidim_array_base.h:115
Euler_direction
void Euler_direction(double alpha, double beta, double gamma, Matrix1D< double > &v)
Definition: geometry.cpp:721
DEG2RAD
#define DEG2RAD(d)
Definition: xmipp_macros.h:312
generateMask
void generateMask(const MultidimArray< double > &I, MultidimArray< unsigned char > &mask, int patchSize)
Definition: tomo_align_tilt_series.cpp:40
Alignment::errorLandmark
MultidimArray< double > errorLandmark
Definition: tomo_align_tilt_series.h:451
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ProgTomographAlignment::maskImg
std::vector< MultidimArray< unsigned char > * > maskImg
Definition: tomo_align_tilt_series.h:240
ProgTomographAlignment::run
void run()
Run: do the real work.
Definition: tomo_align_tilt_series.cpp:2518
BinaryCircularMask
void BinaryCircularMask(MultidimArray< int > &mask, double radius, int mode, double x0, double y0, double z0)
Definition: mask.cpp:193
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
y
static double * y
Definition: numerical_recipes.cpp:8487
ERR_IO_NOWRITE
Couldn&#39;t write to file.
Definition: xmipp_error.h:140
BANDPASS
#define BANDPASS
Definition: fourier_filter.h:79
AffineFitness::Mask2
std::vector< MultidimArray< double > * > Mask2
Definition: tomo_align_tilt_series.cpp:78
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
readFloatList
void readFloatList(const char *str, int N, std::vector< T > &v)
Definition: args.h:104
MDL_SHIFT_X
Shift for the image in the X axis (double)
Definition: metadata_label.h:431
applyGeometry
void applyGeometry(int SplineDegree, MultidimArray< std::complex< double > > &V2, const MultidimArray< std::complex< double > > &V1, const Matrix2D< double > &A, bool inv, bool wrap, std::complex< double > outside, MultidimArray< double > *BcoeffsPtr)
Definition: transformations.cpp:668
ImageBase::write
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
Definition: xmipp_image_base.cpp:292
correlationIndex
double correlationIndex(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)
Definition: multidim_array.h:3974
AffineFitness::minAllowed
Matrix1D< double > minAllowed
Definition: tomo_align_tilt_series.cpp:73
FileName::compose
void compose(const String &str, const size_t no, const String &ext="")
Definition: xmipp_filename.cpp:52
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
MDL_ANGLE_PSI
Special label to be used when gathering MDs in MpiMetadataPrograms.
Definition: metadata_label.h:51
ProgTomographAlignment::maxIterDE
int maxIterDE
Maximum number of iterations in Differential Evolution.
Definition: tomo_align_tilt_series.h:114
MultidimArray::initConstant
void initConstant(T val)
Definition: multidim_array.h:2705
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
ImageBase::name
const FileName & name() const
Definition: xmipp_image_base.h:500
Histogram1D::percentil
double percentil(double percent_mass)
Definition: histogram.cpp:160
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
powellOptimizer
void powellOptimizer(Matrix1D< double > &p, int i0, int n, double(*f)(double *x, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show)
Definition: numerical_tools.cpp:44
Matrix2D::transpose
Matrix2D< T > transpose() const
Definition: matrix2d.cpp:1314
z0
#define z0
M3
#define M3
Definition: numerical_recipes.cpp:53
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
q1
#define q1
ProgTomographAlignment::fnRoot
FileName fnRoot
Output root.
Definition: tomo_align_tilt_series.h:86
AffineFitness::maxAllowed
Matrix1D< double > maxAllowed
Definition: tomo_align_tilt_series.cpp:74
di
double * di
Definition: numerical_recipes.cpp:5168
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
AffineFitness::Powell_affine_fitness_individual
static double Powell_affine_fitness_individual(double *p, void *prm)
Definition: tomo_align_tilt_series.cpp:245
ProgTomographAlignment::refineChain
bool refineChain(LandmarkChain &chain, double &corrChain)
Definition: tomo_align_tilt_series.cpp:1997
bestNonwrappingShift
void bestNonwrappingShift(const MultidimArray< double > &I1, const MultidimArray< double > &I2, double &shiftX, double &shiftY, CorrelationAux &aux)
Definition: filters.cpp:1895
tomo_align_tilt_series.h
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
MDL_ENABLED
Is this image enabled? (int [-1 or 1])
Definition: metadata_label.h:231
selfScaleToSize
void selfScaleToSize(int SplineDegree, MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim)
Definition: transformations.cpp:736
MetaData::setComment
virtual void setComment(const String &newComment="No comment")
Definition: metadata_base.h:315
ProgTomographAlignment::computeAffineTransformations
void computeAffineTransformations(bool globalAffineToUse)
Compute affine transformations.
Definition: tomo_align_tilt_series.cpp:674
Alignment::computeGeometryDependentOfAxis
void computeGeometryDependentOfAxis()
Definition: tomo_align_tilt_series.cpp:2647
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
ERR_MD_OBJECTNUMBER
Incorrect number of objects in Metadata.
Definition: xmipp_error.h:160
ProgTomographAlignment::removeOutlierLandmarks
void removeOutlierLandmarks(const Alignment &alignment)
Remove Outliers.
Definition: tomo_align_tilt_series.cpp:2455
d
doublereal * d
Definition: numerical_recipes.cpp:2230
threadgenerateLandmarkSetCriticalPoints
void * threadgenerateLandmarkSetCriticalPoints(void *args)
Definition: tomo_align_tilt_series.cpp:1371
FOR_ALL_ELEMENTS_IN_ARRAY2D
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
Definition: multidim_array_base.h:390
Matrix1D::transpose
Matrix1D< T > transpose() const
Definition: matrix1d.cpp:644
ProgTomographAlignment::Ncritical
size_t Ncritical
Number of critical points.
Definition: tomo_align_tilt_series.h:98
ProgTomographAlignment::name_list
std::vector< std::string > name_list
Definition: tomo_align_tilt_series.h:249
rotation3DMatrix
void rotation3DMatrix(double ang, char axis, Matrix2D< double > &result, bool homogeneous)
Definition: transformations.cpp:426
ImageBase::tilt
double tilt(const size_t n=0) const
Definition: xmipp_image_base.cpp:378
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
CLIP
#define CLIP(x, x0, xF)
Definition: xmipp_macros.h:260
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
rnd_unif
double rnd_unif()
Definition: xmipp_funcs.cpp:440
ProgTomographAlignment::readParams
void readParams()
Read parameters from argument line.
Definition: tomo_align_tilt_series.cpp:452
FOR_ALL_ELEMENTS_IN_MATRIX1D
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
xmipp_image.h
M2
#define M2
Definition: numerical_recipes.cpp:49
b
doublereal * b
Definition: numerical_recipes.cpp:2230
M1
#define M1
Definition: numerical_recipes.cpp:45
ProgTomographAlignment::generateLandmarkSet
void generateLandmarkSet()
Generate landmark set using a grid.
Definition: tomo_align_tilt_series.cpp:1639
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
MultidimArray::computeMin
T computeMin() const
Definition: multidim_array.h:1568
y0
#define y0
AffineFitness::I2
std::vector< MultidimArray< double > * > I2
Definition: tomo_align_tilt_series.cpp:76
axis
char axis
Definition: image_operate.cpp:257
Alignment::diaxis
std::vector< Matrix1D< double > > diaxis
Definition: tomo_align_tilt_series.h:433
compute_hist
void compute_hist(const MultidimArrayGeneric &array, Histogram1D &hist, int no_steps)
Definition: histogram.cpp:572
x0
#define x0
XX
#define XX(v)
Definition: matrix1d.h:85
erode2D
void erode2D(const MultidimArray< double > &in, MultidimArray< double > &out, int neig, int count, int size)
Definition: morphology.cpp:116
MultidimArray::binarize
void binarize(double val=0, double accuracy=XMIPP_EQUAL_ACCURACY, MultidimArray< int > *mask=NULL)
Definition: multidim_array.h:3252
selfTranslate
void selfTranslate(int SplineDegree, MultidimArray< T > &V1, const Matrix1D< double > &v, bool wrap=xmipp_transformation::WRAP, T outside=0)
Definition: transformations.h:1394
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
ProgTomographAlignment::numThreads
int numThreads
Number of threads to use for parallel computing.
Definition: tomo_align_tilt_series.h:89
reject_outliers
void reject_outliers(T &v, double percentil_out=0.25)
Definition: histogram.h:605
Alignment::psi
Matrix1D< double > psi
Definition: tomo_align_tilt_series.h:301
fourier_filter.h
Alignment::computeGeometryDependentOfRotation
void computeGeometryDependentOfRotation()
Definition: tomo_align_tilt_series.cpp:2687
ProgTomographAlignment::alignImages
void alignImages(const Alignment &alignment)
Align images.
Definition: tomo_align_tilt_series.cpp:2245
XMIPP_EQUAL_ACCURACY
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
LandmarkChain
std::vector< Landmark > LandmarkChain
Definition: tomo_align_tilt_series.h:70
Matrix1D::resize
void resize(size_t Xdim, bool copy=true)
Definition: matrix1d.h:410
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
progress_bar
void progress_bar(long rlen)
Definition: xmipp_funcs.cpp:791
MultidimArray::write
void write(const FileName &fn) const
Definition: multidim_array.cpp:970
ProgTomographAlignment::thresholdAffine
double thresholdAffine
Threshold affine.
Definition: tomo_align_tilt_series.h:147
ProgTomographAlignment::show
void show()
Show parameters.
Definition: tomo_align_tilt_series.cpp:490
AffineFitness::Mask1
std::vector< MultidimArray< double > * > Mask1
Definition: tomo_align_tilt_series.cpp:77
threadComputeTransform
void * threadComputeTransform(void *args)
Definition: tomo_align_tilt_series.cpp:584
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
ProgTomographAlignment::readLandmarkSet
void readLandmarkSet(const FileName &fnLandmark)
Read landmark set.
Definition: tomo_align_tilt_series.cpp:2170
wrapperError
double wrapperError(double *p, void *prm)
Definition: tomo_align_tilt_series.cpp:2508
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
AffineFitness::affine_fitness_individual
double affine_fitness_individual(double *p)
Definition: tomo_align_tilt_series.cpp:100
Alignment::computeError
double computeError() const
Definition: tomo_align_tilt_series.cpp:2712
ERR_IO_NOTEXIST
File or directory does not exist.
Definition: xmipp_error.h:136
ROUND
#define ROUND(x)
Definition: xmipp_macros.h:210
ProgTomographAlignment::produceSideInfo
void produceSideInfo()
Produce side info.
Definition: tomo_align_tilt_series.cpp:742
vectorR2
Matrix1D< double > vectorR2(double x, double y)
Definition: matrix1d.cpp:884
computeAffineTransformation
double computeAffineTransformation(const MultidimArray< unsigned char > &I1, const MultidimArray< unsigned char > &I2, int maxShift, int maxIterDE, Matrix2D< double > &A12, Matrix2D< double > &A21, bool show, double thresholdAffine, bool globalAffine, bool isMirror, bool checkRotation, int pyramidLevel)
Definition: tomo_align_tilt_series.cpp:345
dummy
double dummy
Definition: numerical_recipes.cpp:4635
TomographAlignment::global_prm
const ProgTomographAlignment * global_prm
Definition: tomo_align_tilt_series.cpp:2505
MultidimArray::computeMedian
double computeMedian() const
Definition: multidim_array.h:1951
setupAffineFitness
void setupAffineFitness(AffineFitness &fitness, const MultidimArray< double > &I1, const MultidimArray< double > &I2, int maxShift, bool isMirror, bool checkRotation, int pyramidLevel)
Definition: tomo_align_tilt_series.cpp:251
AffineFitness::I1
std::vector< MultidimArray< double > * > I1
Definition: tomo_align_tilt_series.cpp:75
FileName::exists
bool exists() const
Definition: xmipp_filename.cpp:665
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
operator<<
std::ostream & operator<<(std::ostream &out, Alignment &alignment)
Definition: tomo_align_tilt_series.cpp:2879
MultidimArray::printStats
void printStats(std::ostream &out=std::cout) const
Definition: multidim_array.h:1507
ni
int ni
Definition: image_find_center.cpp:56
sort
void sort(struct DCEL_T *dcel)
Definition: sorting.cpp:18
ProgTomographAlignment::img
std::vector< MultidimArray< unsigned char > * > img
Definition: tomo_align_tilt_series.h:237
q2
#define q2
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
rotation2DMatrix
void rotation2DMatrix(T ang, Matrix2D< T > &result, bool homogeneous)
Definition: transformations.cpp:361
j
#define j
Definition: numerical_recipes.cpp:2493
steps
double steps
Definition: numerical_recipes.cpp:5166
YY
#define YY(v)
Definition: matrix1d.h:93
Alignment::raxis
Matrix1D< double > raxis
Definition: tomo_align_tilt_series.h:304
ImageBase::setEulerAngles
void setEulerAngles(double rot, double tilt, double psi, const size_t n=0)
Definition: xmipp_image_base.cpp:509
ProgTomographAlignment::affineTransformations
std::vector< std::vector< Matrix2D< double > > > affineTransformations
Definition: tomo_align_tilt_series.h:265
ProgTomographAlignment::localSize
double localSize
Local refinement size.
Definition: tomo_align_tilt_series.h:126
Alignment::optimizeGivenAxisDirection
double optimizeGivenAxisDirection()
Definition: tomo_align_tilt_series.cpp:2610
distance
TYPE distance(struct Point_T *p, struct Point_T *q)
Definition: point.cpp:28
error
void error(char *s)
Definition: tools.cpp:107
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
FOR_ALL_ELEMENTS_IN_ARRAY1D
#define FOR_ALL_ELEMENTS_IN_ARRAY1D(v)
Definition: multidim_array_base.h:554
memory_utils.h
MATRIX1D_ARRAY
#define MATRIX1D_ARRAY(v)
Definition: matrix1d.h:58
ERR_VALUE_NOTSET
Value has not been set.
Definition: xmipp_error.h:196
typeCast
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
RAD2DEG
#define RAD2DEG(r)
Definition: xmipp_macros.h:320
MAT_XSIZE
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
String
std::string String
Definition: xmipp_strings.h:34
alglib::psi
double psi(const double x)
Definition: specialfunctions.cpp:2492
ProgTomographAlignment::isCapillar
bool isCapillar
This tilt series comes from a capillar.
Definition: tomo_align_tilt_series.h:132
substractBackgroundRollingBall
void substractBackgroundRollingBall(MultidimArray< double > &I, int radius)
Definition: filters.cpp:75
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
ProgTomographAlignment::identifyOutliers
void identifyOutliers(bool mark)
Identify outliers.
Definition: tomo_align_tilt_series.cpp:701
ProgTomographAlignment::produceInformationFromLandmarks
void produceInformationFromLandmarks()
Produce information from landmarks.
Definition: tomo_align_tilt_series.cpp:2414
ProgTomographAlignment::writeLandmarkSet
void writeLandmarkSet(const FileName &fnLandmark) const
Write landmark set.
Definition: tomo_align_tilt_series.cpp:2146
bestShift
template float bestShift(MultidimArray< float > &, float &, float &, const MultidimArray< int > *, int)
ProgTomographAlignment::writeTransformations
void writeTransformations(const FileName &fnTransformations) const
Write affine transformations.
Definition: tomo_align_tilt_series.cpp:2210
threadgenerateLandmarkSetBlind
void * threadgenerateLandmarkSetBlind(void *args)
Definition: tomo_align_tilt_series.cpp:1216
OUTSIDE_MASK
constexpr int OUTSIDE_MASK
Definition: mask.h:48
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
K
constexpr int K
Definition: integration.cpp:123
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
pi
#define pi
Definition: spherical_harmonics.h:35
MDL_SHIFT_Y
Shift for the image in the Y axis (double)
Definition: metadata_label.h:435
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
MultidimArray::computeAvg
double computeAvg() const
Definition: multidim_array.h:1739
ProgTomographAlignment::showAffine
bool showAffine
Show affine transformations.
Definition: tomo_align_tilt_series.h:144
ERR_VALUE_INCORRECT
Incorrect value received.
Definition: xmipp_error.h:195
FourierFilter::generateMask
void generateMask(MultidimArray< double > &v)
Definition: fourier_filter.cpp:605
Matrix1D::initConstant
void initConstant(T val)
Definition: matrix1d.cpp:83
substractBackgroundPlane
void substractBackgroundPlane(MultidimArray< double > &I)
Definition: filters.cpp:40
n
int * n
Definition: numerical_recipes.cpp:2229
MDL_IMAGE
Name of an image (std::string)
Definition: metadata_label.h:246
fitness
double fitness(double *p)
Definition: volume_align_prog.cpp:88
ProgTomographAlignment::correlationList
std::vector< double > correlationList
Definition: tomo_align_tilt_series.h:268
Matrix2D::resize
void resize(size_t Ydim, size_t Xdim, bool noCopy=false)
Definition: matrix2d.cpp:1022
ZZ
#define ZZ(v)
Definition: matrix1d.h:101
a
doublereal * a
Definition: numerical_recipes.cpp:2230
intWRAP
#define intWRAP(x, x0, xF)
Definition: xmipp_macros.h:272
MultidimArray::indexSort
void indexSort(MultidimArray< int > &indx) const
Definition: multidim_array.cpp:882
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710
xmipp_transformation::BSPLINE3
Definition: transformations_defines.h:32
Generated by   doxygen 1.8.13