Xmipp  v3.23.11-Nereus
filters.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Carlos Oscar S. Sorzano (coss@cnb.csic.es)
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 <queue>
27 #include <list>
28 #include "data/fourier_filter.h"
29 #include "core/histogram.h"
30 #include "core/metadata_extension.h"
31 #include "core/multidim_array.h"
32 #include "core/transformations.h"
33 #include "core/xmipp_program.h"
34 #include "filters.h"
35 #include "mask.h"
36 #include "morphology.h"
37 #include "wavelet.h"
38 
39 /* Subtract background ---------------------------------------------------- */
40 void substractBackgroundPlane(MultidimArray<double> &I)
41 {
42 
43  I.checkDimension(2);
44 
45  Matrix2D<double> A(3, 3);
46  Matrix1D<double> x(3);
47  Matrix1D<double> b(3);
48 
49  // Solve the plane 'x'
50  A.initZeros();
51  b.initZeros();
52  FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
53  {
54  A(0, 0) += j * j;
55  A(0, 1) += j * i;
56  A(0, 2) += j;
57  A(1, 1) += i * i;
58  A(1, 2) += i;
59  A(2, 2) += 1;
60  b(0) += j * A2D_ELEM(I, i, j);
61  b(1) += i * A2D_ELEM(I, i, j);
62  b(2) += A2D_ELEM(I, i, j);
63  }
64  A(1, 0) = A(0, 1);
65  A(2, 0) = A(0, 2);
66  A(2, 1) = A(1, 2);
67  solve(A, b, x);
68 
69  // Now subtract the plane
70  FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
71  A2D_ELEM(I, i, j) -= x(0) * i + x(1) * j + x(2);
72 }
73 
74 /* Subtract background ---------------------------------------------------- */
75 void substractBackgroundRollingBall(MultidimArray<double> &I, int radius)
76 {
77 
78  I.checkDimension(2);
79 
80  // Build the ball
81  int arcTrimPer;
82  int shrinkFactor;
83  if (radius <= 10)
84  {
85  shrinkFactor = 1;
86  arcTrimPer = 24; // trim 24% in x and y
87  }
88  else if (radius <= 30)
89  {
90  shrinkFactor = 2;
91  arcTrimPer = 24; // trim 24% in x and y
92  }
93  else if (radius <= 100)
94  {
95  shrinkFactor = 4;
96  arcTrimPer = 32; // trim 32% in x and y
97  }
98  else
99  {
100  shrinkFactor = 8;
101  arcTrimPer = 40; // trim 40% in x and y
102  }
103 
104  double smallballradius = radius / shrinkFactor;
105  if (smallballradius < 1)
106  smallballradius = 1;
107  double r2 = smallballradius * smallballradius;
108  int xtrim = (int) (arcTrimPer * smallballradius) / 100; // only use a patch of the rolling ball
109  int halfWidth = ROUND(smallballradius - xtrim);
110  int ballWidth = 2 * halfWidth + 1;
111  MultidimArray<double> ball(ballWidth, ballWidth);
112  ball.setXmippOrigin();
113  FOR_ALL_ELEMENTS_IN_ARRAY2D(ball)
114  {
115  double temp = r2 - i * i - j * j;
116  ball(i, j) = temp > 0. ? sqrt(temp) : 0.;
117  }
118 
119  // Shrink the image: each point in the shrinked image is the
120  // minimum of its neighbourhood
121  int sXdim = (XSIZE(I) + shrinkFactor - 1) / shrinkFactor;
122  int sYdim = (YSIZE(I) + shrinkFactor - 1) / shrinkFactor;
123  MultidimArray<double> shrinkI(sYdim, sXdim);
124  shrinkI.setXmippOrigin();
125  for (int ySmall = 0; ySmall < sYdim; ySmall++)
126  {
127  for (int xSmall = 0; xSmall < sXdim; xSmall++)
128  {
129  double minVal = 1e38;
130  for (int j = 0, y = shrinkFactor * ySmall;
131  j < shrinkFactor && y < (int)YSIZE(I); j++, y++)
132  for (int k = 0, x = shrinkFactor * xSmall;
133  k < shrinkFactor && x < (int)XSIZE(I); k++, x++)
134  {
135  double thispixel = DIRECT_A2D_ELEM(I,y,x);
136  if (thispixel < minVal)
137  minVal = thispixel;
138  }
139  DIRECT_A2D_ELEM(shrinkI,ySmall,xSmall) = minVal;
140  }
141  }
142 
143  // Now roll the ball
144  radius = ballWidth / 2;
145  MultidimArray<double> Irolled;
146  Irolled.resizeNoCopy(shrinkI);
147  Irolled.initConstant(-500);
148  for (int yb = -radius; yb < (int)YSIZE(shrinkI) + radius; yb++)
149  {
150  // Limits of the ball
151  int y0 = yb - radius;
152  if (y0 < 0)
153  y0 = 0;
154  int y0b = y0 - yb + radius; //y coordinate in the ball corresponding to y0
155  int yF = yb + radius;
156  if (yF >= (int)YSIZE(shrinkI))
157  yF = (int)YSIZE(shrinkI) - 1;
158 
159  for (int xb = -radius; xb < (int)XSIZE(shrinkI) + radius; xb++)
160  {
161  // Limits of the ball
162  int x0 = xb - radius;
163  if (x0 < 0)
164  x0 = 0;
165  int x0b = x0 - xb + radius;
166  int xF = xb + radius;
167  if (xF >= (int)XSIZE(shrinkI))
168  xF = (int)XSIZE(shrinkI) - 1;
169 
170  double z = 1e38;
171  for (int yp = y0, ybp = y0b; yp <= yF; yp++, ybp++)
172  for (int xp = x0, xbp = x0b; xp <= xF; xp++, xbp++)
173  {
174  double zReduced = DIRECT_A2D_ELEM(shrinkI,yp,xp)
175  - DIRECT_A2D_ELEM(ball,ybp,xbp);
176  if (z > zReduced)
177  z = zReduced;
178  }
179  for (int yp = y0, ybp = y0b; yp <= yF; yp++, ybp++)
180  for (int xp = x0, xbp = x0b; xp <= xF; xp++, xbp++)
181  {
182  double zMin = z + DIRECT_A2D_ELEM(ball,ybp,xbp);
183  if (DIRECT_A2D_ELEM(Irolled,yp,xp) < zMin)
184  DIRECT_A2D_ELEM(Irolled,yp,xp) = zMin;
185  }
186  }
187  }
188 
189  // Now rescale the background
190  MultidimArray<double> bgEnlarged;
191  scaleToSize(xmipp_transformation::LINEAR, bgEnlarged, Irolled, XSIZE(I), YSIZE(I));
192  bgEnlarged.copyShape(I);
193  I -= bgEnlarged;
194 }
195 
196 /* Detect background ------------------------------------------------------ */
197 void detectBackground(const MultidimArray<double> &vol,
198  MultidimArray<double> &mask, double alpha, double &final_mean)
199 {
200 
201  // 2.1.-Background detection------------------------------------------
202  MultidimArray<double> bg; // We create the volumen with
203  bg.resizeNoCopy(vol); // -1:not visited 0:mol 1:background
204  bg.initConstant(-1); // -2:in the list
205 
206  // Ponemos las seis caras de esta variable como visitadas e inicializamos
207  // la cola de pixeles por visitar
208  std::queue<int> list_for_compute; // Lista del modo [x1,y1,z1,...,xi,yi,zi]
209  // que contiene los pixeles por procesar
210  std::vector<double> bg_values; // Vector con los valores del background
211  size_t xdim = XSIZE(bg);
212  size_t ydim = YSIZE(bg);
213  size_t zdim = ZSIZE(bg);
214  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(bg)
215  {
216  if (j == 0 || j == xdim - 1 || i == 0 || i == ydim - 1 || k == 0
217  || k == zdim - 1)
218  { // Visited coord
219  DIRECT_A3D_ELEM(bg,k,i,j) = 1;
220  // We introduce the values of the background
221  bg_values.push_back(DIRECT_A3D_ELEM(vol,k,i,j));
222 
223  }
224  if ((j == 1 || j == xdim - 2) && (i != 0) && (i != ydim - 1) && (k != 0)
225  && (k != zdim - 1))
226  { // Coord for compute for x
227  list_for_compute.push(j);
228  list_for_compute.push(i);
229  list_for_compute.push(k);
230  }
231  if ((i == 1 || i == ydim - 2) && (j > 1) && (j < xdim - 2) && (k != 0)
232  && (k != zdim - 1))
233  { // Coord for compute for y
234  list_for_compute.push(j);
235  list_for_compute.push(i);
236  list_for_compute.push(k);
237  }
238  if ((k == 1 || k == zdim - 2) && (j > 1) && (j < xdim - 2) && (i > 1)
239  && (i < ydim - 2))
240  { // Coord for compute for y
241  list_for_compute.push(j);
242  list_for_compute.push(i);
243  list_for_compute.push(k);
244  }
245  } // end of FOR_ALL_ELEMENTS
246 
247  // We work until the list_for_compute is empty
248  int n = 250; //each 250 pixels renew stats
249  int cont = 250; //We start here for compute stat for first time
250  double A=0; // A and B are numbers such the interval of confidence is [A,B]
251  double B=0; //
252  float z = icdf_gauss(1 - alpha / 2);
253  while (!list_for_compute.empty())
254  {
255 
256  //We compute stat when is needed
257  if (cont == n)
258  {
259  // Compute statistics
260  double avg=0;
261  double stddev=0;
262  computeAvgStddev(bg_values, avg, stddev);
263  final_mean = avg;
264  // Compute confidence interval
265  A = avg - (z * stddev);
266  B = avg + (z * stddev);
267  cont = 0;
268  } // end of if
269  // Now we start to take coords from the list_for_compute
270  int x_coord = list_for_compute.front();
271  list_for_compute.pop();
272  int y_coord = list_for_compute.front();
273  list_for_compute.pop();
274  int z_coord = list_for_compute.front();
275  list_for_compute.pop();
276  // Is visited
277  DIRECT_A3D_ELEM(bg,z_coord,y_coord,x_coord) = -2;
278  //We take the value
279  double value = DIRECT_A3D_ELEM(vol,z_coord,y_coord,x_coord);
280  // We see if is background or not
281  if (A <= value && value <= B)
282  {
283  // We now is background
284  DIRECT_A3D_ELEM(bg,z_coord,y_coord,x_coord) = 1;
285  // We introduce the values of the background
286  bg_values.push_back(value);
287  // We update the cont variable
288  cont++;
289 
290  // We add his neighbours in the list_for_compute
291  for (int xx = x_coord - 1; xx <= x_coord + 1; xx++)
292  for (int yy = y_coord - 1; yy <= y_coord + 1; yy++)
293  for (int zz = z_coord - 1; zz <= z_coord + 1; zz++)
294  //We see if it has been visited
295  if (DIRECT_A3D_ELEM(bg,zz,yy,xx) == -1) // not visited
296  {
297  // So we include it in the list
298  list_for_compute.push(xx);
299  list_for_compute.push(yy);
300  list_for_compute.push(zz);
301  // Is in the list
302  DIRECT_A3D_ELEM(bg,zz,yy,xx) = -2;
303  }
304  } // end of if
305  else
306  {
307  // Isn't background
308  DIRECT_A3D_ELEM(bg,z_coord,y_coord,x_coord) = 0;
309  }
310  } // end of while
311  // Now we change not visited for mol
312  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(bg)
313  if (DIRECT_A3D_ELEM(bg,k,i,j) == -1)
314  DIRECT_A3D_ELEM(bg,k,i,j) = 0;
315 
316  // End of 2.1-----------------------------------------------------------
317  // 2.2.-Matematical Morphology
318  MultidimArray<double> bg_mm; // We create the output volumen
319  bg_mm.initZeros(vol);
320  closing3D(bg, bg_mm, 26, 0, 1);
321  // Output
322  //typeCast(bg_mm,mask);
323  mask = bg_mm;
324 }
325 
326 /* Contranst enhancement --------------------------------------------------- */
327 void contrastEnhancement(Image<double> *I)
328 {
329  (*I)().rangeAdjust(0, 255);
330 }
331 
332 /* Region growing for images ----------------------------------------------- */
333 typedef struct
334 {
335  int ii;
336  int jj;
337 }
338 Coordinate2D;
339 
340 void regionGrowing2D(const MultidimArray<double> &I_in,
341  MultidimArray<double> &I_out, int i, int j, float stop_colour,
342  float filling_colour, bool less, int neighbourhood)
343 {
344  I_in.checkDimension(2);
345 
346  std::list<Coordinate2D> iNeighbours; /* A list for neighbour pixels */
347  int iCurrenti;
348  int iCurrentj; /* Coordinates of the current pixel considered */
349 
350  /* First task is copying the input image into the output one */
351  I_out = I_in;
352 
353  /**** Then the region growing is done ****/
354  /* Insert at the beginning of the list the seed coordinates */
355  Coordinate2D coord;
356  coord.ii = i;
357  coord.jj = j;
358  iNeighbours.push_front(coord);
359 
360  /* Fill the seed coordinates */
361  A2D_ELEM(I_out, i, j) = filling_colour;
362 
363  while (!iNeighbours.empty())
364  {
365  /* Take the current pixel to explore */
366  coord = iNeighbours.front();
367  iNeighbours.pop_front();
368  iCurrenti = coord.ii;
369  iCurrentj = coord.jj;
370 
371 #define CHECK_POINT(i,j) \
372  { \
373  int I=i; \
374  int J=j; \
375  if (INSIDEXY(I_out,J,I)) { \
376  double pixel=A2D_ELEM(I_out,I,J);\
377  if (pixel!=filling_colour) \
378  if ((less && pixel < stop_colour) || \
379  (!less && pixel > stop_colour)) { \
380  coord.ii=I; \
381  coord.jj=J; \
382  A2D_ELEM (I_out,coord.ii,coord.jj)=filling_colour; \
383  iNeighbours.push_front(coord); \
384  } \
385  } \
386  }
387 
388  /* Make the exploration of the pixel's neighbours */
389  CHECK_POINT(iCurrenti, iCurrentj - 1);
390  CHECK_POINT(iCurrenti, iCurrentj + 1);
391  CHECK_POINT(iCurrenti - 1, iCurrentj);
392  CHECK_POINT(iCurrenti + 1, iCurrentj);
393  if (neighbourhood == 8)
394  {
395  CHECK_POINT(iCurrenti - 1, iCurrentj - 1);
396  CHECK_POINT(iCurrenti - 1, iCurrentj + 1);
397  CHECK_POINT(iCurrenti + 1, iCurrentj - 1);
398  CHECK_POINT(iCurrenti + 1, iCurrentj + 1);
399  }
400  }
401 }
402 
403 /* Region growing for volumes ----------------------------------------------- */
404 typedef struct
405 {
406  int ii;
407  int jj;
408  int kk;
409 }
410 Coordinate3D;
411 
412 void regionGrowing3D(const MultidimArray<double> &V_in,
413  MultidimArray<double> &V_out, int k, int i, int j, float stop_colour,
414  float filling_colour, bool less)
415 {
416  V_in.checkDimension(3);
417 
418  std::list<Coordinate3D> iNeighbours; /* A list for neighbour voxels */
419  int iCurrentk;
420  int iCurrenti;
421  int iCurrentj; /* Coordinates of the current voxel considered */
422 
423  /* First task is copying the input volume into the output one */
424  V_out = V_in;
425 
426  /**** Then the region growing is done in output volume ****/
427  /* Insert at the beginning of the list the seed coordinates */
428  Coordinate3D coord;
429  coord.ii = i;
430  coord.jj = j;
431  coord.kk = k;
432  iNeighbours.push_front(coord);
433 
434  /* Fill the seed coordinates */
435  A3D_ELEM(V_out, k, i, j) = filling_colour;
436 
437  while (!iNeighbours.empty())
438  {
439  /* Take the current pixel to explore */
440  coord = iNeighbours.front();
441  iNeighbours.pop_front();
442  iCurrenti = coord.ii;
443  iCurrentj = coord.jj;
444  iCurrentk = coord.kk;
445 
446  /* a macro for doing exploration of a voxel. If the voxel has a value
447  lower than stop_colour, its filled with filling colour and added to the
448  list for exploring its neighbours */
449 #define CHECK_POINT_3D(k,i,j) \
450  {\
451  int I=i; \
452  int J=j; \
453  int K=k; \
454  if (INSIDEXYZ(V_out,J,I,K)) { \
455  double voxel=A3D_ELEM(V_out,K,I,J); \
456  if (voxel!=filling_colour) \
457  if ((less && voxel < stop_colour)|| \
458  (!less &&voxel > stop_colour)) { \
459  coord.ii=I; \
460  coord.jj=J; \
461  coord.kk=K; \
462  A3D_ELEM (V_out,coord.kk,coord.ii,coord.jj)=filling_colour; \
463  iNeighbours.push_front(coord); \
464  } \
465  }\
466  }
467 
468  /* Make the exploration of the pixel neighbours */
469  CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj - 1);
470  CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj + 1);
471  CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj);
472  CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj - 1);
473  CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj + 1);
474  CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj);
475  CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj - 1);
476  CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj + 1);
477  CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj);
478  CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj - 1);
479  CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj + 1);
480  CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj);
481  CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj - 1);
482  CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj + 1);
483  CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj);
484  CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj - 1);
485  CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj + 1);
486  CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj);
487  CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj - 1);
488  CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj + 1);
489  CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
490  CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj - 1);
491  CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
492  CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
493  CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj - 1);
494  CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
495  }
496 }
497 
498 
499 void regionGrowing3DEqualValue(const MultidimArray<double> &V_in, MultidimArray<int> &V_out, int filling_value=0)
500 {
501  V_in.checkDimension(3);
502 
503  std::list<Coordinate3D> iNeighbours; /* A list for neighbour voxels */
504  int iCurrentk;
505  int iCurrenti;
506  int iCurrentj; /* Coordinates of the current voxel considered */
507 
508  /* First task is copying the input volume into the output one */
509  V_out.resizeNoCopy(V_in);
510  V_out.initConstant(1);
511 
512  /**** Then the region growing is done in output volume ****/
513  /* Insert at the beginning of the list the seed coordinates */
514  Coordinate3D coord;
515  coord.ii = STARTINGY(V_in);
516  coord.jj = STARTINGX(V_in);
517  coord.kk = STARTINGZ(V_in);
518  iNeighbours.push_front(coord);
519 
520  /* Fill the seed coordinates */
521  double bgColor = A3D_ELEM(V_in, STARTINGZ(V_in), STARTINGY(V_in), STARTINGX(V_in));
522  A3D_ELEM(V_out, STARTINGZ(V_out), STARTINGY(V_out), STARTINGX(V_out)) = 0;
523 
524  while (!iNeighbours.empty())
525  {
526  /* Take the current pixel to explore */
527  coord = iNeighbours.front();
528  iNeighbours.pop_front();
529  iCurrenti = coord.ii;
530  iCurrentj = coord.jj;
531  iCurrentk = coord.kk;
532 
533  /* a macro for doing exploration of a voxel. If the voxel has a value
534  lower than stop_colour, its filled with filling colour and added to the
535  list for exploring its neighbours */
536 #undef CHECK_POINT_3D
537 #define CHECK_POINT_3D(k,i,j) \
538  {\
539  int I=i; \
540  int J=j; \
541  int K=k; \
542  if (INSIDEXYZ(V_out,J,I,K)) { \
543  if (A3D_ELEM(V_in,K,I,J)==bgColor && A3D_ELEM(V_out,K,I,J)==1) { \
544  coord.ii=I; \
545  coord.jj=J; \
546  coord.kk=K; \
547  A3D_ELEM (V_out,coord.kk,coord.ii,coord.jj)=filling_value; \
548  iNeighbours.push_front(coord); \
549  } \
550  }\
551  }\
552 
553  /* Make the exploration of the pixel neighbours */
554  CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj - 1);
555  CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj + 1);
556  CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj);
557  CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj - 1);
558  CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj + 1);
559  CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj);
560  CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj - 1);
561  CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj + 1);
562  CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj);
563  CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj - 1);
564  CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj + 1);
565  CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj);
566  CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj - 1);
567  CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj + 1);
568  CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj);
569  CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj - 1);
570  CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj + 1);
571  CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj);
572  CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj - 1);
573  CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj + 1);
574  CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
575  CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj - 1);
576  CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
577  CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
578  CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj - 1);
579  CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
580  }
581 }
582 
583 
584 void distanceTransform(const MultidimArray<int> &in, MultidimArray<int> &out,
585  bool wrap)
586 {
587  std::list<int> toExplore; /* A list of points to explore */
588 
589  in.checkDimension(2);
590 
591  out.resize(in);
592  out.initConstant(XSIZE(in) + YSIZE(in));
593 
594 #define CHECK_POINT_DIST(i,j,proposedDistance) \
595  { \
596  int ip=i; \
597  int jp=j; \
598  if (wrap) \
599  { \
600  ip=intWRAP(ip,STARTINGY(out),FINISHINGY(out));\
601  jp=intWRAP(jp,STARTINGX(out),FINISHINGX(out));\
602  } \
603  if (ip>=STARTINGY(out) && ip<=FINISHINGY(out) && \
604  jp>=STARTINGX(out) && jp<=FINISHINGX(out)) \
605  if (out(ip,jp)>proposedDistance) \
606  { \
607  out(ip,jp)=proposedDistance; \
608  toExplore.push_back(ip); \
609  toExplore.push_back(jp); \
610  toExplore.push_back(proposedDistance); \
611  } \
612  }
613 
614  // Look for all elements in the binary mask and set the corresponding
615  // distance to 0
616  FOR_ALL_ELEMENTS_IN_ARRAY2D(in)
617  if (in(i, j))
618  {
619  out(i, j) = 0;
620  CHECK_POINT_DIST(i-1, j, 1);
621  CHECK_POINT_DIST(i+1, j, 1);
622  CHECK_POINT_DIST(i, j-1, 1);
623  CHECK_POINT_DIST(i, j+1, 1);
624  }
625 
626  while (!toExplore.empty())
627  {
628  int i = toExplore.front();
629  toExplore.pop_front();
630  int j = toExplore.front();
631  toExplore.pop_front();
632  int proposedDistance = toExplore.front();
633  toExplore.pop_front();
634 
635  if (proposedDistance == out(i, j))
636  {
637  // If this is the current distance (i.e., it has not
638  // been supersceded by a later value
639  CHECK_POINT_DIST(i-1, j, proposedDistance+1);
640  CHECK_POINT_DIST(i+1, j, proposedDistance+1);
641  CHECK_POINT_DIST(i, j-1, proposedDistance+1);
642  CHECK_POINT_DIST(i, j+1, proposedDistance+1);
643  }
644  }
645 }
646 
647 /* Label image ------------------------------------------------------------ */
648 int labelImage2D(const MultidimArray<double> &I, MultidimArray<double> &label,
649  int neighbourhood)
650 {
651  I.checkDimension(2);
652 
653  label = I;
654  int colour = 32000;
655  FOR_ALL_ELEMENTS_IN_ARRAY2D(label)
656  {
657  if (label(i, j) != 1)
658  continue;
659  regionGrowing2D(label, label, i, j, 0, colour, false, neighbourhood);
660  colour++;
661  }
662  FOR_ALL_ELEMENTS_IN_ARRAY2D(label)
663  if (label(i, j) != 0)
664  label(i, j) = label(i, j) - 31999;
665  return colour - 32000;
666 }
667 
668 /* Label volume ------------------------------------------------------------ */
669 int labelImage3D(const MultidimArray<double> &V, MultidimArray<double> &label)
670 {
671  V.checkDimension(3);
672 
673  label = V;
674  int colour = 32000;
675  FOR_ALL_ELEMENTS_IN_ARRAY3D(label)
676  {
677  if (label(k, i, j) != 1)
678  continue;
679  regionGrowing3D(label, label, k, i, j, 0, colour, false);
680  colour++;
681  }
682  FOR_ALL_ELEMENTS_IN_ARRAY3D(label)
683  if (A3D_ELEM(label,k, i, j) != 0)
684  A3D_ELEM(label,k, i, j) = A3D_ELEM(label,k, i, j) - 31999;
685  return colour - 32000;
686 }
687 
688 /* Remove small components ------------------------------------------------- */
689 void removeSmallComponents(MultidimArray<double> &I, int size,
690  int neighbourhood)
691 {
692  MultidimArray<double> label;
693  int imax;
694  if (ZSIZE(I)==1)
695  imax=labelImage2D(I, label, neighbourhood);
696  else
697  imax=labelImage3D(I, label);
698  MultidimArray<int> nlabel(imax + 1);
699  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
700  {
701  int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
702  DIRECT_A1D_ELEM(nlabel,l)++;
703  }
704  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
705  {
706  int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
707  if (DIRECT_A1D_ELEM(nlabel,l)<size)
708  DIRECT_MULTIDIM_ELEM(I,n)=0;
709  }
710 }
711 
712 /* Keep biggest component -------------------------------------------------- */
713 void keepBiggestComponent(MultidimArray<double> &I, double percentage,
714  int neighbourhood)
715 {
716  MultidimArray<double> label;
717  int imax;
718  if (ZSIZE(I)==1)
719  imax=labelImage2D(I, label, neighbourhood);
720  else
721  imax=labelImage3D(I, label);
722  MultidimArray<int> nlabel(imax + 1);
723  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
724  {
725  int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
726  if (l>0)
727  A1D_ELEM(nlabel,l)++;
728  }
729 
730  MultidimArray <int> best;
731  nlabel.indexSort(best);
732  best -= 1;
733  int nbest = XSIZE(best) - 1;
734  double total = nlabel.sum();
735  double explained = nlabel(best(nbest));
736  while (explained < percentage * total)
737  {
738  nbest--;
739  explained += nlabel(best(nbest));
740  }
741 
742  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
743  {
744  int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
745  bool among_the_best = false;
746  for (int k = nbest; k < imax + 1; k++)
747  if (l == A1D_ELEM(best,k))
748  {
749  among_the_best = true;
750  break;
751  }
752  if (!among_the_best)
753  DIRECT_MULTIDIM_ELEM(I,n)=0;
754  }
755 }
756 
757 /* Fill object ------------------------------------------------------------- */
758 void fillBinaryObject(MultidimArray<double> &I, int neighbourhood)
759 {
760  I.checkDimension(2);
761 
762  MultidimArray<double> label;
763  FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
764  I(i, j) = 1 - I(i, j);
765  labelImage2D(I, label, neighbourhood);
766  double l0 = label(STARTINGY(I), STARTINGX(I));
767  FOR_ALL_ELEMENTS_IN_ARRAY2D(label)
768  if (label(i, j) == l0)
769  I(i, j) = 0;
770  else
771  I(i, j) = 1;
772 }
773 
774 /* Variance filter ----------------------------------------------------------*/
775 void varianceFilter(MultidimArray<double> &I, int kernelSize, bool relative)
776 {
777  int kernelSize_2 = kernelSize/2;
778  MultidimArray<double> kernel;
779  kernel.resize(kernelSize,kernelSize);
780  kernel.setXmippOrigin();
781 
782  // std::cout << " Creating the variance matrix " << std::endl;
783  MultidimArray<double> mVar(YSIZE(I),XSIZE(I));
784  mVar.setXmippOrigin();
785  double stdKernel;
786  double avgKernel;
787  double min_val;
788  double max_val;
789  int x0;
790  int y0;
791  int xF;
792  int yF;
793 
794  // I.computeStats(avgImg, stdImg, min_im, max_im);
795 
796  for (int i=kernelSize_2; i<=(int)YSIZE(I)-kernelSize_2; i+=kernelSize_2)
797  for (int j=kernelSize_2; j<=(int)XSIZE(I)-kernelSize_2; j+=kernelSize_2)
798  {
799  x0 = j-kernelSize_2;
800  y0 = i-kernelSize_2;
801  xF = j+kernelSize_2-1;
802  yF = i+kernelSize_2-1;
803 
804  if (x0 < 0)
805  x0 = 0;
806  if (xF > XSIZE(I))
807  xF = XSIZE(I);
808  if (y0 < 0)
809  y0 = 0;
810  if (yF > YSIZE(I))
811  yF = YSIZE(I);
812 
813  I.window(kernel, y0, x0, yF, xF);
814  kernel.computeStats(avgKernel, stdKernel, min_val, max_val);
815 
816  DIRECT_A2D_ELEM(mVar, i, j) = stdKernel;
817  }
818 
819  // filtering to fill the matrices (convolving with a Gaussian)
820  FourierFilter filter;
821  filter.FilterShape = REALGAUSSIAN;
822  filter.FilterBand = LOWPASS;
823  filter.w1 = kernelSize_2;
824  filter.applyMaskSpace(mVar);
825 
826  if (relative) // normalize to the mean of the variance
827  {
828  double avgVar;
829  double stdVar;
830  double minVar;
831  double maxVar;
832  mVar.computeStats(avgVar, stdVar, minVar, maxVar);
833  mVar = mVar/avgVar;
834  }
835 
836  I = mVar;
837 
838  // filling the borders with the nearest variance value
839  // the corners only are filled once (with Ycoord)
840  for (int i=0; i<YSIZE(I); ++i)
841  for (int j=0; j<kernelSize; j++)
842  {
843  if (i<kernelSize)
844  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, kernelSize, kernelSize);
845  else if (i>YSIZE(I)-kernelSize)
846  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, YSIZE(I)-kernelSize, kernelSize);
847  else
848  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, i, kernelSize);
849  }
850  for (int i=0; i<YSIZE(I); ++i)
851  for (int j=XSIZE(I)-kernelSize; j<XSIZE(I); j++)
852  {
853  if (i<kernelSize)
854  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, kernelSize, XSIZE(I)-kernelSize);
855  else if (i>YSIZE(I)-kernelSize)
856  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, YSIZE(I)-kernelSize, XSIZE(I)-kernelSize);
857  else
858  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, i, XSIZE(I)-kernelSize);
859  }
860  for (int j=kernelSize; j<XSIZE(I)-kernelSize; ++j)
861  for (int i=0; i<kernelSize; i++)
862  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, kernelSize, j);
863  for (int j=kernelSize; j<XSIZE(I)-kernelSize; ++j)
864  for (int i=YSIZE(I)-kernelSize; i<YSIZE(I); i++)
865  DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, YSIZE(I)-kernelSize, j);
866 
867 }
868 
869 /* Noise filter (returns a binary mask where both variance and mean are high)*/
870 void noisyZonesFilter(MultidimArray<double> &I, int kernelSize)
871 {
872  int kernelSize_2 = kernelSize/2;
873  MultidimArray<double> kernel;
874  kernel.resize(kernelSize,kernelSize);
875  kernel.setXmippOrigin();
876 
877  MultidimArray<double> mAvg=I;
878  MultidimArray<double> mVar=I;
879  double stdKernel;
880  double varKernel;
881  double avgKernel;
882  double min_val;
883  double max_val;
884  int x0;
885  int y0;
886  int xF;
887  int yF;
888 
889  for (int i=kernelSize_2; i<(int)YSIZE(I); i+=kernelSize_2)
890  for (int j=kernelSize_2; j<(int)XSIZE(I); j+=kernelSize_2)
891  {
892  x0 = j-kernelSize_2;
893  y0 = i-kernelSize_2;
894  xF = j+kernelSize_2-1;
895  yF = i+kernelSize_2-1;
896 
897  if (x0 < 0)
898  x0 = 0;
899  if (xF > XSIZE(I))
900  xF = XSIZE(I);
901  if (y0 < 0)
902  y0 = 0;
903  if (yF > YSIZE(I))
904  yF = YSIZE(I);
905 
906  I.window(kernel, y0, x0, yF, xF);
907  kernel.computeStats(avgKernel, stdKernel, min_val, max_val);
908  varKernel = stdKernel*stdKernel;
909 
910  DIRECT_A2D_ELEM(mAvg, i, j) = avgKernel*avgKernel;
911  DIRECT_A2D_ELEM(mVar, i, j) = varKernel;
912  }
913 
914  // filtering to fill the matrices (convolving with a Gaussian)
915  FourierFilter filter;
916  filter.FilterShape = REALGAUSSIAN;
917  filter.FilterBand = LOWPASS;
918  filter.w1 = kernelSize_2;
919  filter.applyMaskSpace(mAvg);
920  filter.applyMaskSpace(mVar);
921 
922  // // Draw to debug
923  // Image<double> imAvg(mAvg), imVar(mVar);
924  // imAvg.write("AvgFilter.mrc");
925  // imVar.write("VarFilter.mrc");
926 
927  // Working in a auxilary windows to avoid borders bad defined
928  auto ysize = static_cast<int>YSIZE(I);
929  auto xsize = static_cast<int>XSIZE(I);
930  MultidimArray<double> mAvgAux(ysize-kernelSize,xsize-kernelSize);
931  MultidimArray<double> mVarAux(ysize-kernelSize,xsize-kernelSize);
932  mAvgAux.setXmippOrigin();
933  mVarAux.setXmippOrigin();
934  mAvg.window(mAvgAux,STARTINGY(mAvg)+kernelSize_2, STARTINGX(mAvg)+kernelSize_2,
935  FINISHINGY(mAvg)-kernelSize_2, FINISHINGX(mAvg)-kernelSize_2);
936  mVar.window(mVarAux,STARTINGY(mVar)+kernelSize_2, STARTINGX(mVar)+kernelSize_2,
937  FINISHINGY(mVar)-kernelSize_2, FINISHINGX(mVar)-kernelSize_2);
938 
939  // Refiltering to get a smoother distribution
940  // filter.w1 = XSIZE(I)/40;
941  // filter.applyMaskSpace(mAvgAux);
942  // filter.applyMaskSpace(mVarAux);
943 
944  // Binarization
945  MultidimArray<double> mAvgAuxBin = mAvgAux;
946  MultidimArray<double> mVarAuxBin = mVarAux;
947  EntropySegmentation(mVarAuxBin);
948  // EntropySegmentation(mAvgAuxBin);
949  float th = EntropySegmentation(mAvgAuxBin);
950  mAvgAuxBin.binarize(th*0.92);
951  mAvgAuxBin = 1-mAvgAuxBin;
952  // std::cout << "binarize threshold = " << th << std::endl;
953 
954  // Returning to the previous windows size
955  MultidimArray<double> mAvgBin = mAvg;
956  MultidimArray<double> mVarBin = mVar;
957  mAvgAuxBin.window(mAvgBin, STARTINGY(mVar), STARTINGX(mVar),
958  FINISHINGY(mVar), FINISHINGX(mVar));
959  mVarAuxBin.window(mVarBin, STARTINGY(mVar), STARTINGX(mVar),
960  FINISHINGY(mVar), FINISHINGX(mVar));
961 
962  // // Draw to debug
963  // Image<double> imAvgBin(mAvgBin), imVarBin(mVarBin);
964  // imAvgBin.write("noisyZoneFilter_AVGmask.mrc");
965  // imVarBin.write("noisyZoneFilter_VARmask.mrc");
966 
967  // Combining both masks
968  I = 1-(mVarBin*mAvgBin);
969 }
970 
971 /* Gini Coefficient -- (applies a variance filter to the input Image) ------ */
972 double giniCoeff(MultidimArray<double> &I, int varKernelSize)
973 {
974  MultidimArray<double> im = I;
975 
976  // std::cout << " - Starting fft filtering " << std::endl;
977 
978  FourierFilter filter;
979  filter.FilterShape = REALGAUSSIAN;
980  filter.FilterBand = LOWPASS;
981  filter.w1 = 4;
982  filter.applyMaskSpace(im);
983 
984  // Image<double> imG(im);
985  // imG.write("I_Gauss.mrc");
986 
987  // std::cout << " - Calling varianceFilter() " << std::endl;
988  varianceFilter(I, varKernelSize, true);
989  im = I;
990 
991  // std::cout << " - Starting 2nd fft filtering " << std::endl;
992  filter.w1 = varKernelSize/8;
993  filter.applyMaskSpace(I);
994 
995  // Image<double> imGV(im);
996  // imGV.write("I_Gauss_Var.mrc");
997 
998  im -= im.computeMin();
999  im /= im.computeMax();
1000 
1001  // Image<double> imGVN(im);
1002  // imGVN.write("I_Gauss_Var_Norm.mrc");
1003 
1004  // std::cout << " - Starting histogram analysis " << std::endl;
1005  Histogram1D hist;
1006  hist.clear();
1007  compute_hist(im, hist, 256);
1008 
1009  // std::cout << " - Computing Gini coeff " << std::endl;
1010  MultidimArray<double> sortedList=hist;
1011  hist.sort(sortedList);
1012  double height=0;
1013  double area=0;
1014  for (int i=0; i<XSIZE(sortedList); i++)
1015  {
1016  height += DIRECT_MULTIDIM_ELEM(sortedList,i);
1017  area += height - DIRECT_MULTIDIM_ELEM(sortedList,i)/2.0;
1018  }
1019 
1020  double fair_area = height*XSIZE(hist)/2.0;
1021 
1022  double giniValue = (fair_area-area)/fair_area;
1023 
1024  return giniValue;
1025 }
1026 
1027 /* Otsu Segmentation ------------------------------------------------------- */
1028 double OtsuSegmentation(MultidimArray<double> &V)
1029 {
1030  // V.checkDimension(3);
1031 
1032  // Compute the probability density function
1033  Histogram1D hist;
1034  hist.clear();
1035  compute_hist(V, hist, 200);
1036  hist /= hist.sum();
1037 
1038  // Compute the cumulative 0th and 1st order moments
1039  double x;
1040  MultidimArray<double> mom0;
1041  MultidimArray<double> mom1;
1042  mom0.initZeros(XSIZE(hist));
1043  mom1.initZeros(XSIZE(hist));
1044  mom0(0) = hist(0);
1045  hist.index2val(0, x);
1046  mom1(0) = hist(0) * x;
1047  for (size_t i = 1; i < XSIZE(mom0); i++)
1048  {
1049  mom0(i) = mom0(i - 1) + hist(i);
1050  hist.index2val(i, x);
1051  mom1(i) = mom1(i - 1) + hist(i) * x;
1052  }
1053 
1054  // Maximize sigma2B
1055  double bestSigma2B = -1;
1056  int ibestSigma2B = -1;
1057  for (size_t i = 0; i < XSIZE(hist) - 1; i++)
1058  {
1059  double w1 = mom0(i);
1060  double w2 = 1 - mom0(i);
1061  double mu1 = mom1(i);
1062  double mu2 = mom1(XSIZE(mom1) - 1) - mom1(i);
1063  double sigma2B = w1 * w2 * (mu1 - mu2) * (mu1 - mu2);
1064  if (sigma2B > bestSigma2B)
1065  {
1066  bestSigma2B = sigma2B;
1067  ibestSigma2B = i;
1068  }
1069  }
1070 
1071  hist.index2val(ibestSigma2B, x);
1072  V.binarize(x);
1073 
1074  return x;
1075 }
1076 
1077 /* Entropy Segmentation ---------------------------------------------------- */
1078 double EntropySegmentation(MultidimArray<double> &V)
1079 {
1080  // V.checkDimension(3);
1081 
1082  // Compute the probability density function
1083  Histogram1D hist;
1084  hist.clear();
1085  compute_hist(V, hist, 200);
1086  hist /= hist.sum();
1087 
1088  // Compute the cumulative 0th and 1st order moments
1089  double x;
1090  MultidimArray<double> mom0;
1091  mom0.initZeros(XSIZE(hist));
1092  mom0(0) = hist(0);
1093  for (size_t i = 1; i < XSIZE(mom0); i++)
1094  mom0(i) = mom0(i - 1) + hist(i);
1095 
1096  // Entropy for black and white parts of the histogram
1097  const double epsilon = 1e-15;
1098  MultidimArray<double> h1;
1099  MultidimArray<double> h2;
1100  h1.initZeros(XSIZE(hist));
1101  h2.initZeros(XSIZE(hist));
1102  for (size_t i = 0; i < XSIZE(hist); i++)
1103  {
1104  // Entropy h1
1105  double w1 = mom0(i);
1106  if (w1 > epsilon)
1107  for (size_t ii = 0; ii <= i; ii++)
1108  if (hist(ii) > epsilon)
1109  {
1110  double aux = hist(ii) / w1;
1111  h1(i) -= aux * log10(aux);
1112  }
1113 
1114  // Entropy h2
1115  double w2 = 1 - mom0(i);
1116  if (w2 > epsilon)
1117  for (size_t ii = i + 1; ii < XSIZE(hist); ii++)
1118  if (hist(ii) > epsilon)
1119  {
1120  double aux = hist(ii) / w2;
1121  h2(i) -= aux * log10(aux);
1122  }
1123  }
1124 
1125  // Find histogram index with maximum entropy
1126  double Hmax = h1(0) + h2(0);
1127  size_t iHmax = 0;
1128  for (size_t i = 1; i < XSIZE(hist) - 1; i++)
1129  {
1130  double H = h1(i) + h2(i);
1131  if (H > Hmax)
1132  {
1133  Hmax = H;
1134  iHmax = i;
1135  }
1136  }
1137 
1138  hist.index2val(iHmax, x);
1139  V.binarize(x);
1140 
1141  return x;
1142 }
1143 
1144 /* Otsu+Entropy Segmentation ----------------------------------------------- */
1145 double EntropyOtsuSegmentation(MultidimArray<double> &V, double percentil,
1146  bool binarizeVolume)
1147 {
1148  //V.checkDimension(3);
1149 
1150  // Compute the probability density function
1151  Histogram1D hist;
1152  hist.clear();
1153  compute_hist(V, hist, 300);
1154  hist /= hist.sum();
1155 
1156  // Compute the cumulative 0th and 1st order moments
1157  double x;
1158  MultidimArray<double> mom0;
1159  MultidimArray<double> mom1;
1160  mom0.initZeros(XSIZE(hist));
1161  mom1.initZeros(XSIZE(hist));
1162  mom0(0) = hist(0);
1163  hist.index2val(0, x);
1164  mom1(0) = hist(0) * x;
1165  for (size_t i = 1; i < XSIZE(mom0); i++)
1166  {
1167  mom0(i) = mom0(i - 1) + hist(i);
1168  hist.index2val(i, x);
1169  mom1(i) = mom1(i - 1) + hist(i) * x;
1170  }
1171 
1172  // Entropy for black and white parts of the histogram
1173  const double epsilon = 1e-15;
1174  MultidimArray<double> h1;
1175  MultidimArray<double> h2;
1176  h1.initZeros(XSIZE(hist));
1177  h2.initZeros(XSIZE(hist));
1178  for (size_t i = 0; i < XSIZE(hist); i++)
1179  {
1180  // Entropy h1
1181  double w1 = mom0(i);
1182  if (w1 > epsilon)
1183  for (size_t ii = 0; ii <= i; ii++)
1184  if (hist(ii) > epsilon)
1185  {
1186  double aux = hist(ii) / w1;
1187  h1(i) -= aux * log10(aux);
1188  }
1189 
1190  // Entropy h2
1191  double w2 = 1 - mom0(i);
1192  if (w2 > epsilon)
1193  for (size_t ii = i + 1; ii < XSIZE(hist); ii++)
1194  if (hist(ii) > epsilon)
1195  {
1196  double aux = hist(ii) / w2;
1197  h2(i) -= aux * log10(aux);
1198  }
1199  }
1200 
1201  // Compute sigma2B and H
1202  MultidimArray<double> sigma2B;
1203  MultidimArray<double> H;
1204  MultidimArray<double> HSigma2B;
1205  sigma2B.initZeros(XSIZE(hist) - 1);
1206  H.initZeros(XSIZE(hist) - 1);
1207  HSigma2B.initZeros(XSIZE(hist) - 1);
1208  for (size_t i = 0; i < XSIZE(hist) - 1; i++)
1209  {
1210  double w1 = mom0(i);
1211  double w2 = 1 - mom0(i);
1212  double mu1 = mom1(i);
1213  double mu2 = mom1(XSIZE(mom1) - 1) - mom1(i);
1214  sigma2B(i) = w1 * w2 * (mu1 - mu2) * (mu1 - mu2);
1215  H(i) = h1(i) + h2(i);
1216  HSigma2B(i) = -log10(sigma2B(i)) / H(i);
1217  // The logic behind this expression is
1218  // Otsu: max sigma2B -> max log10(sigma2B) -> min -log10(sigma2B)
1219  // Entropy: max H -> max H -> min 1/H
1220  }
1221 
1222  // Sort HSigma2B and take a given percentage of it
1223  MultidimArray<double> HSigma2Bsorted;
1224  HSigma2B.sort(HSigma2Bsorted);
1225  int iTh = ROUND(XSIZE(HSigma2B)*percentil);
1226  double threshold = HSigma2Bsorted(iTh);
1227 
1228  // Find the first value within HSigma2B falling below this threshold
1229  iTh = 0;
1230  while (A1D_ELEM(HSigma2B,iTh) >= threshold)
1231  iTh++;
1232  iTh--;
1233  if (iTh <= 0)
1234  x = threshold;
1235  else
1236  hist.index2val(iTh, x);
1237 
1238  if (binarizeVolume)
1239  V.binarize(x);
1240  return x;
1241 }
1242 
1243 /* Fast correntropy -------------------------------------------------------- */
1244 double fastCorrentropy(const MultidimArray<double> &x,
1245  const MultidimArray<double> &y, double sigma,
1246  const GaussianInterpolator &G, const MultidimArray<int> &mask)
1247 {
1248  double retvalxy = 0;
1249  double isigma = 1.0 / sigma;
1250  int maskSum = 0;
1251  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
1252  {
1253  if (DIRECT_MULTIDIM_ELEM(mask,n))
1254  {
1255  retvalxy += G.getValue(
1256  isigma
1257  * (DIRECT_MULTIDIM_ELEM(x,n)
1258  - DIRECT_MULTIDIM_ELEM(y,n)));
1259  ++maskSum;
1260  }
1261  }
1262  if (0 == maskSum) {
1263  REPORT_ERROR(ERR_NUMERICAL, "maskSum is zero (0), which would lead to division by zero");
1264  }
1265  return (retvalxy / maskSum);
1266 }
1267 
1268 /* Imed distance ---------------------------------------------------------- */
1269 double imedDistance(const MultidimArray<double>& I1, const MultidimArray<double>& I2)
1270 {
1271  // [x,y]=meshgrid([-3:1:3],[-3:1:3])
1272  // format long // w=1/sqrt(2*pi)*exp(-0.5*(x.*x+y.*y))
1273  double *refW;
1274  double w[49]={
1275  0.000049233388666, 0.000599785460091, 0.002688051941039, 0.004431848411938, 0.002688051941039, 0.000599785460091, 0.000049233388666,
1276  0.000599785460091, 0.007306882745281, 0.032747176537767, 0.053990966513188, 0.032747176537767, 0.007306882745281, 0.000599785460091,
1277  0.002688051941039, 0.032747176537767, 0.146762663173740, 0.241970724519143, 0.146762663173740, 0.032747176537767, 0.002688051941039,
1278  0.004431848411938, 0.053990966513188, 0.241970724519143, 0.398942280401433, 0.241970724519143, 0.053990966513188, 0.004431848411938,
1279  0.002688051941039, 0.032747176537767, 0.146762663173740, 0.241970724519143, 0.146762663173740, 0.032747176537767, 0.002688051941039,
1280  0.000599785460091, 0.007306882745281, 0.032747176537767, 0.053990966513188, 0.032747176537767, 0.007306882745281, 0.000599785460091,
1281  0.000049233388666, 0.000599785460091, 0.002688051941039, 0.004431848411938, 0.002688051941039, 0.000599785460091, 0.000049233388666};
1282 
1283  MultidimArray<double> diffImage = I1 - I2;
1284  double imed = 0;
1285  int imiddle=YSIZE(I1)/2;
1286  int jmiddle=XSIZE(I1)/2;
1287  int R2max=imiddle*imiddle;
1288  int ysize=(int)YSIZE(I1);
1289  int xsize=(int)XSIZE(I1);
1290  for (int i=3; i<ysize-3; ++i)
1291  {
1292  int i2=(i-imiddle)*(i-imiddle);
1293  for (int j=3; j<xsize-3; ++j)
1294  {
1295  int j2=(j-jmiddle)*(j-jmiddle);
1296  if (i2+j2>R2max) // Measure only within the maximum circle
1297  continue;
1298  double diffi=DIRECT_A2D_ELEM(diffImage,i,j);
1299  int index=0;
1300  for (int ii=-3; ii<=3; ++ii)
1301  {
1302  refW = &w[index];
1303  index = index + 7;
1304  double *diffj=&DIRECT_A2D_ELEM(diffImage,i+ii,j-3);
1305  double imedAux=(*refW++)*(*diffj++);
1306  imedAux+=(*refW++)*(*diffj++);
1307  imedAux+=(*refW++)*(*diffj++);
1308  imedAux+=(*refW++)*(*diffj++);
1309  imedAux+=(*refW++)*(*diffj++);
1310  imedAux+=(*refW++)*(*diffj++);
1311  imedAux+=(*refW++)*(*diffj++);
1312 
1313  imed+=imedAux*diffi;
1314  }
1315  }
1316  }
1317  return sqrt(imed);
1318 }
1319 
1320 /* Imed distance ---------------------------------------------------------- */
1321 double imedNormalizedDistance(const MultidimArray<double>& I1, const MultidimArray<double>& I2)
1322 {
1323  // [x,y]=meshgrid([-3:1:3],[-3:1:3])
1324  // format long
1325  // w=1/sqrt(2*pi)*exp(-0.5*(x.*x+y.*y))
1326  const double w[7][7]={
1327  {0.000049233388666, 0.000599785460091, 0.002688051941039, 0.004431848411938, 0.002688051941039, 0.000599785460091, 0.000049233388666},
1328  {0.000599785460091, 0.007306882745281, 0.032747176537767, 0.053990966513188, 0.032747176537767, 0.007306882745281, 0.000599785460091},
1329  {0.002688051941039, 0.032747176537767, 0.146762663173740, 0.241970724519143, 0.146762663173740, 0.032747176537767, 0.002688051941039},
1330  {0.004431848411938, 0.053990966513188, 0.241970724519143, 0.398942280401433, 0.241970724519143, 0.053990966513188, 0.004431848411938},
1331  {0.002688051941039, 0.032747176537767, 0.146762663173740, 0.241970724519143, 0.146762663173740, 0.032747176537767, 0.002688051941039},
1332  {0.000599785460091, 0.007306882745281, 0.032747176537767, 0.053990966513188, 0.032747176537767, 0.007306882745281, 0.000599785460091},
1333  {0.000049233388666, 0.000599785460091, 0.002688051941039, 0.004431848411938, 0.002688051941039, 0.000599785460091, 0.000049233388666}
1334  };
1335  int imiddle=YSIZE(I1)/2;
1336  int jmiddle=XSIZE(I1)/2;
1337  int R2max=imiddle*imiddle;
1338  // Measure the average of both images
1339  double avg1=0;
1340  double avg2=0;
1341  double N=0;
1342  int ydim=YSIZE(I1);
1343  int xdim=XSIZE(I1);
1344  for (int i=3; i<ydim-3; ++i)
1345  {
1346  int i2=(i-imiddle)*(i-imiddle);
1347  for (int j=3; j<xdim-3; ++j)
1348  {
1349  int j2=(j-jmiddle)*(j-jmiddle);
1350  if (i2+j2>R2max) // Measure only within the maximum circle
1351  continue;
1352  avg1+=DIRECT_A2D_ELEM(I1,i,j);
1353  avg2+=DIRECT_A2D_ELEM(I2,i,j);
1354  ++N;
1355  }
1356  }
1357  avg1/=N;
1358  avg2/=N;
1359  double diffAvg=avg1-avg2;
1360 
1361  double imed = 0;
1362  double imed1=0;
1363  double imed2=0;
1364  for (int i=3; i<ydim-3; ++i)
1365  {
1366  int i2=(i-imiddle)*(i-imiddle);
1367  for (int j=3; j<xdim-3; ++j)
1368  {
1369  int j2=(j-jmiddle)*(j-jmiddle);
1370  if (i2+j2>R2max) // Measure only within the maximum circle
1371  continue;
1372  double subtractedI1ij=DIRECT_A2D_ELEM(I1,i,j)-avg1;
1373  double subtractedI2ij=DIRECT_A2D_ELEM(I2,i,j)-avg2;
1374  double diffi=DIRECT_A2D_ELEM(I1,i,j)-DIRECT_A2D_ELEM(I2,i,j)-diffAvg;
1375  for (int ii=-3; ii<=3; ++ii)
1376  {
1377  int i_ii=i+ii;
1378  for (int jj=-3; jj<=3; ++jj)
1379  {
1380  int j_jj=j+jj;
1381  //double diffj=DIRECT_A2D_ELEM(I1,i_ii,j_jj)-DIRECT_A2D_ELEM(I2,i_ii,j_jj)-diffAvg;
1382  double diffj=DIRECT_A2D_ELEM(I1,i,j)-DIRECT_A2D_ELEM(I2,i_ii,j_jj)-diffAvg;
1383  double wiijj=w[ii+3][jj+3];
1384  imed+=wiijj*diffi*diffj;
1385  imed1+=wiijj*subtractedI1ij*(DIRECT_A2D_ELEM(I1,i_ii,j_jj)-avg1);
1386  imed2+=wiijj*subtractedI2ij*(DIRECT_A2D_ELEM(I2,i_ii,j_jj)-avg2);
1387  }
1388  }
1389  }
1390  }
1391  return imed/sqrt((imed1+imed2)/2);
1392 }
1393 
1394 
1395 
1396 /* Correlation measured inside a mask defined by the standard deviation value of the pixels of first image ---- */
1397 double correlationMasked(const MultidimArray<double>& I1, const MultidimArray<double>& I2)
1398 {
1399  double mean1;
1400  double std1;
1401  double th1;
1402  I1.computeAvgStdev(mean1,std1);
1403  th1 = std1;
1404 
1405  // Estimate the mean and stddev within the mask
1406  double N1=0;
1407  double sumMI1=0;
1408  double sumMI2=0;
1409  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I1)
1410  {
1411  double p1=DIRECT_MULTIDIM_ELEM(I1,n);
1412  double p2=DIRECT_MULTIDIM_ELEM(I2,n);
1413  if(p1>=th1){
1414  sumMI1+=p1;
1415  sumMI2+=p2;
1416  N1+=1.0;
1417  }
1418  }
1419 
1420  double sumMI1I2=0.0;
1421  double sumMI1I1=0.0;
1422  double sumMI2I2=0.0;
1423  double iN1=1;
1424  double avgM1=0;
1425  double avgM2=0;
1426  double corrM1M2;
1427  if (N1>0){
1428  iN1=1.0/N1;
1429  avgM1=sumMI1*iN1;
1430  avgM2=sumMI2*iN1;
1431  }
1432  else
1433  return 0;
1434 
1435  double p1a;
1436  double p2a;
1437  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I1)
1438  {
1439  double p1=DIRECT_MULTIDIM_ELEM(I1,n);
1440  double p2=DIRECT_MULTIDIM_ELEM(I2,n);
1441  if (p1>th1){
1442  p1a=p1-avgM1;
1443  p2a=p2-avgM2;
1444  sumMI1I1 +=p1a*p1a;
1445  sumMI2I2 +=p2a*p2a;
1446  sumMI1I2+=p1a*p2a;
1447  }
1448  }
1449  corrM1M2=sumMI1I2/sqrt(sumMI1I1*sumMI2I2);
1450 
1451  return corrM1M2;
1452 }
1453 
1454 
1455 
1456 /* Weighted correlation measured with a weight defined by the difference between both images ---- */
1457 double correlationWeighted(MultidimArray<double>& I1, MultidimArray<double>& I2)
1458 {
1459 
1460  MultidimArray<double> Idiff;
1461  MultidimArray<double> I2Aligned;
1462  Matrix2D<double> M;
1463 
1464  I1.setXmippOrigin();
1465  I2.setXmippOrigin();
1466  I2Aligned=I2;
1467  alignImages(I1, I2Aligned, M, false);
1468  Idiff=I1;
1469  Idiff-=I2Aligned;
1470 
1471  double mean;
1472  double std;
1473  Idiff.computeAvgStdev(mean,std);
1474  Idiff.selfABS();
1475  double threshold=std;
1476 
1477  // Estimate the mean and stddev within the mask
1478  double N=0;
1479  double sumWI1=0;
1480  double sumWI2=0;
1481  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Idiff)
1482  {
1483  if (DIRECT_MULTIDIM_ELEM(Idiff,n)>threshold)
1484  {
1485  double p1=DIRECT_MULTIDIM_ELEM(I1,n);
1486  double p2Alg=DIRECT_MULTIDIM_ELEM(I2Aligned,n);
1487  sumWI1+=p1;
1488  sumWI2+=p2Alg;
1489  N+=1.0;
1490  }
1491  }
1492 
1493  double sumWI1I2=0.0;
1494  double sumWI1I1=0.0;
1495  double sumWI2I2=0.0;
1496  double iN;
1497  double avgW1;
1498  double avgW2;
1499  double corrW1W2;
1500  if(N>0){
1501  iN=1.0/N;
1502  avgW1=sumWI1*iN;
1503  avgW2=sumWI2*iN;
1504  }
1505  else
1506  return -1;
1507  double p1a;
1508  double p2a;
1509  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I1)
1510  {
1511  if (DIRECT_MULTIDIM_ELEM(Idiff,n)>threshold)
1512  {
1513  double p1=DIRECT_MULTIDIM_ELEM(I1,n);
1514  double p2Alg=DIRECT_MULTIDIM_ELEM(I2Aligned,n);
1515  p1a=p1-avgW1;
1516  p2a=p2Alg-avgW2;
1517 
1518  double w=DIRECT_MULTIDIM_ELEM(Idiff,n);
1519  double wp1a=w*p1a;
1520  double wp2a=w*p2a;
1521 
1522  sumWI1I1 +=wp1a*p1a;
1523  sumWI2I2 +=wp2a*p2a;
1524  sumWI1I2 +=wp1a*p2a;
1525  }
1526  }
1527  corrW1W2=sumWI1I2/sqrt(sumWI1I1*sumWI2I2);
1528 
1529  return corrW1W2;
1530 }
1531 
1532 
1533 
1534 
1535 double svdCorrelation(const MultidimArray<double>& I1, const MultidimArray<double>& I2, const MultidimArray< int >* mask)
1536 {
1537  // Copy input images into matrix2Ds
1538  Matrix2D<double> mI1;
1539  Matrix2D<double> mI2;
1540  I1.copy(mI1);
1541  I2.copy(mI2);
1542 
1543  // SVD Decomposition of I1
1544  Matrix2D<double> U1;
1545  Matrix2D<double> V1;
1546  Matrix1D<double> W1;
1547  svdcmp(mI1,U1,W1,V1);
1548 
1549  // Express I2 in the basis of I1
1550  Matrix2D<double> aux;
1551  Matrix2D<double> W2;
1552  matrixOperation_AB(mI2,V1,aux);
1553  matrixOperation_AtB(U1,aux,W2);
1554 
1555  // Make W2 to be diagonal
1556  FOR_ALL_ELEMENTS_IN_MATRIX2D(W2)
1557  if (i!=j)
1558  MAT_ELEM(W2,i,j)=0.0;
1559 
1560  // Recover I2p
1561  matrixOperation_ABt(W2,V1,aux);
1562  matrixOperation_AB(U1,aux,mI2);
1563 
1564  // Measure correlation
1565  MultidimArray<double> I2p;
1566  I2p=mI2;
1567 // Image<double> save;
1568 // save()=I1;
1569 // save.write("PPPI1.xmp");
1570 // save()=I2;
1571 // save.write("PPPI2.xmp");
1572 // save()=I2p;
1573 // save.write("PPPI2p.xmp");
1574 // std::cout << "Correlation= " << correlationIndex(I1,I2p,mask) << std::endl;
1575 // std::cout << "Press any key\n";
1576 // char c; std::cin >> c;
1577  // double ratio=I2p.computeStddev()/I2.computeStddev();
1578  return correlationIndex(I1,I2p,mask); //*ratio;
1579 }
1580 
1581 /* Covariance matrix ------------------------------------------------------ */
1582 void covarianceMatrix(const MultidimArray<double> &I, Matrix2D<double> &C)
1583 {
1584  Matrix2D<double> mI;
1585  I.copy(mI);
1586  subtractColumnMeans(mI);
1587  matrixOperation_AtA(mI,C);
1588  C*=1.0/(YSIZE(I)-1.0);
1589 }
1590 
1591 /* Best shift -------------------------------------------------------------- */
1592 template float bestShift(MultidimArray<float>&, float&, float&, const MultidimArray<int>*, int);
1593 template<typename T>
1594 T bestShift(MultidimArray<T> &Mcorr,
1595  T &shiftX, T &shiftY, const MultidimArray<int> *mask, int maxShift)
1596 {
1597  int imax = INT_MIN;
1598  int jmax;
1599  int i_actual;
1600  int j_actual;
1601  double xmax;
1602  double ymax;
1603  double avecorr;
1604  double stdcorr;
1605  double dummy;
1606  bool neighbourhood = true;
1607 
1608  /*
1609  Warning: for masks with a small number of non-zero pixels, this routine is NOT reliable...
1610  Anyway, maybe using a mask is not a good idea at al...
1611  */
1612 
1613  // Adjust statistics within shiftmask to average 0 and stddev 1
1614  if (mask != nullptr)
1615  {
1616  if ((*mask).sum() < 2)
1617  {
1618  shiftX = shiftY = 0.;
1619  return -1;
1620  }
1621  else
1622  {
1623  computeStats_within_binary_mask(*mask, Mcorr, dummy, dummy, avecorr,
1624  stdcorr);
1625  double istdcorr = 1.0 / stdcorr;
1626  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Mcorr)
1627  if (DIRECT_MULTIDIM_ELEM(*mask, n))
1628  DIRECT_MULTIDIM_ELEM(Mcorr, n) =
1629  (DIRECT_MULTIDIM_ELEM(Mcorr, n) - avecorr)
1630  * istdcorr;
1631  else
1632  DIRECT_MULTIDIM_ELEM(Mcorr, n) = 0.;
1633  }
1634  }
1635  else
1636  Mcorr.statisticsAdjust((T)0, (T)1);
1637 
1638  // Look for maximum shift
1639  if (maxShift==-1)
1640  Mcorr.maxIndex(imax, jmax);
1641  else
1642  {
1643  int maxShift2=maxShift*maxShift;
1644  double bestCorr=std::numeric_limits<T>::lowest();
1645  for (int i=-maxShift; i<=maxShift; i++)
1646  for (int j=-maxShift; j<=maxShift; j++)
1647  {
1648  if (i*i+j*j>maxShift2) // continue if the Euclidean distance is too far
1649  continue;
1650  else if (A2D_ELEM(Mcorr, i, j)>bestCorr)
1651  {
1652  imax=i;
1653  jmax=j;
1654  bestCorr=A2D_ELEM(Mcorr, imax, jmax);
1655  }
1656  }
1657  }
1658  double max = A2D_ELEM(Mcorr, imax, jmax);
1659 
1660  // Estimate n_max around the maximum
1661  int n_max = -1;
1662  while (neighbourhood)
1663  {
1664  n_max++;
1665  for (int i = -n_max; i <= n_max && neighbourhood; i++)
1666  {
1667  i_actual = i + imax;
1668  if (i_actual < STARTINGY(Mcorr) || i_actual > FINISHINGY(Mcorr))
1669  {
1670  neighbourhood = false;
1671  break;
1672  }
1673  for (int j = -n_max; j <= n_max && neighbourhood; j++)
1674  {
1675  j_actual = j + jmax;
1676  if (j_actual < STARTINGX(Mcorr) || j_actual > FINISHINGX(Mcorr))
1677  {
1678  neighbourhood = false;
1679  break;
1680  }
1681  else if (max / 1.414 > A2D_ELEM(Mcorr, i_actual, j_actual))
1682  {
1683  neighbourhood = false;
1684  break;
1685  }
1686  }
1687  }
1688  }
1689 
1690  // We have the neighbourhood => looking for the gravity centre
1691  xmax = ymax = 0.;
1692  double sumcorr = 0.;
1693  if (imax-n_max<STARTINGY(Mcorr))
1694  n_max=std::min(imax-STARTINGY(Mcorr),n_max);
1695  if (imax+n_max>FINISHINGY(Mcorr))
1696  n_max=std::min(FINISHINGY(Mcorr)-imax,n_max);
1697  if (jmax-n_max<STARTINGY(Mcorr))
1698  n_max=std::min(jmax-STARTINGX(Mcorr),n_max);
1699  if (jmax+n_max>FINISHINGY(Mcorr))
1700  n_max=std::min(FINISHINGX(Mcorr)-jmax,n_max);
1701  for (int i = -n_max; i <= n_max; i++)
1702  {
1703  i_actual = i + imax;
1704  for (int j = -n_max; j <= n_max; j++)
1705  {
1706  j_actual = j + jmax;
1707  double val = A2D_ELEM(Mcorr, i_actual, j_actual);
1708  ymax += i_actual * val;
1709  xmax += j_actual * val;
1710  sumcorr += val;
1711  }
1712  }
1713  if (sumcorr != 0)
1714  {
1715  shiftX = xmax / sumcorr;
1716  shiftY = ymax / sumcorr;
1717  }
1718  return max;
1719 }
1720 
1721 double bestShift(const MultidimArray<double> &I1, const MultidimArray< std::complex<double> > &FFTI1,
1722  const MultidimArray<double> &I2,
1723  double &shiftX, double &shiftY, CorrelationAux &aux,
1724  const MultidimArray<int> *mask, int maxShift)
1725 {
1726  I1.checkDimension(2);
1727  I2.checkDimension(2);
1728 
1729  MultidimArray<double> Mcorr;
1730 
1731  correlation_matrix(FFTI1, I2, Mcorr, aux);
1732  return bestShift(Mcorr, shiftX, shiftY, mask, maxShift);
1733 }
1734 
1735 double bestShift(const MultidimArray<double> &I1, const MultidimArray<double> &I2,
1736  double &shiftX, double &shiftY, CorrelationAux &aux,
1737  const MultidimArray<int> *mask, int maxShift)
1738 {
1739  aux.transformer1.FourierTransform((MultidimArray<double> &)I1, aux.FFT1, false);
1740  return bestShift(I1,aux.FFT1,I2,shiftX,shiftY,aux,mask,maxShift);
1741 }
1742 
1743 
1744 double bestShift(const MultidimArray< std::complex<double> > &FFTI1,
1745  const MultidimArray< std::complex<double> > &FFTI2,
1746  MultidimArray<double> &Mcorr,
1747  double &shiftX, double &shiftY, CorrelationAux &aux,
1748  const MultidimArray<int> *mask, int maxShift)
1749 {
1750  correlation_matrix(FFTI1, FFTI2, Mcorr, aux);
1751  return bestShift(Mcorr, shiftX, shiftY, mask, maxShift);
1752 }
1753 
1754 /* Best shift -------------------------------------------------------------- */
1755 void bestShift(const MultidimArray<double> &I1, const MultidimArray<double> &I2,
1756  double &shiftX, double &shiftY, double &shiftZ, CorrelationAux &aux,
1757  const MultidimArray<int> *mask)
1758 {
1759  I1.checkDimension(3);
1760  I2.checkDimension(3);
1761 
1762  int imax;
1763  int jmax;
1764  int kmax;
1765  int i_actual;
1766  int j_actual;
1767  int k_actual;
1768  double max;
1769  double xmax;
1770  double ymax;
1771  double zmax;
1772  double sumcorr;
1773  double avecorr;
1774  double stdcorr;
1775  double dummy;
1776  bool neighbourhood = true;
1777  MultidimArray<double> Mcorr;
1778 
1779  correlation_matrix(I1, I2, Mcorr, aux);
1780 
1781  /*
1782  Warning: for masks with a small number of non-zero pixels, this routine is NOT reliable...
1783  Anyway, maybe using a mask is not a good idea at al...
1784  */
1785 
1786  // Adjust statistics within shiftmask to average 0 and stddev 1
1787  if (mask != nullptr)
1788  {
1789  if ((*mask).sum() < 2)
1790  {
1791  shiftX = shiftY = shiftZ = 0.;
1792  return;
1793  }
1794  else
1795  {
1796  computeStats_within_binary_mask(*mask, Mcorr, dummy, dummy, avecorr,
1797  stdcorr);
1798  double istdcorr = 1.0 / stdcorr;
1799  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Mcorr)
1800  if (DIRECT_MULTIDIM_ELEM(*mask, n))
1801  DIRECT_MULTIDIM_ELEM(Mcorr, n) =
1802  (DIRECT_MULTIDIM_ELEM(Mcorr, n) - avecorr)
1803  * istdcorr;
1804  else
1805  DIRECT_MULTIDIM_ELEM(Mcorr, n) = 0.;
1806  }
1807  }
1808  else
1809  Mcorr.statisticsAdjust(0., 1.);
1810  Mcorr.maxIndex(kmax, imax, jmax);
1811  max = A3D_ELEM(Mcorr, kmax, imax, jmax);
1812 
1813  // Estimate n_max around the maximum
1814  int n_max = -1;
1815  while (neighbourhood)
1816  {
1817  n_max++;
1818  for (int k = -n_max; k <= n_max && neighbourhood; k++)
1819  {
1820  k_actual = k + kmax;
1821  if (k_actual < STARTINGZ(Mcorr) || k_actual > FINISHINGZ(Mcorr))
1822  {
1823  neighbourhood = false;
1824  break;
1825  }
1826  for (int i = -n_max; i <= n_max && neighbourhood; i++)
1827  {
1828  i_actual = i + imax;
1829  if (i_actual < STARTINGY(Mcorr) || i_actual > FINISHINGY(Mcorr))
1830  {
1831  neighbourhood = false;
1832  break;
1833  }
1834  for (int j = -n_max; j <= n_max && neighbourhood; j++)
1835  {
1836  j_actual = j + jmax;
1837  if (j_actual < STARTINGX(Mcorr)
1838  || j_actual > FINISHINGX(Mcorr))
1839  {
1840  neighbourhood = false;
1841  break;
1842  }
1843  else if (max
1844  / 1.414 > A3D_ELEM(Mcorr, k_actual, i_actual, j_actual))
1845  {
1846  neighbourhood = false;
1847  break;
1848  }
1849  }
1850  }
1851  }
1852  }
1853 
1854  // We have the neighbourhood => looking for the gravity centre
1855  zmax = xmax = ymax = sumcorr = 0.;
1856  if (kmax-n_max<STARTINGZ(Mcorr))
1857  n_max=std::min(kmax-STARTINGZ(Mcorr),n_max);
1858  if (kmax+n_max>FINISHINGZ(Mcorr))
1859  n_max=std::min(FINISHINGZ(Mcorr)-kmax,n_max);
1860  if (imax-n_max<STARTINGY(Mcorr))
1861  n_max=std::min(imax-STARTINGY(Mcorr),n_max);
1862  if (imax+n_max>FINISHINGY(Mcorr))
1863  n_max=std::min(FINISHINGY(Mcorr)-imax,n_max);
1864  if (jmax-n_max<STARTINGY(Mcorr))
1865  n_max=std::min(jmax-STARTINGX(Mcorr),n_max);
1866  if (jmax+n_max>FINISHINGY(Mcorr))
1867  n_max=std::min(FINISHINGX(Mcorr)-jmax,n_max);
1868  for (int k = -n_max; k <= n_max; k++)
1869  {
1870  k_actual = k + kmax;
1871  for (int i = -n_max; i <= n_max; i++)
1872  {
1873  i_actual = i + imax;
1874  for (int j = -n_max; j <= n_max; j++)
1875  {
1876  j_actual = j + jmax;
1877  double val = A3D_ELEM(Mcorr, k_actual, i_actual, j_actual);
1878  zmax += k_actual * val;
1879  ymax += i_actual * val;
1880  xmax += j_actual * val;
1881  sumcorr += val;
1882  }
1883  }
1884  }
1885  if (sumcorr != 0)
1886  {
1887  shiftX = xmax / sumcorr;
1888  shiftY = ymax / sumcorr;
1889  shiftZ = zmax / sumcorr;
1890  }
1891 }
1892 
1893 /* Best non-wrapping shift ------------------------------------------------- */
1894 //#define DEBUG
1895 void bestNonwrappingShift(const MultidimArray<double> &I1,
1896  const MultidimArray<double> &I2, double &shiftX, double &shiftY,
1897  CorrelationAux &aux)
1898 {
1899  aux.transformer1.FourierTransform((MultidimArray<double> &)I1, aux.FFT1, false);
1900  bestNonwrappingShift(I1,aux.FFT1,I2,shiftX,shiftY,aux);
1901 }
1902 
1903 void bestNonwrappingShift(const MultidimArray<double> &I1, const MultidimArray< std::complex<double> >&FFTI1,
1904  const MultidimArray<double> &I2, double &shiftX, double &shiftY,
1905  CorrelationAux &aux)
1906 {
1907  I1.checkDimension(2);
1908  I2.checkDimension(2);
1909 
1910  bestShift(I1, FFTI1, I2, shiftX, shiftY, aux);
1911  double bestCorr;
1912  double corr;
1913  MultidimArray<double> Iaux;
1914 
1915  translate(1, Iaux, I1, vectorR2(-shiftX, -shiftY), xmipp_transformation::DONT_WRAP);
1916  //I1.translate(vectorR2(-shiftX,-shiftY),Iaux, DONT_WRAP);
1917  bestCorr = corr = fastCorrelation(I2, Iaux);
1918  double finalX = shiftX;
1919  double finalY = shiftY;
1920 #ifdef DEBUG
1921 
1922  std::cout << "shiftX=" << shiftX << " shiftY=" << shiftY
1923  << " corr=" << corr << std::endl;
1924  ImageXmipp save;
1925  save()=I1;
1926  save.write("PPPI1.xmp");
1927  save()=I2;
1928  save.write("PPPI2.xmp");
1929  save()=Iaux;
1930  save.write("PPPpp.xmp");
1931 #endif
1932 
1933  Iaux.initZeros();
1934  double testX = (shiftX > 0) ? (shiftX - XSIZE(I1)) : (shiftX + XSIZE(I1));
1935  double testY = shiftY;
1936  translate(1, Iaux, I1, vectorR2(-testX, -testY), xmipp_transformation::DONT_WRAP);
1937  //I1.translate(vectorR2(-testX,-testY),Iaux,DONT_WRAP);
1938  corr = fastCorrelation(I2, Iaux);
1939  if (corr > bestCorr)
1940  finalX = testX;
1941 #ifdef DEBUG
1942 
1943  std::cout << "shiftX=" << testX << " shiftY=" << testY
1944  << " corr=" << corr << std::endl;
1945  save()=Iaux;
1946  save.write("PPPmp.xmp");
1947 #endif
1948 
1949  Iaux.initZeros();
1950  testX = shiftX;
1951  testY = (shiftY > 0) ? (shiftY - YSIZE(I1)) : (shiftY + YSIZE(I1));
1952  translate(1, Iaux, I1, vectorR2(-testX, -testY), xmipp_transformation::DONT_WRAP);
1953  //I1.translate(vectorR2(-testX,-testY),Iaux,DONT_WRAP);
1954  corr = fastCorrelation(I2, Iaux);
1955  if (corr > bestCorr)
1956  finalY = testY;
1957 #ifdef DEBUG
1958 
1959  std::cout << "shiftX=" << testX << " shiftY=" << testY
1960  << " corr=" << corr << std::endl;
1961  save()=Iaux;
1962  save.write("PPPpm.xmp");
1963 #endif
1964 
1965  Iaux.initZeros();
1966  testX = (shiftX > 0) ? (shiftX - XSIZE(I1)) : (shiftX + XSIZE(I1));
1967  testY = (shiftY > 0) ? (shiftY - YSIZE(I1)) : (shiftY + YSIZE(I1));
1968  translate(1, Iaux, I1, vectorR2(-testX, -testY), xmipp_transformation::DONT_WRAP);
1969  //I1.translate(vectorR2(-testX,-testY),Iaux,DONT_WRAP);
1970  corr = fastCorrelation(I2, Iaux);
1971  if (corr > bestCorr)
1972  {
1973  finalX = testX;
1974  finalY = testY;
1975  }
1976 #ifdef DEBUG
1977  std::cout << "shiftX=" << testX << " shiftY=" << testY
1978  << " corr=" << corr << std::endl;
1979  save()=Iaux;
1980  save.write("PPPmm.xmp");
1981 #endif
1982 
1983  shiftX = finalX;
1984  shiftY = finalY;
1985 }
1986 #undef DEBUG
1987 
1988 /* Best shift -------------------------------------------------------------- */
1989 double bestShiftRealSpace(const MultidimArray<double> &I1, MultidimArray<double> &I2,
1990  double &shiftX, double &shiftY,
1991  const MultidimArray<int> *mask, int maxShift, double shiftStep)
1992 {
1993  I1.checkDimension(2);
1994  I2.checkDimension(2);
1995 
1996  double bestCorr=-1e38;
1997 
1998  MultidimArray<double> alignedI2;
1999  MultidimArray<double> bestI2;
2000  int maxShift2=maxShift*maxShift;
2001  Matrix1D<double> shift(2);
2002  for (double y=-maxShift; y<=maxShift; y+=shiftStep)
2003  for (double x=-maxShift; x<=maxShift; x+=shiftStep)
2004  {
2005  if (y*y+x*x>maxShift2)
2006  continue;
2007  YY(shift)=y;
2008  XX(shift)=x;
2009  translate(xmipp_transformation::LINEAR,alignedI2,I2,shift,xmipp_transformation::DONT_WRAP,0.0);
2010 
2011  double corr=correlationIndex(I1,alignedI2,mask);
2012  if (corr>bestCorr)
2013  {
2014  bestI2=alignedI2;
2015  bestCorr=corr;
2016  shiftY=y;
2017  shiftX=x;
2018  }
2019  }
2020  I2=bestI2;
2021 
2022  return bestCorr;
2023 }
2024 
2025 /* Align two images -------------------------------------------------------- */
2026 AlignmentAux::AlignmentAux()
2027 {
2028  plans = nullptr;
2029 }
2030 AlignmentAux::~AlignmentAux()
2031 {
2032  delete plans;
2033 }
2034 
2035 void computeAlignmentTransforms(const MultidimArray<double>& I, AlignmentTransforms &ITransforms,
2036  AlignmentAux &aux, CorrelationAux &aux2)
2037 {
2038  aux2.transformer1.FourierTransform((MultidimArray<double> &)I, ITransforms.FFTI, false);
2039  polarFourierTransform<true>(I, ITransforms.polarFourierI, false, XSIZE(I) / 5, XSIZE(I) / 2, aux.plans, 1);
2040 }
2041 
2042 constexpr double SHIFT_THRESHOLD = 0.95; // Shift threshold in pixels.
2043 constexpr float ROTATE_THRESHOLD = 1.0 ; // Rotate threshold in degrees.
2044 constexpr double INITIAL_SHIFT_THRESHOLD = SHIFT_THRESHOLD + 1.0; // Shift threshold in pixels.
2045 constexpr float INITIAL_ROTATE_THRESHOLD = ROTATE_THRESHOLD + 1.0 ; // Rotate threshold in degrees.
2046 
2047 double alignImages(const MultidimArray<double>& Iref, const AlignmentTransforms& IrefTransforms, MultidimArray<double>& I,
2048  Matrix2D<double>&M, bool wrap, AlignmentAux &aux, CorrelationAux &aux2,
2049  RotationalCorrelationAux &aux3)
2050 {
2051  I.checkDimension(2);
2052 
2053  aux.ARS.initIdentity(3);
2054  aux.ASR.initIdentity(3);
2055  aux.IauxSR = I;
2056  aux.IauxRS = I;
2057 
2058  aux.rotationalCorr.resize(2 * IrefTransforms.polarFourierI.getSampleNoOuterRing() - 1);
2059  aux3.local_transformer.setReal(aux.rotationalCorr);
2060 
2061  double shiftXSR=INITIAL_SHIFT_THRESHOLD;
2062  double shiftYSR=INITIAL_SHIFT_THRESHOLD;
2063  double bestRotSR=INITIAL_ROTATE_THRESHOLD;
2064  double shiftXRS=INITIAL_SHIFT_THRESHOLD;
2065  double shiftYRS=INITIAL_SHIFT_THRESHOLD;
2066  double bestRotRS=INITIAL_ROTATE_THRESHOLD;
2067 
2068  // Align the image with the reference
2069  for (int i = 0; i < 3; i++)
2070  {
2071  // Shift then rotate
2072  if (((shiftXSR > SHIFT_THRESHOLD) || (shiftXSR < (-SHIFT_THRESHOLD))) ||
2073  ((shiftYSR > SHIFT_THRESHOLD) || (shiftYSR < (-SHIFT_THRESHOLD))))
2074  {
2075  bestNonwrappingShift(Iref, IrefTransforms.FFTI, aux.IauxSR, shiftXSR, shiftYSR, aux2);
2076  MAT_ELEM(aux.ASR,0,2) += shiftXSR;
2077  MAT_ELEM(aux.ASR,1,2) += shiftYSR;
2078  applyGeometry(xmipp_transformation::LINEAR, aux.IauxSR, I, aux.ASR, xmipp_transformation::IS_NOT_INV, wrap);
2079  }
2080 
2081  if (bestRotSR > ROTATE_THRESHOLD)
2082  {
2083  polarFourierTransform<true>(aux.IauxSR, aux.polarFourierI, true,
2084  XSIZE(Iref) / 5, XSIZE(Iref) / 2, aux.plans, 1);
2085 
2086  bestRotSR = best_rotation(IrefTransforms.polarFourierI, aux.polarFourierI, aux3);
2087  rotation2DMatrix(bestRotSR, aux.R);
2088  aux.ASR = aux.R * aux.ASR;
2089  applyGeometry(xmipp_transformation::LINEAR, aux.IauxSR, I, aux.ASR, xmipp_transformation::IS_NOT_INV, wrap);
2090  }
2091 
2092  // Rotate then shift
2093  if (bestRotRS > ROTATE_THRESHOLD)
2094  {
2095  polarFourierTransform<true>(aux.IauxRS, aux.polarFourierI, true,
2096  XSIZE(Iref) / 5, XSIZE(Iref) / 2, aux.plans, 1);
2097  bestRotRS = best_rotation(IrefTransforms.polarFourierI, aux.polarFourierI, aux3);
2098  rotation2DMatrix(bestRotRS, aux.R);
2099  aux.ARS = aux.R * aux.ARS;
2100  applyGeometry(xmipp_transformation::LINEAR, aux.IauxRS, I, aux.ARS, xmipp_transformation::IS_NOT_INV, wrap);
2101  }
2102 
2103  if (((shiftXRS > SHIFT_THRESHOLD) || (shiftXRS < (-SHIFT_THRESHOLD))) ||
2104  ((shiftYRS > SHIFT_THRESHOLD) || (shiftYRS < (-SHIFT_THRESHOLD))))
2105  {
2106  bestNonwrappingShift(Iref, IrefTransforms.FFTI, aux.IauxRS, shiftXRS, shiftYRS, aux2);
2107  MAT_ELEM(aux.ARS,0,2) += shiftXRS;
2108  MAT_ELEM(aux.ARS,1,2) += shiftYRS;
2109  applyGeometry(xmipp_transformation::LINEAR, aux.IauxRS, I, aux.ARS, xmipp_transformation::IS_NOT_INV, wrap);
2110  }
2111  }
2112 
2113  double corrRS = correlationIndex(aux.IauxRS, Iref);
2114  double corrSR = correlationIndex(aux.IauxSR, Iref);
2115  double corr;
2116  if (corrRS > corrSR)
2117  {
2118  I = aux.IauxRS;
2119  M = aux.ARS;
2120  corr = corrRS;
2121  }
2122  else
2123  {
2124  I = aux.IauxSR;
2125  M = aux.ASR;
2126  corr = corrSR;
2127  }
2128  return corr;
2129 }
2130 
2131 double alignImages(const MultidimArray<double>& Iref, MultidimArray<double>& I,
2132  Matrix2D<double>&M, bool wrap, AlignmentAux &aux, CorrelationAux &aux2,
2133  RotationalCorrelationAux &aux3)
2134 {
2135  Iref.checkDimension(2);
2136  AlignmentTransforms IrefTransforms;
2137  computeAlignmentTransforms(Iref,IrefTransforms, aux, aux2);
2138  return alignImages(Iref, IrefTransforms, I, M, wrap, aux, aux2, aux3);
2139 }
2140 
2141 double alignImages(const MultidimArray<double>& Iref, MultidimArray<double>& I,
2142  Matrix2D<double>&M, bool wrap)
2143 {
2144  AlignmentAux aux;
2145  CorrelationAux aux2;
2146  RotationalCorrelationAux aux3;
2147  return alignImages(Iref, I, M, wrap, aux, aux2, aux3);
2148 }
2149 
2150 double alignImagesConsideringMirrors(const MultidimArray<double>& Iref, const AlignmentTransforms& IrefTransforms,
2151  MultidimArray<double>& I, Matrix2D<double> &M, AlignmentAux& aux,
2152  CorrelationAux& aux2, RotationalCorrelationAux &aux3, bool wrap,
2153  const MultidimArray<int>* mask)
2154 {
2155  MultidimArray<double> Imirror;
2156  Matrix2D<double> Mmirror;
2157  Imirror = I;
2158  Imirror.selfReverseX();
2159  Imirror.setXmippOrigin();
2160 
2161  double corr=alignImages(Iref, IrefTransforms, I, M, wrap, aux, aux2, aux3);
2162  double corrMirror=alignImages(Iref, IrefTransforms, Imirror, Mmirror, wrap, aux, aux2, aux3);
2163  if (mask!=nullptr)
2164  {
2165  corr = correlationIndex(Iref, I, mask);
2166  corrMirror = correlationIndex(Iref, Imirror, mask);
2167  }
2168  double bestCorr = corr;
2169  if (corrMirror > bestCorr)
2170  {
2171  bestCorr = corrMirror;
2172  I = Imirror;
2173  M = Mmirror;
2174  MAT_ELEM(M,0,0) *= -1;
2175  MAT_ELEM(M,1,0) *= -1;
2176  }
2177  return bestCorr;
2178 }
2179 
2180 double alignImagesConsideringMirrors(const MultidimArray<double>& Iref, MultidimArray<double>& I,
2181  Matrix2D<double>&M, bool wrap)
2182 {
2183  AlignmentAux aux;
2184  CorrelationAux aux2;
2185  RotationalCorrelationAux aux3;
2186  return alignImagesConsideringMirrors(Iref, I, M, aux, aux2, aux3, wrap, nullptr);
2187 }
2188 
2189 double alignImagesConsideringMirrors(const MultidimArray<double>& Iref,
2190  MultidimArray<double>& I, Matrix2D<double> &M, AlignmentAux& aux,
2191  CorrelationAux& aux2, RotationalCorrelationAux &aux3, bool wrap,
2192  const MultidimArray<int>* mask)
2193 {
2194  AlignmentTransforms IrefTransforms;
2195  computeAlignmentTransforms(Iref,IrefTransforms, aux, aux2);
2196  return alignImagesConsideringMirrors(Iref, IrefTransforms, I, M, aux, aux2, aux3, wrap, mask);
2197 }
2198 
2199 void alignSetOfImages(MetaData &MD, MultidimArray<double>& Iavg, int Niter,
2200  bool considerMirror)
2201 {
2202  Image<double> I;
2203  MultidimArray<double> InewAvg;
2204  FileName fnImg;
2205  AlignmentAux aux;
2206  CorrelationAux aux2;
2207  RotationalCorrelationAux aux3;
2208  Matrix2D<double> M;
2209  size_t Nimgs;
2210  size_t Xdim;
2211  size_t Ydim;
2212  size_t Zdim;
2213  getImageSize(MD, Xdim, Ydim, Zdim, Nimgs);
2214  for (int n = 0; n < Niter; ++n)
2215  {
2216  bool lastIteration = (n == (Niter - 1));
2217  InewAvg.initZeros(Ydim, Xdim);
2218  InewAvg.setXmippOrigin();
2219  for (size_t objId : MD.ids())
2220  {
2221  MD.getValue(MDL_IMAGE, fnImg, objId);
2222  I.read(fnImg);
2223  I().setXmippOrigin();
2224  double corr;
2225  if (considerMirror)
2226  corr = alignImagesConsideringMirrors(Iavg, I(), M, aux, aux2,
2227  aux3, xmipp_transformation::WRAP);
2228  else
2229  corr = alignImages(Iavg, I(), M, xmipp_transformation::WRAP, aux, aux2, aux3);
2230  InewAvg += I();
2231  if (n == 0)
2232  Iavg = InewAvg;
2233  if (lastIteration)
2234  {
2235  double scale;
2236  double shiftx;
2237  double shifty;
2238  double psi;
2239  bool flip;
2240  transformationMatrix2Parameters2D(M, flip, scale, shiftx,
2241  shifty, psi);
2242  MD.setValue(MDL_FLIP, flip, objId);
2243  MD.setValue(MDL_SHIFT_X, shiftx, objId);
2244  MD.setValue(MDL_SHIFT_Y, shifty, objId);
2245  MD.setValue(MDL_ANGLE_PSI, psi, objId);
2246  MD.setValue(MDL_MAXCC, corr, objId);
2247  }
2248  }
2249  InewAvg /= Nimgs;
2250  Iavg = InewAvg;
2251  centerImage(Iavg, aux2, aux3, 4);
2252  }
2253 }
2254 
2255 double fastBestRotation(const MultidimArray<double>& IrefCyl,
2256  const MultidimArray<double>& Icyl, CorrelationAux &aux,
2257  VolumeAlignmentAux &aux2, double deltaAng)
2258 {
2259  correlation_matrix(IrefCyl, Icyl, aux2.corr, aux, false);
2260  STARTINGZ(aux2.corr) = STARTINGY(aux2.corr) = STARTINGX(aux2.corr) = 0;
2261  double bestCorr = A3D_ELEM(aux2.corr,0,STARTINGY(aux2.corr),0);
2262  double bestAngle = 0;
2263  for (int i = STARTINGY(aux2.corr) + 1; i <= FINISHINGY(aux2.corr); i++)
2264  {
2265  double corr = A3D_ELEM(aux2.corr,0,i,0);
2266  if (corr > bestCorr)
2267  {
2268  bestCorr = corr;
2269  bestAngle = i;
2270  }
2271  }
2272  bestAngle *= deltaAng * 180.0 / PI;
2273  return -bestAngle;
2274 }
2275 
2276 double fastBestRotationAroundZ(const MultidimArray<double>& IrefCyl,
2277  const MultidimArray<double>& I, CorrelationAux &aux,
2278  VolumeAlignmentAux &aux2)
2279 {
2280  double deltaAng = atan(2.0 / XSIZE(I));
2281  Matrix1D<double> v(3);
2282  XX(v) = 0;
2283  YY(v) = 0;
2284  ZZ(v) = 1;
2285  volume_convertCartesianToCylindrical(I, aux2.Icyl, 3, XSIZE(I) / 2, 1, 0,
2286  2 * PI, deltaAng, v);
2287  return fastBestRotation(IrefCyl,aux2.Icyl,aux,aux2,deltaAng);
2288 }
2289 
2290 double fastBestRotationAroundY(const MultidimArray<double>& IrefCyl,
2291  const MultidimArray<double>& I, CorrelationAux &aux,
2292  VolumeAlignmentAux &aux2)
2293 {
2294  double deltaAng = atan(2.0 / XSIZE(I));
2295  Matrix1D<double> v(3);
2296  XX(v) = 0;
2297  YY(v) = 1;
2298  ZZ(v) = 0;
2299  volume_convertCartesianToCylindrical(I, aux2.Icyl, 3, XSIZE(I) / 2, 1, 0,
2300  2 * PI, deltaAng, v);
2301  return -fastBestRotation(IrefCyl,aux2.Icyl,aux,aux2,deltaAng);
2302 }
2303 
2304 double fastBestRotationAroundX(const MultidimArray<double>& IrefCyl,
2305  const MultidimArray<double>& I, CorrelationAux &aux,
2306  VolumeAlignmentAux &aux2)
2307 {
2308  double deltaAng = atan(2.0 / XSIZE(I));
2309  Matrix1D<double> v(3);
2310  XX(v) = 1;
2311  YY(v) = 0;
2312  ZZ(v) = 0;
2313  volume_convertCartesianToCylindrical(I, aux2.Icyl, 3, XSIZE(I) / 2, 1, 0,
2314  2 * PI, deltaAng, v);
2315  return fastBestRotation(IrefCyl,aux2.Icyl,aux,aux2,deltaAng);
2316 }
2317 
2318 void fastBestRotation(const MultidimArray<double>& IrefCylZ,
2319  const MultidimArray<double>& IrefCylY,
2320  const MultidimArray<double>& IrefCylX,
2321  const MultidimArray<double>& I,
2322  const MultidimArray<double>& Icurrent,
2323  MultidimArray<double>& Ifinal,
2324  char axis,
2325  Matrix2D<double> &R,
2326  CorrelationAux &aux, VolumeAlignmentAux &aux2)
2327 {
2328  double bestAngle;
2329  if (axis=='Z')
2330  bestAngle = fastBestRotationAroundZ(IrefCylZ, Icurrent, aux, aux2);
2331  else if (axis=='Y')
2332  bestAngle = fastBestRotationAroundY(IrefCylY, Icurrent, aux, aux2);
2333  else
2334  bestAngle = fastBestRotationAroundX(IrefCylX, Icurrent, aux, aux2);
2335  Matrix2D<double> Raux;
2336  rotation3DMatrix(bestAngle, axis, Raux);
2337  R=Raux*R;
2338  applyGeometry(xmipp_transformation::LINEAR, Ifinal, I, R, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
2339 }
2340 
2341 void fastBestRotation(const MultidimArray<double>& IrefCylZ,
2342  const MultidimArray<double>& IrefCylY,
2343  const MultidimArray<double>& IrefCylX,
2344  MultidimArray<double>& I,
2345  const String &eulerAngles,
2346  Matrix2D<double> &R,
2347  CorrelationAux &aux, VolumeAlignmentAux &aux2)
2348 {
2349  R.initIdentity(4);
2350  fastBestRotation(IrefCylZ,IrefCylY,IrefCylX,I,I,aux2.I1,eulerAngles[0],R, aux,aux2);
2351  fastBestRotation(IrefCylZ,IrefCylY,IrefCylX,I,aux2.I1,aux2.I12,eulerAngles[1],R,aux,aux2);
2352  fastBestRotation(IrefCylZ,IrefCylY,IrefCylX,I,aux2.I12,aux2.I123,eulerAngles[2],R,aux,aux2);
2353  I=aux2.I123;
2354 }
2355 
2356 double bestRotationAroundZ(const MultidimArray<double>& Iref,
2357  const MultidimArray<double>& I, CorrelationAux &aux,
2358  VolumeAlignmentAux &aux2)
2359 {
2360  Iref.checkDimension(3);
2361  I.checkDimension(3);
2362  if (!I.sameShape(Iref))
2363  REPORT_ERROR(ERR_MULTIDIM_SIZE,
2364  "Both volumes should be of the same shape");
2365 
2366  double deltaAng = atan(2.0 / XSIZE(I));
2367  Matrix1D<double> v(3);
2368  XX(v) = 0;
2369  YY(v) = 0;
2370  ZZ(v) = 1;
2371  volume_convertCartesianToCylindrical(Iref, aux2.IrefCyl, 3, XSIZE(I) / 2, 1,
2372  0, 2 * PI, deltaAng, v);
2373  return fastBestRotationAroundZ(aux2.IrefCyl, I, aux, aux2);
2374 }
2375 
2376 /* Estimate 2D Gaussian ---------------------------------------------------- */
2377 /* See Brandle, Chen, Bischof, Lapp. Robust parametric and semi-parametric
2378  spot fitting for spot array images. 2000 */
2379 double unnormalizedGaussian2D(const Matrix1D<double> &r,
2380  const Matrix1D<double> &mu, const Matrix2D<double> &sigmainv)
2381 {
2382  double x = XX(r) - XX(mu);
2383  double y = YY(r) - YY(mu);
2384  return exp(
2385  -0.5
2386  * (sigmainv(0, 0) * x * x + 2 * sigmainv(0, 1) * x * y
2387  + sigmainv(1, 1) * y * y));
2388 }
2389 
2390 void estimateGaussian2D(const MultidimArray<double> &I, double &a, double &b,
2391  Matrix1D<double> &mu, Matrix2D<double> &sigma, bool estimateMu,
2392  int iterations)
2393 {
2394  I.checkDimension(2);
2395 
2396  MultidimArray<double> z(I);
2397 
2398  // Estimate b
2399  Histogram1D hist;
2400  compute_hist(z, hist, 100);
2401  b = hist.percentil(5);
2402 
2403  // Iteratively estimate all parameters
2404  for (int n = 0; n < iterations; n++)
2405  {
2406  // Reestimate z
2407  FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
2408  z(i, j) = XMIPP_MAX(I(i,j)-b,0);
2409 
2410  // Sum of z
2411  double T = z.sum();
2412 
2413  // Estimate center
2414  mu.initZeros(2);
2415  if (estimateMu)
2416  {
2417  FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
2418  {
2419  double val = z(i, j);
2420  XX(mu) += val * j;
2421  YY(mu) += val * i;
2422  }
2423  mu /= T;
2424  }
2425 
2426  // Estimate sigma
2427  sigma.initZeros(2, 2);
2428  FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
2429  {
2430  double val = z(i, j);
2431  double j_mu = j - XX(mu);
2432  double i_mu = i - YY(mu);
2433  sigma(0, 0) += val * j_mu * j_mu;
2434  sigma(0, 1) += val * i_mu * j_mu;
2435  sigma(1, 1) += val * i_mu * i_mu;
2436  }
2437  sigma(1, 0) = sigma(0, 1);
2438  sigma /= T;
2439 
2440  // Estimate amplitude
2441  Matrix2D<double> sigmainv = sigma.inv();
2442  Matrix1D<double> r(2);
2443  double G2 = 0;
2444  a = 0;
2445  FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
2446  {
2447  XX(r) = j;
2448  YY(r) = i;
2449  double G = unnormalizedGaussian2D(r, mu, sigmainv);
2450  a += z(i, j) * G;
2451  G2 += G * G;
2452  }
2453  a /= G2;
2454 
2455  // Reestimate b
2456  FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
2457  {
2458  XX(r) = j;
2459  YY(r) = i;
2460  double G = unnormalizedGaussian2D(r, mu, sigmainv);
2461  z(i, j) = I(i, j) - a * G;
2462  }
2463  compute_hist(z, hist, 100);
2464  b = hist.percentil(5);
2465  }
2466 }
2467 
2468 /* Fourier-Bessel decomposition. ------------------------------------------- */
2469 void fourierBesselDecomposition(const MultidimArray<double> &img_in,
2470  double r2, int k1, int k2)
2471 {
2472  img_in.checkDimension(2);
2473 
2474  for (int k = k1; k <= k2; k++)
2475  {
2476  int k_1 = k - 1;
2477 
2478  // Compute h and a,b coefficients
2479  double h = 0;
2480  double my5 = 0;
2481  if (k_1 != 0)
2482  {
2483  double my = 1 + PI * r2 / 2 / k_1;
2484  double my4 = my * k_1;
2485  double ntot = 4 * my4;
2486  h = 2 * PI / ntot;
2487  }
2488  else
2489  {
2490  double my = 1 + PI * r2 / 2;
2491  double my4 = my;
2492  double ntot = 4 * my4;
2493  h = 2 * PI / ntot;
2494  }
2495 
2496  MultidimArray<double> sine(CEIL(my5));
2497  FOR_ALL_ELEMENTS_IN_ARRAY1D(sine)
2498  sine(i) = sin((i + 1) * h);
2499 
2500  }
2501 }
2502 
2503 /* Harmonic decomposition. ------------------------------------------------- */
2504 //void harmonicDecomposition(const MultidimArray<double> &img_in,
2505 // MultidimArray<double> &v_out)
2506 //{}
2507 
2508 
2509 /* Shah energy ------------------------------------------------------------- */
2510 /* This function computes the current functional energy */
2511 double Shah_energy(const MultidimArray<double> &img,
2512  const MultidimArray<double> &surface_strength,
2513  const MultidimArray<double> &edge_strength, double K,
2514  const Matrix1D<double> &W)
2515 {
2516  img.checkDimension(2);
2517 
2518  int Ydim1 = YSIZE(img) - 1;
2519  int Xdim1 = XSIZE(img) - 1;
2520 
2521  double Kinv = 1.0 / K;
2522 
2523  /* Calculate surface energy */
2524  double E1 = 0.0;
2525  double E2 = 0.0;
2526  double E3 = 0.0;
2527  double E4 = 0.0;
2528  double w0 = VEC_ELEM(W,0);
2529  double w1 = VEC_ELEM(W,1);
2530  double w2 = VEC_ELEM(W,2);
2531  double w3 = VEC_ELEM(W,3);
2532  for (int i = 1; i < Ydim1; i++)
2533  {
2534  int ip1 = i + 1;
2535  int im1 = i - 1;
2536  for (int j = 1; j < Xdim1; j++)
2537  {
2538  int jp1 = j + 1;
2539  int jm1 = j - 1;
2540  /* Calculate data matching terms */
2541  double D = dAij(img, i, j);
2542  double F = dAij(surface_strength, i, j);
2543  double S = dAij(edge_strength, i, j);
2544  double diff = D - F;
2545  E1 += w0 * diff * diff;
2546  E3 += w2 * K * S * S;
2547 
2548  /* Calculate first derivative terms */
2549  double Fx = 0.5
2550  * (dAij(surface_strength, i, jp1)
2551  - dAij(surface_strength, i, jm1));
2552  double Fy = 0.5
2553  * (dAij(surface_strength, ip1, j)
2554  - dAij(surface_strength, im1, j));
2555  double Sx =
2556  0.5
2557  * (dAij(edge_strength, i, jp1)
2558  - dAij(edge_strength, i, jm1));
2559  double Sy =
2560  0.5
2561  * (dAij(edge_strength, ip1, j)
2562  - dAij(edge_strength, im1, j));
2563  double S_1 = 1 - S;
2564  E2 += w1 * S_1 * S_1 * (Fx * Fx + Fy * Fy);
2565  E4 += w3 * Kinv * (Sx * Sx + Sy * Sy);
2566  }
2567  }
2568 
2569  return E1 + E2 + E3 + E4; // Total energy
2570 }
2571 
2572 /* Update Surface Shah ----------------------------------------------------- */
2573 /* This routine performs one update to the edge estimate based
2574  on a finite differences solution to the following equation:
2575  0 = dE/df = dF/df - d(dF/dfx)/dx - d(dF/dfy)/dy
2576  + dd(dF/dfxx)/dxx + dd(dF/dfxy)/dxy + dd(dF/dfyy)/dyy */
2577 double Update_surface_Shah(MultidimArray<double> &img,
2578  MultidimArray<double> &surface_strength,
2579  MultidimArray<double> &edge_strength, const Matrix1D<double> &W)
2580 {
2581  img.checkDimension(2);
2582 
2583  double Diff = 0.0;
2584  double Norm = 0.0;
2585  size_t Ydim1 = YSIZE(img) - 1;
2586  size_t Xdim1 = XSIZE(img) - 1;
2587 
2588  /* Update surface estimate */
2589  double w0 = VEC_ELEM(W,0);
2590  double w1 = VEC_ELEM(W,1);
2591  for (size_t i = 1; i < Ydim1; i++)
2592  {
2593  int ip1 = i + 1;
2594  int im1 = i - 1;
2595  for (size_t j = 1; j < Xdim1; j++)
2596  {
2597  int jp1 = j + 1;
2598  int jm1 = j - 1;
2599  /* Calculate edge partial derivative terms */
2600  double S = dAij(edge_strength, i, j);
2601  double Sx =
2602  0.5
2603  * (dAij(edge_strength, i, jp1)
2604  - dAij(edge_strength, i, jm1));
2605  double Sy =
2606  0.5
2607  * (dAij(edge_strength, ip1, j)
2608  - dAij(edge_strength, im1, j));
2609 
2610  double nS = 1 - S;
2611  double nS2 = nS * nS;
2612 
2613  /* Calculate surface partial derivative terms (excluding central pixel) */
2614  double F;
2615  double D;
2616  F = D = dAij(img, i, j);
2617  double SS_i_jp1 = dAij(surface_strength, i, jp1);
2618  double SS_i_jm1 = dAij(surface_strength, i, jm1);
2619  double SS_ip1_j = dAij(surface_strength, ip1, j);
2620  double SS_im1_j = dAij(surface_strength, im1, j);
2621  double Fx = 0.5 * (SS_i_jp1 - SS_i_jm1);
2622  double Fy = 0.5 * (SS_ip1_j - SS_im1_j);
2623  double Fxx = SS_i_jp1 + SS_i_jm1;
2624  double Fyy = SS_ip1_j + SS_im1_j;
2625 
2626  /* Calculate surface partial derivative weights */
2627  double wFx = 4 * w1 * nS * Sx;
2628  double wFy = 4 * w1 * nS * Sy;
2629  double wFxx = -2 * w1 * nS2;
2630  double wFyy = -2 * w1 * nS2;
2631 
2632  /* Calculate new surface value */
2633  double Constant = -2 * w0 * D;
2634  double Central = -2 * w0 + 2 * wFxx + 2 * wFyy;
2635  double Neighbors = wFx * Fx + wFy * Fy + wFxx * Fxx + wFyy * Fyy;
2636 
2637  if (fabs(Central) > XMIPP_EQUAL_ACCURACY)
2638  F = (Constant + Neighbors) / Central;
2639  F = CLIP(F, 0, 1);
2640 
2641  // Compute the difference.
2642  double SS_i_j = dAij(surface_strength, i, j);
2643  Diff += fabs(SS_i_j - F);
2644  Norm += fabs(SS_i_j);
2645 
2646  // Update the new value.
2647  dAij(surface_strength, i, j) = F;
2648  }
2649  }
2650  return Diff / Norm; // Return the relative difference.
2651 }
2652 
2653 /* Update Edge Shah -------------------------------------------------------- */
2654 /* This routine performs one update to the edge estimate based
2655  on a finite differences solution to the following equation:
2656  0 = dE/ds = dF/ds - d(dF/dsx)/dx - d(dF/dsy)/dy */
2657 double Update_edge_Shah(MultidimArray<double> &img,
2658  MultidimArray<double> &surface_strength,
2659  MultidimArray<double> &edge_strength, double K,
2660  const Matrix1D<double> &W)
2661 {
2662  img.checkDimension(2);
2663 
2664  double Diff = 0.0;
2665  double Norm = 0.0;
2666  size_t Ydim1 = YSIZE(img) - 1;
2667  size_t Xdim1 = XSIZE(img) - 1;
2668  double Kinv = 1.0 / K;
2669 
2670  /* Update edge estimate */
2671  double w1 = VEC_ELEM(W,1);
2672  double w2 = VEC_ELEM(W,2);
2673  double w3 = VEC_ELEM(W,3);
2674  for (size_t i = 1; i < Ydim1; i++)
2675  {
2676  int ip1 = i + 1;
2677  int im1 = i - 1;
2678  for (size_t j = 1; j < Xdim1; j++)
2679  {
2680  int jp1 = j + 1;
2681  int jm1 = j - 1;
2682 
2683  /* Calculate first and second derivative terms */
2684  double Fx = 0.5
2685  * (dAij(surface_strength, i, jp1)
2686  - dAij(surface_strength, i, jm1));
2687  double Fy = 0.5
2688  * (dAij(surface_strength, ip1, j)
2689  - dAij(surface_strength, im1, j));
2690  double Constant = w1 * (Fx * Fx + Fy * Fy);
2691 
2692  /* Calculate weights for central pixel and neighbors */
2693  double Central = w2 * K + w3 * Kinv * 4;
2694  double Neighbors = w3 * Kinv
2695  * (dAij(edge_strength, im1, j) + dAij(edge_strength, ip1, j)
2696  + dAij(edge_strength, i, jm1)
2697  + dAij(edge_strength, i, jp1));
2698 
2699  /* Calculate new S value */
2700  double Old_edge_strength = dAij(edge_strength, i, j);
2701  double S = (Constant + Neighbors) / (Constant + Central);
2702  if (S < 0)
2703  dAij(edge_strength, i, j) *= 0.5;
2704  else if (S > 1)
2705  dAij(edge_strength, i, j) = 0.5
2706  * (dAij(edge_strength, i, j) + 1);
2707  else
2708  dAij(edge_strength, i, j) = S;
2709 
2710  // Compute the difference.
2711  Diff += fabs(dAij(edge_strength, i, j) - Old_edge_strength);
2712  Norm += fabs(Old_edge_strength);
2713  }
2714  }
2715  return Diff / Norm; // Return the relative difference.
2716 }
2717 
2718 /* Smoothing Shah ---------------------------------------------------------- */
2719 constexpr double SHAH_CONVERGENCE_THRESHOLD = 0.0001;
2720 void smoothingShah(MultidimArray<double> &img,
2721  MultidimArray<double> &surface_strength,
2722  MultidimArray<double> &edge_strength, const Matrix1D<double> &W,
2723  int OuterLoops, int InnerLoops, int RefinementLoops,
2724  bool adjust_range)
2725 {
2726 
2727  img.checkDimension(2);
2728 
2729  typeCast(img, surface_strength);
2730  if (adjust_range)
2731  surface_strength.rangeAdjust(0, 1);
2732  edge_strength.resizeNoCopy(img);
2733 
2734  for (int k = 1; k <= RefinementLoops; k++)
2735  {
2736  // Initialize Edge Image.
2737  edge_strength.initZeros();
2738 
2739  double diffsurface = MAXFLOAT; // Reset surface difference
2740  for (int i = 0;
2741  ((i < OuterLoops) && OuterLoops)
2742  || ((diffsurface > SHAH_CONVERGENCE_THRESHOLD)
2743  && !OuterLoops); i++)
2744  {
2745 
2746  /* std::cout << "Iteration ..." << i+1;*/
2747  /* Iteratively update surface estimate */
2748  for (int j = 0; j < InnerLoops; j++)
2749  diffsurface = Update_surface_Shah(img, surface_strength,
2750  edge_strength, W);
2751 
2752  /* Iteratively update edge estimate */
2753  for (int j = 0; j < InnerLoops; j++)
2754  Update_edge_Shah(img, surface_strength, edge_strength, k, W);
2755 
2756  /* Calculate new functional energy */
2757  Shah_energy(img, surface_strength, edge_strength, k, W);
2758  }
2759  }
2760 }
2761 
2762 /* Tomographic diffusion --------------------------------------------------- */
2763 //#define DEBUG
2764 double tomographicDiffusion(MultidimArray<double>& V,
2765  const Matrix1D<double>& alpha, double lambda)
2766 {
2767  V.checkDimension(3);
2768 
2769  double alphax = XX(alpha);
2770  double alphay = YY(alpha);
2771  double alphaz = ZZ(alpha);
2772  double diffx;
2773  double diffy;
2774  double diffz;
2775 
2776  // Compute regularization error
2777  double regError = 0;
2778  for (size_t z = 1; z < ZSIZE(V) - 1; z++)
2779  for (size_t y = 1; y < YSIZE(V) - 1; y++)
2780  for (size_t x = 1; x < XSIZE(V) - 1; x++)
2781  {
2782  diffx = DIRECT_A3D_ELEM(V,z,y,x+1) - DIRECT_A3D_ELEM(V,z,y,x-1);
2783  diffy = DIRECT_A3D_ELEM(V,z,y+1,x) - DIRECT_A3D_ELEM(V,z,y-1,x);
2784  diffz = DIRECT_A3D_ELEM(V,z+1,y,x) - DIRECT_A3D_ELEM(V,z-1,y,x);
2785  regError += sqrt(
2786  alphax * diffx * diffx + alphay * diffy * diffy
2787  + alphaz * diffz * diffz);
2788  }
2789  regError *= 0.5;
2790 
2791  // Compute the gradient of the regularization error
2792  MultidimArray<double> gradient;
2793  gradient.initZeros(V);
2794  for (size_t z = 2; z < ZSIZE(V) - 2; z++)
2795  for (size_t y = 2; y < YSIZE(V) - 2; y++)
2796  for (size_t x = 2; x < XSIZE(V) - 2; x++)
2797  {
2798  // First term
2799  double V000 = DIRECT_A3D_ELEM(V,z,y,x);
2800  double V_200 = DIRECT_A3D_ELEM(V,z,y,x-2);
2801  double V_110 = DIRECT_A3D_ELEM(V,z,y+1,x-1);
2802  double V_1_10 = DIRECT_A3D_ELEM(V,z,y-1,x-1);
2803  double V_101 = DIRECT_A3D_ELEM(V,z+1,y,x-1);
2804  double V_10_1 = DIRECT_A3D_ELEM(V,z-1,y,x-1);
2805  diffx = V000 - V_200;
2806  diffy = V_110 - V_1_10;
2807  diffz = V_101 - V_10_1;
2808  double t1 = diffx
2809  / sqrt(
2810  alphax * diffx * diffx + alphay * diffy * diffy
2811  + alphaz * diffz * diffz);
2812 
2813  // Second term
2814  double V200 = DIRECT_A3D_ELEM(V,z,y,x+2);
2815  double V110 = DIRECT_A3D_ELEM(V,z,y+1,x+1);
2816  double V1_10 = DIRECT_A3D_ELEM(V,z,y-1,x+1);
2817  double V101 = DIRECT_A3D_ELEM(V,z+1,y,x+1);
2818  double V10_1 = DIRECT_A3D_ELEM(V,z-1,y,x+1);
2819  diffx = V200 - V000;
2820  diffy = V110 - V1_10;
2821  diffz = V101 - V10_1;
2822  double t2 = diffx
2823  / sqrt(
2824  alphax * diffx * diffx + alphay * diffy * diffy
2825  + alphaz * diffz * diffz);
2826 
2827  // Third term
2828  double V0_20 = DIRECT_A3D_ELEM(V,z,y-2,x);
2829  double V0_11 = DIRECT_A3D_ELEM(V,z+1,y-1,x);
2830  double V0_1_1 = DIRECT_A3D_ELEM(V,z-1,y-1,x);
2831  diffx = V1_10 - V_1_10;
2832  diffy = V000 - V0_20;
2833  diffz = V0_11 - V0_1_1;
2834  double t3 = diffy
2835  / sqrt(
2836  alphax * diffx * diffx + alphay * diffy * diffy
2837  + alphaz * diffz * diffz);
2838 
2839  // Fourth term
2840  double V020 = DIRECT_A3D_ELEM(V,z,y+2,x);
2841  double V011 = DIRECT_A3D_ELEM(V,z+1,y+1,x);
2842  double V01_1 = DIRECT_A3D_ELEM(V,z-1,y+1,x);
2843  diffx = V110 - V_110;
2844  diffy = V020 - V000;
2845  diffz = V011 - V01_1;
2846  double t4 = diffy
2847  / sqrt(
2848  alphax * diffx * diffx + alphay * diffy * diffy
2849  + alphaz * diffz * diffz);
2850 
2851  // Fifth term
2852  double V00_2 = DIRECT_A3D_ELEM(V,z-2,y,x);
2853  diffx = V10_1 - V_10_1;
2854  diffy = V01_1 - V0_1_1;
2855  diffz = V000 - V00_2;
2856  double t5 = diffz
2857  / sqrt(
2858  alphax * diffx * diffx + alphay * diffy * diffy
2859  + alphaz * diffz * diffz);
2860 
2861  // Sixth term
2862  double V002 = DIRECT_A3D_ELEM(V,z+2,y,x);
2863  diffx = V101 - V_101;
2864  diffy = V011 - V0_11;
2865  diffz = V002 - V000;
2866  double t6 = diffz
2867  / sqrt(
2868  alphax * diffx * diffx + alphay * diffy * diffy
2869  + alphaz * diffz * diffz);
2870 
2871  // Compute gradient
2872  DIRECT_A3D_ELEM(gradient,z,y,x) = 0.5
2873  * (alphax * (t1 - t2) + alphay * (t3 - t4)
2874  + alphaz * (t5 - t6));
2875  }
2876 #ifdef DEBUG
2877  VolumeXmipp save;
2878  save()=V;
2879  save.write("PPPvolume.vol");
2880  save()=gradient;
2881  save.write("PPPgradient.vol");
2882  std::cout << "Press any key\n";
2883  char c;
2884  std::cin >> c;
2885 #endif
2886 
2887  // Update volume
2888  for (size_t z = 2; z < ZSIZE(V) - 2; z++)
2889  for (size_t y = 2; y < YSIZE(V) - 2; y++)
2890  for (size_t x = 2; x < XSIZE(V) - 2; x++)
2891  DIRECT_A3D_ELEM(V,z,y,x) -= lambda
2892  * DIRECT_A3D_ELEM(gradient,z,y,x);
2893 
2894  // Finish
2895  return regError;
2896 }
2897 #undef DEBUG
2898 
2899 /* Rotational invariant moments -------------------------------------------- */
2900 void rotationalInvariantMoments(const MultidimArray<double> &img,
2901  const MultidimArray<int> *mask, MultidimArray<double> &v_out)
2902 {
2903  img.checkDimension(2);
2904 
2905  // Prepare some variables
2906  double m_11 = 0;
2907  double m_02 = 0;
2908  double m_20 = 0;
2909  double m_12 = 0;
2910  double m_21 = 0;
2911  double m_03 = 0;
2912  double m_30 = 0;
2913  double normalize_x = 2.0 / XSIZE(img);
2914  double normalize_y = 2.0 / YSIZE(img);
2915 
2916  // Compute low-level moments
2917  FOR_ALL_ELEMENTS_IN_ARRAY2D(img)
2918  {
2919  if (mask != nullptr)
2920  if (!(*mask)(i, j))
2921  continue;
2922  double val = img(i, j);
2923  double dx = j * normalize_x;
2924  double dy = i * normalize_y;
2925  double dx2 = dx * dx;
2926  double dx3 = dx2 * dx;
2927  double dy2 = dy * dy;
2928  double dy3 = dy2 * dy;
2929  // m_00+=val;
2930  m_11 += val * dx * dy;
2931  m_02 += val * dy2;
2932  m_20 += val * dx2;
2933  m_12 += val * dx * dy2;
2934  m_21 += val * dx2 * dy;
2935  m_03 += val * dy3;
2936  m_30 += val * dx3;
2937  }
2938  //m_11/=m_00; m_02/=m_00; m_20/=m_00;
2939  //m_12/=m_00; m_21/=m_00; m_03/=m_00; m_30/=m_00;
2940 
2941  // Compute high-level rotational invariant moments
2942  v_out.resize(5);
2943  v_out(0) = m_20 + m_02;
2944  v_out(1) = (m_20 - m_02) * (m_20 - m_02) + 4 * m_11 * m_11;
2945  v_out(2) = (m_30 - 3 * m_12) * (m_30 - 3 * m_12)
2946  + (3 * m_21 - m_03) * (3 * m_21 - m_03);
2947  v_out(3) = (m_30 + m_12) * (m_30 + m_12) + (m_21 + m_03) * (m_21 + m_03);
2948  v_out(4) =
2949  (m_30 - 3 * m_12) * (m_30 + m_12)
2950  * ((m_30 + m_12) * (m_30 + m_12)
2951  - 3 * (m_21 + m_03) * (m_21 + m_03))
2952  + (3 * m_21 - m_03) * (m_21 + m_03)
2953  * (3 * (m_30 + m_12) * (m_30 + m_12)
2954  - (m_21 + m_03) * (m_21 + m_03));
2955  /*
2956  v_out( 5)=(m_20+m_02)*(
2957  (m_30+m_12)*(m_30+m_12)
2958  -3*(m_21+m_03)*(m_21+m_03))
2959  +4*m_11*(m_30+m_12)*(m_03+m_21);
2960  v_out( 6)=(3*m_21-m_03)*(m_12+m_30)*(
2961  (m_30+m_12)*(m_30+m_12)
2962  -3*(m_21+m_03)*(m_21+m_03))
2963  -(m_30-3*m_12)*(m_12+m_03)*(
2964  3*(m_30+m_12)*(m_30+m_12)
2965  -(m_21+m_03)*(m_21+m_03));
2966  */
2967 }
2968 
2969 /* Inertia moments --------------------------------------------------------- */
2970 void inertiaMoments(const MultidimArray<double> &img,
2971  const MultidimArray<int> *mask, Matrix1D<double> &v_out,
2972  Matrix2D<double> &u)
2973 {
2974  img.checkDimension(2);
2975 
2976  // Prepare some variables
2977  double m_11 = 0;
2978  double m_02 = 0;
2979  double m_20 = 0;
2980  double normalize_x = 2.0 / XSIZE(img);
2981  double normalize_y = 2.0 / YSIZE(img);
2982 
2983  // Compute low-level moments
2984  FOR_ALL_ELEMENTS_IN_ARRAY2D(img)
2985  {
2986  if (mask != nullptr)
2987  if (!(*mask)(i, j))
2988  continue;
2989  double val = img(i, j);
2990  double dx = j * normalize_x;
2991  double dy = i * normalize_y;
2992  double dx2 = dx * dx;
2993  double dy2 = dy * dy;
2994  m_11 += val * dx * dy;
2995  m_02 += val * dy2;
2996  m_20 += val * dx2;
2997  }
2998 
2999  // Compute the eigen values of the inertia matrix
3000  // [m_02 m_11
3001  // m_11 m_20]
3002  Matrix2D<double> A(2, 2);
3003  A(0, 0) = m_02;
3004  A(0, 1) = A(1, 0) = m_11;
3005  A(1, 1) = m_20;
3006  Matrix2D<double> v;
3007  svdcmp(A, u, v_out, v);
3008  v_out = v_out.sort();
3009 }
3010 
3011 /* Fill triangle ----------------------------------------------------------- */
3012 void fillTriangle(MultidimArray<double> &img, int *tx, int *ty, double color)
3013 {
3014  img.checkDimension(2);
3015 
3016  /*
3017  * Order in y values
3018  */
3019  int y1 = 0;
3020  while (!y1)
3021  {
3022  y1 = 1;
3023  for (int y2 = 0; y2 < 2; y2++)
3024  {
3025  if (ty[y2] > ty[y2 + 1]
3026  || (ty[y2] == ty[y2 + 1] && tx[y2] < tx[y2 + 1]))
3027  {
3028  int x1 = ty[y2];
3029  ty[y2] = ty[y2 + 1];
3030  ty[y2 + 1] = x1;
3031  x1 = tx[y2];
3032  tx[y2] = tx[y2 + 1];
3033  tx[y2 + 1] = x1;
3034  y1 = 0;
3035  }
3036  }
3037  }
3038 
3039  int dx1 = tx[1] - tx[0];
3040  int dx2 = tx[2] - tx[0];
3041  int dy1 = ty[1] - ty[0];
3042  int dy2 = ty[2] - ty[0];
3043 
3044  int sx1 = SGN0(dx1);
3045  int sx2 = SGN0(dx2);
3046  int sy1 = SGN0(dy1);
3047 
3048  int ix1 = ABS(dx1);
3049  int ix2 = ABS(dx2);
3050  int iy1 = ABS(dy1);
3051  int iy2 = ABS(dy2);
3052 
3053  int inc1 = XMIPP_MAX(ix1, iy1);
3054  int inc2 = XMIPP_MAX(ix2, iy2);
3055 
3056  int x1;
3057  int x2;
3058  int y2;
3059  int xl;
3060  int xr;
3061  x1 = x2 = y1 = y2 = 0;
3062  xl = xr = tx[0];
3063  int y = ty[0];
3064 
3065  while (y != ty[1])
3066  {
3067  int doit1 = 0;
3068  int doit2 = 0;
3069 
3070  while (!doit1)
3071  { /* Wait until y changes */
3072  x1 += ix1;
3073  y1 += iy1;
3074  if (x1 > inc1)
3075  {
3076  x1 -= inc1;
3077  xl += sx1;
3078  }
3079  if (y1 > inc1)
3080  {
3081  y1 -= inc1;
3082  y += sy1;
3083  doit1 = 1;
3084  }
3085  }
3086 
3087  while (!doit2)
3088  { /* Wait until y changes */
3089  x2 += ix2;
3090  y2 += iy2;
3091  if (x2 > inc2)
3092  {
3093  x2 -= inc2;
3094  xr += sx2;
3095  }
3096  if (y2 > inc2)
3097  {
3098  y2 -= inc2;
3099  doit2 = 1;
3100  }
3101  }
3102 
3103  for (int myx = xl; myx <= xr; myx++)
3104  img(y, myx) = color;
3105  }
3106 
3107  dx1 = tx[2] - tx[1];
3108  dy1 = ty[2] - ty[1];
3109 
3110  sx1 = SGN0(dx1);
3111  sy1 = SGN0(dy1);
3112 
3113  ix1 = ABS(dx1);
3114  iy1 = ABS(dy1);
3115 
3116  inc1 = XMIPP_MAX(ix1, iy1);
3117  xl = tx[1];
3118  x1 = 0;
3119 
3120  while (y != ty[2])
3121  {
3122  int doit1 = 0;
3123  int doit2 = 0;
3124 
3125  while (!doit1)
3126  { /* Wait until y changes */
3127  x1 += ix1;
3128  y1 += iy1;
3129  if (x1 > inc1)
3130  {
3131  x1 -= inc1;
3132  xl += sx1;
3133  }
3134  if (y1 > inc1)
3135  {
3136  y1 -= inc1;
3137  y += sy1;
3138  doit1 = 1;
3139  }
3140  }
3141 
3142  while (!doit2)
3143  { /* Wait until y changes */
3144  x2 += ix2;
3145  y2 += iy2;
3146  if (x2 > inc2)
3147  {
3148  x2 -= inc2;
3149  xr += sx2;
3150  }
3151  if (y2 > inc2)
3152  {
3153  y2 -= inc2;
3154  doit2 = 1;
3155  }
3156  }
3157 
3158  for (int myx = xl; myx <= xr; myx++)
3159  img(y, myx) = color;
3160  }
3161 }
3162 
3163 /* Local thresholding ------------------------------------------------------ */
3164 void localThresholding(MultidimArray<double> &img, double C, double dimLocal,
3165  MultidimArray<int> &result, MultidimArray<int> *mask)
3166 {
3167 
3168  // Convolve the input image with the kernel
3169  MultidimArray<double> convolved;
3170  convolved.initZeros(img);
3171  FOR_ALL_ELEMENTS_IN_ARRAY2D(convolved)
3172  {
3173  if (mask != nullptr)
3174  if (!(*mask)(i, j))
3175  continue;
3176  int ii0 = XMIPP_MAX(STARTINGY(convolved), FLOOR(i - dimLocal));
3177  int jj0 = XMIPP_MAX(STARTINGX(convolved), FLOOR(j - dimLocal));
3178  int iiF = XMIPP_MIN(FINISHINGY(convolved), CEIL(i + dimLocal));
3179  int jjF = XMIPP_MIN(FINISHINGX(convolved), CEIL(j + dimLocal));
3180  double N = 0;
3181  for (int ii = ii0; ii <= iiF; ii++)
3182  for (int jj = jj0; jj <= jjF; jj++)
3183  {
3184  if (mask == nullptr)
3185  {
3186  convolved(i, j) += img(ii, jj);
3187  ++N;
3188  }
3189  else if ((*mask)(i, j))
3190  {
3191  convolved(i, j) += img(ii, jj);
3192  ++N;
3193  }
3194  }
3195  if (N != 0)
3196  convolved(i, j) /= N;
3197  }
3198 
3199  // Subtract the original from the convolved image and threshold
3200  result.initZeros(img);
3201  FOR_ALL_ELEMENTS_IN_ARRAY2D(convolved)
3202  {
3203  if (mask != nullptr)
3204  if (!(*mask)(i, j))
3205  continue;
3206  if (img(i, j) - convolved(i, j) > C)
3207  result(i, j) = 1;
3208  }
3209 }
3210 
3211 /* Center translationally -------------------------------------------------- */
3212 void centerImageTranslationally(MultidimArray<double> &I, CorrelationAux &aux)
3213 {
3214  I.checkDimension(2);
3215 
3216  MultidimArray<double> Ix = I;
3217  Ix.selfReverseX();
3218  Ix.setXmippOrigin();
3219  MultidimArray<double> Iy = I;
3220  Iy.selfReverseY();
3221  Iy.setXmippOrigin();
3222  MultidimArray<double> Ixy = Ix;
3223  Ixy.selfReverseY();
3224  Ixy.setXmippOrigin();
3225 
3226  double meanShiftX = 0;
3227  double meanShiftY = 0;
3228  double shiftX;
3229  double shiftY;
3230  bestNonwrappingShift(I, Ix, meanShiftX, meanShiftY, aux);
3231  bestNonwrappingShift(I, Iy, shiftX, shiftY, aux);
3232  meanShiftX += shiftX;
3233  meanShiftY += shiftY;
3234  bestNonwrappingShift(I, Ixy, shiftX, shiftY, aux);
3235  meanShiftX += shiftX;
3236  meanShiftY += shiftY;
3237  meanShiftX /= 3;
3238  meanShiftY /= 3;
3239 
3240  Matrix1D<double> shift(2);
3241  VECTOR_R2(shift, -meanShiftX, -meanShiftY);
3242  MultidimArray<double> aux2 = I;
3243  translate(3, I, aux2, shift);
3244  //I.selfTranslateBSpline(3,shift);
3245 }
3246 
3247 /* Center rotationally ----------------------------------------------------- */
3248 void centerImageRotationally(MultidimArray<double> &I,
3249  RotationalCorrelationAux &aux)
3250 {
3251  I.checkDimension(2);
3252 
3253  MultidimArray<double> Ix = I;
3254  Ix.selfReverseX();
3255  Ix.setXmippOrigin();
3256 
3257  Polar_fftw_plans *plans = nullptr;
3258  Polar<std::complex<double> > polarFourierI;
3259  Polar<std::complex<double> > polarFourierIx;
3260  polarFourierTransform<true>(Ix, polarFourierIx, false, XSIZE(Ix) / 5,
3261  XSIZE(Ix) / 2, plans);
3262  polarFourierTransform<true>(I, polarFourierI, true, XSIZE(I) / 5,
3263  XSIZE(I) / 2, plans);
3264 
3265  MultidimArray<double> rotationalCorr;
3266  rotationalCorr.resize(2 * polarFourierI.getSampleNoOuterRing() - 1);
3267  aux.local_transformer.setReal(rotationalCorr);
3268  double bestRot = best_rotation(polarFourierIx, polarFourierI, aux);
3269 
3270  MultidimArray<double> auxI = I;
3271  rotate(3, I, auxI, -bestRot / 2, xmipp_transformation::WRAP);
3272  //I.selfRotateBSpline(3,-bestRot/2,WRAP);
3273 }
3274 
3275 /* Center both rotationally and translationally ---------------------------- */
3276 //#define DEBUG
3277 Matrix2D<double> centerImage(MultidimArray<double> &I, CorrelationAux &aux,
3278  RotationalCorrelationAux &aux2, int Niter, bool limitShift)
3279 {
3280  I.checkDimension(2);
3281 
3282  I.setXmippOrigin();
3283  double avg = I.computeAvg();
3284  I -= avg;
3285 
3286  MultidimArray<double> Ix;
3287  MultidimArray<double> Iy;
3288  MultidimArray<double> Ixy;
3289  MultidimArray<double> Iaux;
3290  Matrix2D<double> A;
3291  A.initIdentity(3);
3292  Iaux = I;
3293 
3294  MultidimArray<int> mask;
3295  mask.initZeros(I);
3296  BinaryCircularMask(mask, XSIZE(I) / 2);
3297 
3298  MultidimArray<double> lineY;
3299  MultidimArray<double> lineX;
3300  lineY.initZeros(YSIZE(I));
3301  STARTINGX(lineY) = STARTINGY(I);
3302  lineX.initZeros(XSIZE(I));
3303  STARTINGX(lineX) = STARTINGX(I);
3304 
3305  Polar_fftw_plans *plans = nullptr;
3306  Polar<std::complex<double> > polarFourierI;
3307  Polar<std::complex<double> > polarFourierIx;
3308  MultidimArray<double> rotationalCorr;
3309  Matrix2D<double> R;
3310  for (int i = 0; i < Niter; i++)
3311  {
3312  // Mask Iaux
3313  FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
3314  if (!A2D_ELEM(mask,i,j))
3315  A2D_ELEM(Iaux,i,j) = 0;
3316 
3317  // Center translationally
3318  Ix = Iaux;
3319  Ix.selfReverseX();
3320  Ix.setXmippOrigin();
3321  Iy = Iaux;
3322  Iy.selfReverseY();
3323  Iy.setXmippOrigin();
3324  Ixy = Ix;
3325  Ixy.selfReverseY();
3326  Ixy.setXmippOrigin();
3327 
3328  double meanShiftX = 0;
3329  double meanShiftY = 0;
3330  double shiftX;
3331  double shiftY;
3332  double Nx = 0;
3333  double Ny = 0;
3334  bestNonwrappingShift(Iaux, Ix, shiftX, shiftY, aux);
3335 #ifdef DEBUG
3336 
3337  ImageXmipp save;
3338  save()=Ix;
3339  save.write("PPPx.xmp");
3340  std::cout << "con Ix: " << shiftX << " " << shiftY << std::endl;
3341 #endif
3342 
3343  if (fabs(shiftX) < XSIZE(I) / 3 || !limitShift)
3344  {
3345  meanShiftX += shiftX;
3346  Nx++;
3347  }
3348  if (fabs(shiftY) < YSIZE(I) / 3 || !limitShift)
3349  {
3350  meanShiftY += shiftY;
3351  Ny++;
3352  }
3353  bestNonwrappingShift(Iaux, Iy, shiftX, shiftY, aux);
3354 #ifdef DEBUG
3355 
3356  save()=Iy;
3357  save.write("PPPy.xmp");
3358  std::cout << "con Iy: " << shiftX << " " << shiftY << std::endl;
3359 #endif
3360 
3361  if (fabs(shiftX) < XSIZE(I) / 3 || !limitShift)
3362  {
3363  meanShiftX += shiftX;
3364  Nx++;
3365  }
3366  if (fabs(shiftY) < YSIZE(I) / 3 || !limitShift)
3367  {
3368  meanShiftY += shiftY;
3369  Ny++;
3370  }
3371  bestNonwrappingShift(Iaux, Ixy, shiftX, shiftY, aux);
3372 #ifdef DEBUG
3373 
3374  save()=Ixy;
3375  save.write("PPPxy.xmp");
3376  std::cout << "con Ixy: " << shiftX << " " << shiftY << std::endl;
3377 #endif
3378 
3379  if (fabs(shiftX) < XSIZE(I) / 3 || !limitShift)
3380  {
3381  meanShiftX += shiftX;
3382  Nx++;
3383  }
3384  if (fabs(shiftY) < YSIZE(I) / 3 || !limitShift)
3385  {
3386  meanShiftY += shiftY;
3387  Ny++;
3388  }
3389  if (Nx > 0)
3390  meanShiftX /= Nx;
3391  if (Ny > 0)
3392  meanShiftY /= Ny;
3393 
3394  MAT_ELEM(A,0,2) += -meanShiftX / 2;
3395  MAT_ELEM(A,1,2) += -meanShiftY / 2;
3396  Iaux.initZeros();
3397  applyGeometry(xmipp_transformation::LINEAR, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
3398  FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
3399  if (!A2D_ELEM(mask,i,j))
3400  A2D_ELEM(Iaux,i,j) = 0;
3401 
3402 #ifdef DEBUG
3403 
3404  std::cout << "Iter " << i << std::endl;
3405  std::cout << "shift=" << -meanShiftX << "," << -meanShiftY << std::endl;
3406  save()=I;
3407  save.write("PPP.xmp");
3408  save()=Iaux;
3409  save.write("PPPshift.xmp");
3410 #endif
3411 
3412  // Center rotationally
3413  Ix = Iaux;
3414  Ix.selfReverseX();
3415  Ix.setXmippOrigin();
3416 
3417  polarFourierTransform<true>(Iaux, polarFourierI, true, XSIZE(I) / 5,
3418  XSIZE(I) / 2, plans);
3419  rotationalCorr.resizeNoCopy(
3420  2 * polarFourierI.getSampleNoOuterRing() - 1);
3421  aux2.local_transformer.setReal(rotationalCorr);
3422 
3423  polarFourierTransform<true>(Ix, polarFourierIx, false,
3424  XSIZE(Ix) / 5, XSIZE(Ix) / 2, plans);
3425  double bestRot = best_rotation(polarFourierIx, polarFourierI, aux2);
3426  bestRot = realWRAP(bestRot,0,180);
3427  if (bestRot > 90)
3428  bestRot = bestRot - 180;
3429 
3430  rotation2DMatrix(bestRot / 2, R);
3431  A = R * A;
3432  Iaux.initZeros();
3433  applyGeometry(xmipp_transformation::LINEAR, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
3434  FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
3435  if (!A2D_ELEM(mask,i,j))
3436  A2D_ELEM(Iaux,i,j) = 0;
3437 
3438 #ifdef DEBUG
3439 
3440  std::cout << "rot=" << -bestRot/2 << std::endl;
3441  save()=Iaux;
3442  save.write("PPProt.xmp");
3443 #endif
3444 
3445  // Remove horizontal/vertical ambiguity
3446  lineX.initZeros();
3447  lineY.initZeros();
3448  FOR_ALL_ELEMENTS_IN_ARRAY2D(Iaux)
3449  {
3450  double val = A2D_ELEM(Iaux,i,j);
3451  if (j == 0)
3452  lineY(i) = val;
3453  else if (lineY(i) < val)
3454  lineY(i) = val;
3455  if (i == 0)
3456  lineX(j) = val;
3457  else if (lineX(j) < val)
3458  lineX(j) = val;
3459  }
3460 
3461  double thX = lineX.computeMin()
3462  + 0.75 * (lineX.computeMax() - lineX.computeMin());
3463  double thY = lineY.computeMin()
3464  + 0.75 * (lineY.computeMax() - lineY.computeMin());
3465  int x0 = STARTINGX(lineX);
3466  while (lineX(x0) < thX)
3467  x0++;
3468  int y0 = STARTINGX(lineY);
3469  while (lineY(y0) < thY)
3470  y0++;
3471  int xF = FINISHINGX(lineX);
3472  while (lineX(xF) < thX)
3473  xF--;
3474  int yF = FINISHINGX(lineY);
3475  while (lineY(yF) < thY)
3476  yF--;
3477  if ((xF - x0) > (yF - y0))
3478  {
3479  rotation2DMatrix(-90., R);
3480  A = R * A;
3481  }
3482  applyGeometry(xmipp_transformation::LINEAR, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
3483 #ifdef DEBUG
3484 
3485  lineX.write("PPPlineX.txt");
3486  lineY.write("PPPlineY.txt");
3487  std::cout << "dev X=" << xF-x0 << std::endl;
3488  std::cout << "dev Y=" << yF-y0 << std::endl;
3489  save()=Iaux;
3490  save.write("PPPhorver.xmp");
3491  std::cout << "Press any key\n";
3492  char c;
3493  std::cin >> c;
3494 #endif
3495 
3496  }
3497  applyGeometry(xmipp_transformation::BSPLINE3, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
3498  I = Iaux;
3499  I += avg;
3500  delete plans;
3501  return A;
3502 }
3503 #undef DEBUG
3504 
3506 void forcePositive(MultidimArray<double> &V)
3507 {
3508  MultidimArray<char> mask(ZSIZE(V), YSIZE(V), XSIZE(V));
3509  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V)
3510  {
3511  double x = DIRECT_MULTIDIM_ELEM(V, n);
3512  DIRECT_MULTIDIM_ELEM(mask, n) = (x <= 0);
3513  }
3514  boundMedianFilter(V, mask);
3515 }
3516 
3517 void computeEdges(const MultidimArray<double>& vol,
3518  MultidimArray<double> &vol_edge)
3519 {
3520  MultidimArray<double> BSpline_coefs;
3521  BSpline_coefs.initZeros(vol);
3522  produceSplineCoefficients(3, BSpline_coefs, vol);
3523  vol_edge.initZeros(vol);
3524 
3525  FOR_ALL_ELEMENTS_IN_ARRAY3D(vol)
3526  {
3527  double V_dx;
3528  double V_dy;
3529  double V_dz;
3530  V_dx = interpolatedElementBSplineDiffX(BSpline_coefs, j, i, k, 3);
3531  V_dy = interpolatedElementBSplineDiffY(BSpline_coefs, j, i, k, 3);
3532  V_dz = interpolatedElementBSplineDiffZ(BSpline_coefs, j, i, k, 3);
3533  A3D_ELEM(vol_edge,k,i,j) = sqrt(
3534  (V_dx * V_dx) + (V_dy * V_dy) + (V_dz * V_dz));
3535 
3536  }
3537 }
3538 
3539 void forceDWTSparsity(MultidimArray<double> &V, double eps)
3540 {
3541  int size0=XSIZE(V);
3542  int sizeF=(int)NEXT_POWER_OF_2(size0);
3543  selfScaleToSize(xmipp_transformation::BSPLINE3,V,sizeF,sizeF,sizeF);
3544  MultidimArray<double> vol_wavelets;
3545  MultidimArray<double> vol_wavelets_abs;
3546  set_DWT_type(DAUB12);
3547  DWT(V,vol_wavelets);
3548  vol_wavelets_abs=vol_wavelets;
3549  vol_wavelets_abs.selfABS();
3550  double *begin=MULTIDIM_ARRAY(vol_wavelets_abs);
3551  double *end=MULTIDIM_ARRAY(vol_wavelets_abs)+MULTIDIM_SIZE(vol_wavelets_abs);
3552  std::sort(begin,end);
3553  double threshold1=DIRECT_MULTIDIM_ELEM(vol_wavelets_abs,
3554  (long int)((1.0-eps)*MULTIDIM_SIZE(vol_wavelets_abs)));
3555  vol_wavelets.threshold("abs_below", threshold1, 0.0);
3556  IDWT(vol_wavelets,V);
3557  selfScaleToSize(xmipp_transformation::BSPLINE3,V,size0, size0, size0);
3558 }
3559 
3561 
3563 void BadPixelFilter::defineParams(XmippProgram * program)
3564 {
3565  program->addParamsLine("== Bad pixels ==");
3566  program->addParamsLine(
3567  " [ --bad_pixels <type>] : Applied filters on bad pixels of the image.");
3568  program->addParamsLine(" where <type> ");
3569  program->addParamsLine(
3570  " negative : Applied at those negative values. Positive values are untouched.");
3571  program->addParamsLine(
3572  " mask <mask_file> : Applied at those pixels given by mask.");
3573  program->addParamsLine(
3574  " outliers <factor> : Applied at those pixels out of the range [mean - factor*std, mean + factor*std].");
3575  program->addParamsLine(" alias -b; ");
3576 }
3577 
3579 void BadPixelFilter::readParams(XmippProgram * program)
3580 {
3581  type = NEGATIVE;
3582  // Check operation to do
3583  String typeStr = program->getParam("--bad_pixels");
3584  if (typeStr == "negative")
3585  ; //nothing to do type already equal to NEGATIVE
3586  else if (typeStr == "mask")
3587  {
3588  mask = new Image<char>;
3589  mask->read(program->getParam("--bad_pixels", "mask"));
3590  type = MASK;
3591  }
3592  else if (typeStr == "outliers")
3593  {
3594  factor = program->getDoubleParam("--bad_pixels", "outliers");
3595  type = OUTLIER;
3596  }
3597 }
3598 
3600 void BadPixelFilter::apply(MultidimArray<double> &img)
3601 {
3602  switch (type)
3603  {
3604  case NEGATIVE:
3605  forcePositive(img);
3606  break;
3607  case MASK:
3608  boundMedianFilter(img, mask->data);
3609  break;
3610  case OUTLIER:
3611  pixelDesvFilter(img, factor);
3612  break;
3613 
3614  }
3615 }
3616 
3617 void LogFilter::defineParams(XmippProgram * program)
3618 {
3619  program->addParamsLine("== Log filter (for scanners) uses equation fa-fb*log(x+fc) ==");
3620  program->addParamsLine(" [--log]");
3621  program->addParamsLine(" [--fa <a>]");
3622  program->addParamsLine(" [--fb <b>]");
3623  program->addParamsLine(" [--fc <c>]");
3624 }
3625 
3627 void LogFilter::readParams(XmippProgram * program)
3628 {
3629  a = program->getDoubleParam("--fa");
3630  b = program->getDoubleParam("--fb");
3631  c = program->getDoubleParam("--fc");
3632 }
3633 
3635 void LogFilter::apply(MultidimArray<double> &img)
3636 {
3637  logFilter(img, a,b,c);
3638 }
3639 
3640 void DenoiseTVFilter::defineParams(XmippProgram * program)
3641 {
3642  program->addParamsLine("== Total Variation Denoising method for micrographs ==");
3643  program->addParamsLine(" [--denoiseTV]");
3644  program->addParamsLine(" [--maxIterTV <maxIter>]");
3645 }
3646 
3648 void DenoiseTVFilter::readParams(XmippProgram * program)
3649 {
3650  maxIter = program->getIntParam("--maxIterTV");
3651 }
3652 
3654 void DenoiseTVFilter::apply(MultidimArray<double> &img)
3655 {
3656  denoiseTVFilter(img, maxIter);
3657 }
3658 
3660 void BackgroundFilter::defineParams(XmippProgram * program)
3661 {
3662  program->addParamsLine("== Background removal ==");
3663  program->addParamsLine(
3664  " [ --background <type=plane> ] : Filters to remove the background.");
3665  program->addParamsLine(" where <type> ");
3666  program->addParamsLine(
3667  " plane : Remove the plane that best fits the pixels.");
3668  program->addParamsLine(
3669  " rollingball <radius> : The background is computed as a rolling ball operation.");
3670  program->addParamsLine(" alias -g; ");
3671 }
3672 
3674 void BackgroundFilter::readParams(XmippProgram * program)
3675 {
3676  type = PLANE;
3677  // Check operation to do
3678  String typeStr = program->getParam("--background");
3679  if (typeStr == "plane") //Nothing to do, plane by default
3680  ;
3681  else if (typeStr == "rollingball")
3682  {
3683  type = ROLLINGBALL;
3684  radius = program->getIntParam("--background", "rollingball");
3685  }
3686 }
3687 
3689 void BackgroundFilter::apply(MultidimArray<double> &img)
3690 {
3691  switch (type)
3692  {
3693  case PLANE:
3694  substractBackgroundPlane(img);
3695  break;
3696  case ROLLINGBALL:
3697  substractBackgroundRollingBall(img, radius);
3698  break;
3699 
3700  }
3701 }
3702 
3704 void MedianFilter::defineParams(XmippProgram * program)
3705 {
3706  program->addParamsLine("== Median ==");
3707  program->addParamsLine(" [ --median ] : Use median filter.");
3708  program->addParamsLine(" alias -m; ");
3709 }
3710 
3712 void MedianFilter::readParams(XmippProgram * program)
3713 { //Do nothing by now
3714 }
3715 
3717 void MedianFilter::apply(MultidimArray<double> &img)
3718 {
3719  static MultidimArray<double> tmp;
3720  tmp = img;
3721  medianFilter3x3(tmp, img);
3722 }
3723 
3725 void DiffusionFilter::defineParams(XmippProgram * program)
3726 {
3727  program->addParamsLine("== Anisotropic diffusion ==");
3728  program->addParamsLine(
3729  " [--diffusion] : Use anisotropic diffusion filter.");
3730  program->addParamsLine(
3731  " [--shah_iter+ <outer=10> <inner=1> <refinement=1>] : Diffusion outer, inner and refinement iterations");
3732  program->addParamsLine(" requires --diffusion;");
3733  program->addParamsLine(
3734  " [--shah_weight+ <w0=0> <w1=50> <w2=50> <w3=0.02>]:Diffusion weights");
3735  program->addParamsLine(
3736  " : w0 = data matching ");
3737  program->addParamsLine(
3738  " : w1 = 1st derivative smooth ");
3739  program->addParamsLine(
3740  " : w2 = edge strength ");
3741  program->addParamsLine(
3742  " : w3 = edge smoothness ");
3743  program->addParamsLine(" requires --diffusion;");
3744  program->addParamsLine(
3745  " [--shah_only_edge+] : Produce the edge image of the diffusion");
3746  program->addParamsLine(" requires --diffusion;");
3747 }
3748 
3750 void DiffusionFilter::readParams(XmippProgram * program)
3751 {
3752  Shah_weight.resizeNoCopy(4);
3753  Shah_outer = program->getIntParam("--shah_iter", 0);
3754  Shah_inner = program->getIntParam("--shah_iter", 1);
3755  Shah_refinement = program->getIntParam("--shah_iter", 2);
3756  Shah_weight(0) = program->getDoubleParam("--shah_weight", 0);
3757  Shah_weight(1) = program->getDoubleParam("--shah_weight", 1);
3758  Shah_weight(2) = program->getDoubleParam("--shah_weight", 2);
3759  Shah_weight(3) = program->getDoubleParam("--shah_weight", 3);
3760  Shah_edge = program->checkParam("--shah_only_edge");
3761 }
3762 
3763 void DiffusionFilter::show()
3764 {
3765  std::cout << " Shah difussion\n" << " Outer iterations " << Shah_outer
3766  << std::endl << " Inner iterations " << Shah_inner << std::endl
3767  << " Refinement interations " << Shah_refinement << std::endl
3768  << " Weight " << Shah_weight.transpose() << std::endl;
3769  if (Shah_edge)
3770  std::cout << " Generating edge image\n";
3771 
3772 }
3773 
3775 void DiffusionFilter::apply(MultidimArray<double> &img)
3776 {
3777  MultidimArray<double> surface_strength;
3778  MultidimArray<double> edge_strength;
3779  smoothingShah(img, surface_strength, edge_strength, Shah_weight, Shah_outer,
3780  Shah_inner, Shah_refinement);
3781  if (Shah_edge)
3782  img = edge_strength;
3783  else
3784  img = surface_strength;
3785 }
3786 
3788 void BasisFilter::defineParams(XmippProgram * program)
3789 {
3790  program->addParamsLine("== Basis filter ==");
3791  program->addParamsLine(
3792  " [--basis <file> <N=-1>] : Stack file with the basis, N is the number of elements to consider");
3793 }
3794 
3796 void BasisFilter::readParams(XmippProgram * program)
3797 {
3798  fnBasis = program->getParam("--basis");
3799  Nbasis = program->getIntParam("--basis", 1);
3800  basis.read(fnBasis);
3801  if (Nbasis > 0)
3802  basis().resize(Nbasis, 1, YSIZE(basis()), XSIZE(basis()));
3803 }
3804 
3805 void BasisFilter::show()
3806 {
3807  std::cout << " Basis filter\n" << " Basis file " << fnBasis << std::endl
3808  << " Number of basis " << Nbasis << std::endl;
3809 }
3810 
3812 void BasisFilter::apply(MultidimArray<double> &img)
3813 {
3814  const MultidimArray<double> &mBasis = basis();
3815  if (XSIZE(img) != XSIZE(mBasis) || YSIZE(img) != YSIZE(mBasis))
3816  REPORT_ERROR(ERR_MULTIDIM_SIZE,
3817  "Images and basis are of different size");
3818 
3819  MultidimArray<double> result;
3820  result.initZeros(img);
3821  for (size_t nn = 0; nn < NSIZE(mBasis); ++nn)
3822  {
3823  double cnn = 0;
3824  double *ptrBasis = &NZYX_ELEM(mBasis,nn,0,0,0);
3825  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img)
3826  cnn += DIRECT_MULTIDIM_ELEM(img,n) * (*ptrBasis++);
3827  ptrBasis = &NZYX_ELEM(mBasis,nn,0,0,0);
3828  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img)
3829  DIRECT_MULTIDIM_ELEM(result,n) += cnn * (*ptrBasis++);
3830  }
3831  img = result;
3832 }
3833 
3835 void RetinexFilter::defineParams(XmippProgram * program)
3836 {
3837  program->addParamsLine("== Retinex ==");
3838  program->addParamsLine(
3839  " [ --retinex <percentile=0.9> <mask_file=\"\"> <eps=1>] : Retinex filter, remove a given percentile of the");
3840  program->addParamsLine(" : Laplacian within the mask, epsilon controls the behaviour at low frequency");
3841 }
3842 
3844 void RetinexFilter::readParams(XmippProgram * program)
3845 {
3846  percentile = program->getDoubleParam("--retinex",0);
3847  String fnMask = program->getParam("--retinex",1);
3848  if (fnMask=="")
3849  mask = nullptr;
3850  else
3851  {
3852  mask = new Image<int>;
3853  mask->read(fnMask);
3854  }
3855  eps = program->getDoubleParam("--retinex",2);
3856 }
3857 void RetinexFilter::laplacian(const MultidimArray<double> &img,
3858  MultidimArray< std::complex<double> > &fimg, bool direct)
3859 {
3860  if (ZSIZE(img)>1)
3861  {
3862  Matrix1D<double> w(3);
3863  for (size_t k=0; k<ZSIZE(fimg); k++)
3864  {
3865  FFT_IDX2DIGFREQ(k,ZSIZE(img),ZZ(w));
3866  double filter_k=6+eps-2*cos(TWOPI*ZZ(w));
3867  for (size_t i=0; i<YSIZE(fimg); i++)
3868  {
3869  FFT_IDX2DIGFREQ(i,YSIZE(img),YY(w));
3870  double filter_ki=filter_k-2*cos(TWOPI*YY(w));
3871  for (size_t j=0; j<XSIZE(fimg); j++)
3872  {
3873  FFT_IDX2DIGFREQ(j,XSIZE(img),XX(w));
3874  double filter_kij=filter_ki-2*cos(TWOPI*XX(w));
3875  if (!direct && filter_kij>0)
3876  filter_kij=1.0/filter_kij;
3877  DIRECT_A3D_ELEM(fimg,k,i,j)*=filter_kij;
3878  }
3879  }
3880  }
3881  }
3882  else
3883  {
3884  Matrix1D<double> w(2);
3885  for (size_t i=0; i<YSIZE(fimg); i++)
3886  {
3887  FFT_IDX2DIGFREQ(i,YSIZE(img),YY(w));
3888  double filter_i=4+eps-2*cos(TWOPI*YY(w));
3889  for (size_t j=0; j<XSIZE(fimg); j++)
3890  {
3891  FFT_IDX2DIGFREQ(j,XSIZE(img),XX(w));
3892  double filter_ij=filter_i-2*cos(TWOPI*XX(w));
3893  if (!direct && filter_ij>0)
3894  filter_ij=1.0/filter_ij;
3895  DIRECT_A2D_ELEM(fimg,i,j)*=filter_ij;
3896  }
3897  }
3898  }
3899 }
3900 
3902 void RetinexFilter::apply(MultidimArray<double> &img)
3903 {
3904  FourierTransformer fft;
3905  MultidimArray< std::complex<double> > fimg;
3906 
3907  // Apply forward Laplacian
3908  fft.FourierTransform(img,fimg,false);
3909  laplacian(img,fimg,true);
3910  fft.inverseFourierTransform();
3911 
3912  // Compute the threshold outside the mask
3913  double *sortedValues=nullptr;
3914  size_t Nsorted=0;
3915  if (mask==nullptr)
3916  {
3917  Nsorted=MULTIDIM_SIZE(img);
3918  sortedValues=new double[Nsorted];
3919  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img)
3920  sortedValues[n]=fabs(DIRECT_MULTIDIM_ELEM(img,n));
3921  }
3922  else
3923  {
3924  const MultidimArray<int> &mmask=(*mask)();
3925  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img)
3926  if (DIRECT_MULTIDIM_ELEM(mmask,n)==0)
3927  Nsorted+=1;
3928 
3929  sortedValues=new double[Nsorted];
3930  size_t i=0;
3931  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img)
3932  if (DIRECT_MULTIDIM_ELEM(mmask,n)==0)
3933  sortedValues[i++]=fabs(DIRECT_MULTIDIM_ELEM(img,n));
3934  }
3935  std::sort(sortedValues,sortedValues+Nsorted);
3936  double threshold=sortedValues[(size_t)(percentile*Nsorted)];
3937  delete[] sortedValues;
3938 
3939  // Apply threshold
3940  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(img)
3941  if (fabs(DIRECT_MULTIDIM_ELEM(img,n))<threshold)
3942  DIRECT_MULTIDIM_ELEM(img,n)=0.0;
3943 
3944  fft.FourierTransform(img,fimg,false);
3945  laplacian(img,fimg,false);
3946  fft.inverseFourierTransform();
3947 }
3948 
3949 // Denoising using Spectral Projected Gradient (SPG) optimization
3959 double denoiseTVenergy(double mu,
3960  const MultidimArray<double>& X,
3961  const MultidimArray<double>& Y,
3962  double s,
3963  double q,
3964  double lambda,
3965  double sigmag,
3966  double g)
3967 {
3968  // parameter beta which role is to prevent division with zeros
3969  const double beta = 0.00001;
3970  const double beta2 = beta*beta;
3971 
3972  double m = 0;
3973  double tv = 0;
3974  const double K1 = ((3.0/8.0)*lambda*lambda + sigmag*sigmag - lambda*g) / (s*s);
3975  const double K2 = lambda * (q/(s*s));
3976  const double K3 = 2.0 / lambda;
3977  double Xij;
3978  double Yij;
3979  double dXx;
3980  double dXy;
3981  double d;
3982  double msqrt;
3983  int i1;
3984  int j1;
3985 
3986  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(X)
3987  {
3988  Xij = DIRECT_A2D_ELEM(X,i,j);
3989  Yij = DIRECT_A2D_ELEM(Y,i,j);
3990  i1 = i + 1;
3991  if (i1==YSIZE(X))
3992  i1 = 0;
3993  j1 = j + 1;
3994  if (j1==XSIZE(X))
3995  j1 = 0;
3996 
3997  dXx = DIRECT_A2D_ELEM(X,i,j1) - Xij;
3998  dXy = DIRECT_A2D_ELEM(X,i1,j) - Xij;
3999  tv += sqrt(dXx*dXx + dXy*dXy + beta2);
4000 
4001  d = K2*Xij + K1;
4002  msqrt = K3 * sqrt(std::max(0.0, d)) - Yij;
4003  m += msqrt*msqrt;
4004  }
4005  return 0.5*m + mu*tv;
4006 }
4007 
4017 void denoiseTVgradient(double mu,
4018  const MultidimArray<double>& X,
4019  const MultidimArray<double>& Y,
4020  double s,
4021  double q,
4022  double lambda,
4023  double sigmag,
4024  double g,
4025  MultidimArray<double>& gradientI
4026  )
4027 {
4028  // parameter beta which role is to prevent division with zeros
4029  const double beta = 0.00001;
4030  const double beta2 = beta*beta;
4031 
4032  MultidimArray<double> d = X;
4033  double K1 = ((3.0/8.0)*lambda*lambda + sigmag*sigmag - lambda*g) / (s*s);
4034  double K2 = lambda * (q/s*s);
4035  double K3 = (2.0 / (lambda*lambda)) * (q / (s*s)) * lambda;
4036  double Xij;
4037  double Yij;
4038  double dXx;
4039  double dXy;
4040  double dij;
4041  double d_left;
4042  double d_up;
4043  double X_left;
4044  double X_right;
4045  double X_up;
4046  double X_down;
4047  double dTV;
4048  double dE;
4049  int i1;
4050  int j1;
4051  int i_1;
4052  int j_1;
4053 
4054  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(X)
4055  {
4056  Xij = DIRECT_A2D_ELEM(X,i,j);
4057  i1 = i + 1;
4058  if (i1==YSIZE(X))
4059  i1 = 0;
4060  j1 = j + 1;
4061  if (j1==XSIZE(X))
4062  j1 = 0;
4063 
4064  dXx = DIRECT_A2D_ELEM(X,i,j1) - Xij;
4065  dXy = DIRECT_A2D_ELEM(X,i1,j) - Xij;
4066  DIRECT_A2D_ELEM(d,i,j) = 1.0 / sqrt(dXx*dXx + dXy*dXy + beta2);
4067  }
4068 
4069  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(d)
4070  {
4071  dij = DIRECT_A2D_ELEM(d,i,j);
4072  Xij = DIRECT_A2D_ELEM(X,i,j);
4073  Yij = DIRECT_A2D_ELEM(Y,i,j);
4074  i1 = i + 1;
4075  if (i1==YSIZE(d))
4076  i1 = 0;
4077  j1 = j + 1;
4078  if (j1==XSIZE(d))
4079  j1 = 0;
4080  i_1 = i - 1;
4081  if (i_1==-1)
4082  i_1 = YSIZE(d) - 1;
4083  j_1 = j - 1;
4084  if (j_1==-1)
4085  j_1 = XSIZE(d) - 1;
4086 
4087  d_left = DIRECT_A2D_ELEM(d,i,j_1);
4088  d_up = DIRECT_A2D_ELEM(d,i_1,j);
4089  X_left = DIRECT_A2D_ELEM(X,i,j_1);
4090  X_right = DIRECT_A2D_ELEM(X,i,j1);
4091  X_up = DIRECT_A2D_ELEM(X,i_1,j);
4092  X_down = DIRECT_A2D_ELEM(X,i1,j);
4093 
4094  dTV = Xij * (2.0*dij + d_left + d_up) - X_left*d_left - X_up*d_up - dij*(X_right + X_down);
4095 
4096  if (K2*Xij + K1 > 0)
4097  dE = K3 - (q / (s*s)) * Yij / sqrt(Xij * (q / (s*s)) * lambda + K1);
4098  else
4099  dE = 0;
4100 
4101  DIRECT_A2D_ELEM(gradientI,i,j) = dE + mu*dTV;
4102  }
4103 }
4104 
4108 void denoiseTVproj(const MultidimArray<double>& X,
4109  const MultidimArray<double>& Y,
4110  double theta,
4111  MultidimArray<double>& dold)
4112 {
4113  double div;
4114  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(X)
4115  {
4116  div = DIRECT_A2D_ELEM(X,i,j) - (DIRECT_A2D_ELEM(Y,i,j)*theta);
4117  if (div < 0)
4118  DIRECT_A2D_ELEM(dold,i,j) = 0 - DIRECT_A2D_ELEM(X,i,j);
4119  else if (div > 1)
4120  DIRECT_A2D_ELEM(dold,i,j) = 1 - DIRECT_A2D_ELEM(X,i,j);
4121  else
4122  DIRECT_A2D_ELEM(dold,i,j) = div - DIRECT_A2D_ELEM(X,i,j);
4123  }
4124 }
4125 
4129 void denoiseTVFilter(MultidimArray<double> &xnew, int maxIter)
4130 {
4131  // parameters for denoising
4132  double lambda = 1.0; // Poisson scaling factor
4133  double sigmag = 5.8; // Gaussian component N(g,sigma^2)
4134  double g = 0.0;
4135 
4136  // *** Anscombe transformation of input image ***
4137  double q = 255.0;
4138 
4139  double xnewmin = xnew.computeMin();
4140  double xnewscale = 255.0 / (xnew.computeMax() - xnewmin);
4141 
4142  // apply generalized Anscombe VST tranformation
4143  double K1 = ((3.0/8.0)*lambda*lambda + sigmag*sigmag - lambda*g);
4144  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(xnew)
4145  {
4146  DIRECT_MULTIDIM_ELEM(xnew, n) = (DIRECT_MULTIDIM_ELEM(xnew, n) - xnewmin) * xnewscale;
4147  DIRECT_MULTIDIM_ELEM(xnew, n) = 2.0 / lambda * sqrt(std::max(0.0, lambda*DIRECT_MULTIDIM_ELEM(xnew, n) + K1));
4148  }
4149 
4150  double mx = xnew.computeMax();
4151 
4152  // xold = initial degraded image
4153  MultidimArray<double> xold = xnew;
4154  MultidimArray<double> grold = xnew;
4155  MultidimArray<double> grnew = xnew;
4156  MultidimArray<double> dold = xnew;
4157  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
4158  {
4159  DIRECT_A2D_ELEM(xold,i,j) = DIRECT_A2D_ELEM(xnew,i,j) / mx;
4160  }
4161  MultidimArray<double> origInput = xold;
4162 
4163  // *** Spectral Projected Gradient (SPG) optimization ***
4164  // parameters of optimization
4165  double tol = 0;
4166  double theta;
4167  double thetamin = 0.001;
4168  double thetamax = 1000;
4169  double gamma = 0.0001;
4170  double sigma1 = 0.1;
4171  double sigma2 = 0.9;
4172  double mu = 0.03;
4173 
4174  // initial objective function value of xold
4175  double fold = denoiseTVenergy(mu, xold, origInput, mx, q, lambda, sigmag, g);
4176 
4177  // gradient of objective function
4178  denoiseTVgradient(mu, xold, origInput, mx, q, lambda, sigmag, g, grold);
4179 
4180  denoiseTVproj(xold, grold, 1.0, dold);
4181 
4182  double delta;
4183  double ksi;
4184  double fnew;
4185  double s_norm;
4186  double p;
4187  double s2;
4188  double xij;
4189  for (int kk=1; kk <= maxIter; kk++)
4190  {
4191  delta = 0;
4192 
4193  // make a step
4194  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
4195  {
4196  DIRECT_A2D_ELEM(xnew,i,j) = DIRECT_A2D_ELEM(xold,i,j) + DIRECT_A2D_ELEM(dold,i,j);
4197  delta += DIRECT_A2D_ELEM(grold,i,j) * DIRECT_A2D_ELEM(dold,i,j);
4198  }
4199 
4200  // objective function of xnew
4201  fnew = denoiseTVenergy(mu, xnew, origInput, mx, q, lambda, sigmag, g);
4202 
4203  ksi = 1.0;
4204 
4205  while (fnew > fold + gamma*ksi*delta)
4206  {
4207  double ksitsl = -0.5 * (ksi*ksi) * delta / (fnew - fold - ksi*delta);
4208  if ((ksitsl >= sigma1) && (ksitsl <= sigma2*ksi))
4209  ksi = ksitsl;
4210  else
4211  ksi = ksi / 2.0;
4212 
4213  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
4214  {
4215  DIRECT_A2D_ELEM(xnew,i,j) = DIRECT_A2D_ELEM(xold,i,j) + ksi*DIRECT_A2D_ELEM(dold,i,j);
4216  }
4217 
4218  // because we updated xnew, update also function value
4219  fnew = denoiseTVenergy(mu, xnew, origInput, mx, q, lambda, sigmag, g);
4220  }
4221 
4222  denoiseTVgradient(mu, xnew, origInput, mx, q, lambda, sigmag, g, grnew);
4223 
4224  s_norm = -100;
4225  p = 0;
4226  s2 = 0;
4227  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
4228  {
4229  xij = (DIRECT_A2D_ELEM(xnew,i,j) - DIRECT_A2D_ELEM(xold,i,j));
4230  p += xij * (DIRECT_A2D_ELEM(grnew,i,j) - DIRECT_A2D_ELEM(grold,i,j));
4231  s2 += xij * xij;
4232  if (xij > s_norm) s_norm = xij;
4233  }
4234 
4235  if (s_norm < tol)
4236  break;
4237 
4238  if (p <= 0)
4239  theta = thetamax;
4240  else
4241  theta = std::min(thetamax, std::max(thetamin, s2/p));
4242 
4243  denoiseTVproj(xnew, grnew, theta, dold);
4244 
4245  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xold)
4246  {
4247  DIRECT_A2D_ELEM(xold,i,j) = DIRECT_A2D_ELEM(xnew,i,j);
4248  DIRECT_A2D_ELEM(grold,i,j) = DIRECT_A2D_ELEM(grnew,i,j);
4249  }
4250 
4251  // we already computed it, let's use it for next round
4252  fold = fnew;
4253  }
4254 }
4255 
4256 void matlab_filter(Matrix1D<double> &B, Matrix1D<double> &A,
4257  const MultidimArray<double> &X, MultidimArray<double> &Y, MultidimArray<double> &Z)
4258 {
4259  if (VEC_ELEM(A,0)!=1.0)
4260  {
4261  A/=VEC_ELEM(A,0);
4262  B/=VEC_ELEM(B,0);
4263  }
4264 
4265  size_t input_size = XSIZE(X);
4266  size_t filter_order = std::max(VEC_XSIZE(A),VEC_XSIZE(B));
4267  if (VEC_XSIZE(B)!=filter_order)
4268  B.resize(filter_order);
4269  if (VEC_XSIZE(A)!=filter_order)
4270  A.resize(filter_order);
4271  Y.resize(input_size);
4272  Z.initZeros(filter_order);
4273 
4274  for (size_t i = 0; i < input_size; ++i)
4275  {
4276  size_t order = filter_order - 1;
4277  while (order)
4278  {
4279  if (i >= order)
4280  {
4281  DIRECT_A1D_ELEM(Z,order - 1) = VEC_ELEM(B,order) * DIRECT_A1D_ELEM(X,i - order) -
4282  VEC_ELEM(A,order) * DIRECT_A1D_ELEM(Y,i - order) + DIRECT_A1D_ELEM(Z,order);
4283  }
4284  --order;
4285  }
4286  DIRECT_A1D_ELEM(Y,i) = VEC_ELEM(B,0) * DIRECT_A1D_ELEM(X,i) + DIRECT_A1D_ELEM(Z,0);
4287  }
4288 }
4289 
4290 template <typename T>
4291 double mutualInformation(const MultidimArray< T >& x,
4292  const MultidimArray< T >& y,
4293  int nx,
4294  int ny,
4295  const MultidimArray< int >* mask)
4296 {
4297  SPEED_UP_temps;
4298 
4299  long n = 0;
4300  Histogram1D histx;
4301  Histogram1D histy;
4302  Histogram2D histxy;
4303  MultidimArray< T > aux_x;
4304  MultidimArray< T > aux_y;
4305  MultidimArray< double > mx;
4306  MultidimArray< double > my;
4307  MultidimArray< double > mxy;
4308  int xdim;
4309  int ydim;
4310  int zdim;
4311  double retval = 0.0;
4312 
4313  xdim=XSIZE(x);
4314  ydim=YSIZE(x);
4315  zdim=ZSIZE(x);
4316  aux_x.resize(xdim * ydim * zdim);
4317  xdim=XSIZE(y);
4318  ydim=YSIZE(y);
4319  zdim=ZSIZE(y);
4320  aux_y.resize(xdim * ydim * zdim);
4321 
4322  FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
4323  {
4324  if (mask != nullptr)
4325  if (!(*mask)(k, i, j))
4326  continue;
4327 
4328  aux_x(n) = A3D_ELEM(x, k, i, j);
4329  aux_y(n) = A3D_ELEM(y, k, i, j);
4330  n++;
4331  }
4332 
4333  aux_x.resize(n);
4334  aux_y.resize(n);
4335 
4336  if (n != 0)
4337  {
4338  if (nx == 0)
4339  //Assume Gaussian distribution
4340  nx = (int)((log((double) n) / LOG2) + 1);
4341 
4342  if (ny == 0)
4343  //Assume Gaussian distribution
4344  ny = (int)((log((double) n) / LOG2) + 1);
4345 
4346  compute_hist(aux_x, histx, nx);
4347  compute_hist(aux_y, histy, ny);
4348  compute_hist(aux_x, aux_y, histxy, nx, ny);
4349 
4350  mx = histx;
4351  my = histy;
4352  mxy = histxy;
4353  for (int i = 0; i < nx; i++)
4354  {
4355  double histxi = (histx(i)) / n;
4356  for (int j = 0; j < ny; j++)
4357  {
4358  double histyj = (histy(j)) / n;
4359  double histxyij = (histxy(i, j)) / n;
4360  if (histxyij > 0)
4361  retval += histxyij * log(histxyij / (histxi * histyj)) /
4362  LOG2;
4363  }
4364  }
4365 
4366  return retval;
4367  }
4368  else
4369  return 0;
4370 }
DenoiseTVFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3640
FOR_ALL_ELEMENTS_IN_MATRIX2D
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
NSIZE
#define NSIZE(v)
Definition: multidim_array_base.h:103
VECTOR_R2
#define VECTOR_R2(v, x, y)
Definition: matrix1d.h:112
DiffusionFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3725
Shah_energy
double Shah_energy(const MultidimArray< double > &img, const MultidimArray< double > &surface_strength, const MultidimArray< double > &edge_strength, double K, const Matrix1D< double > &W)
Definition: filters.cpp:2511
denoiseTVFilter
void denoiseTVFilter(MultidimArray< double > &xnew, int maxIter)
Definition: filters.cpp:4129
Histogram1D::clear
void clear()
Definition: histogram.cpp:40
centerImageRotationally
void centerImageRotationally(MultidimArray< double > &I, RotationalCorrelationAux &aux)
Definition: filters.cpp:3248
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
imedNormalizedDistance
double imedNormalizedDistance(const MultidimArray< double > &I1, const MultidimArray< double > &I2)
Definition: filters.cpp:1321
min
void min(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:118
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
CorrelationAux::transformer1
FourierTransformer transformer1
Definition: xmipp_fftw.h:555
set_DWT_type
void set_DWT_type(int DWT_type)
Definition: wavelet.cpp:154
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
giniCoeff
double giniCoeff(MultidimArray< double > &I, int varKernelSize)
Definition: filters.cpp:972
AlignmentAux::IauxRS
MultidimArray< double > IauxRS
Definition: filters.h:512
centerImageTranslationally
void centerImageTranslationally(MultidimArray< double > &I, CorrelationAux &aux)
Definition: filters.cpp:3212
nn
nn
Definition: numerical_recipes.cpp:4700
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
MultidimArray::threshold
void threshold(const String &type, T a, T b=0, MultidimArray< int > *mask=NULL)
Definition: multidim_array.h:3089
AlignmentTransforms::FFTI
MultidimArray< std::complex< double > > FFTI
Definition: filters.h:526
subtractColumnMeans
void subtractColumnMeans(Matrix2D< double > &A)
Definition: matrix2d.cpp:230
INITIAL_ROTATE_THRESHOLD
constexpr float INITIAL_ROTATE_THRESHOLD
Definition: filters.cpp:2045
LogFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3627
FourierTransformer::inverseFourierTransform
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
forcePositive
void forcePositive(MultidimArray< double > &V)
Definition: filters.cpp:3506
transformationMatrix2Parameters2D
void transformationMatrix2Parameters2D(const Matrix2D< T > &A, bool &flip, T &scale, T &shiftX, T &shiftY, T &psi)
Definition: transformations.cpp:272
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
MultidimArray::sort
void sort(MultidimArray< T > &result) const
Definition: multidim_array.cpp:1092
LogFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3635
fastBestRotationAroundY
double fastBestRotationAroundY(const MultidimArray< double > &IrefCyl, const MultidimArray< double > &I, CorrelationAux &aux, VolumeAlignmentAux &aux2)
Definition: filters.cpp:2290
FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
Definition: multidim_array_base.h:509
AlignmentAux::ASR
Matrix2D< double > ASR
Definition: filters.h:509
VolumeAlignmentAux::Icyl
MultidimArray< double > Icyl
Definition: filters.h:571
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
c
doublereal * c
Definition: numerical_recipes.cpp:2230
g
doublereal * g
Definition: numerical_recipes.cpp:2473
MULTIDIM_SIZE
#define MULTIDIM_SIZE(v)
Definition: multidim_array_base.h:115
interpolatedElementBSplineDiffZ
double interpolatedElementBSplineDiffZ(MultidimArray< double > &vol, double x, double y, double z, int SplineDegree)
Definition: transformations.cpp:1164
matrixOperation_AB
void matrixOperation_AB(const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
Definition: matrix2d.cpp:411
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
correlationMasked
double correlationMasked(const MultidimArray< double > &I1, const MultidimArray< double > &I2)
Definition: filters.cpp:1397
dAij
#define dAij(M, i, j)
Definition: multidim_array_base.h:345
OtsuSegmentation
double OtsuSegmentation(MultidimArray< double > &V)
Definition: filters.cpp:1028
BasisFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3796
rotationalInvariantMoments
void rotationalInvariantMoments(const MultidimArray< double > &img, const MultidimArray< int > *mask, MultidimArray< double > &v_out)
Definition: filters.cpp:2900
xl
doublereal * xl
Definition: numerical_recipes.cpp:2230
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
SGN0
#define SGN0(x)
Definition: xmipp_macros.h:169
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
BinaryCircularMask
void BinaryCircularMask(MultidimArray< int > &mask, double radius, int mode, double x0, double y0, double z0)
Definition: mask.cpp:193
AlignmentAux::IauxSR
MultidimArray< double > IauxSR
Definition: filters.h:511
Polar::getSampleNoOuterRing
int getSampleNoOuterRing() const
Definition: polar.h:386
BasisFilter::show
void show()
Definition: filters.cpp:3805
MultidimArray::computeStats
void computeStats(double &avg, double &stddev, T &minval, T &maxval) const
Definition: multidim_array.h:1791
y
static double * y
Definition: numerical_recipes.cpp:8487
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
BadPixelFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3600
MDL_SHIFT_X
Shift for the image in the X axis (double)
Definition: metadata_label.h:430
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
DIRECT_A2D_ELEM
#define DIRECT_A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:341
RetinexFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3844
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
A1D_ELEM
#define A1D_ELEM(v, i)
Definition: multidim_array_base.h:539
AlignmentAux::plans
Polar_fftw_plans * plans
Definition: filters.h:514
getImageSize
void getImageSize(const MetaData &md, size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim, MDLabel image_label)
Definition: metadata_extension.cpp:273
MultidimArray::computeAvgStdev
void computeAvgStdev(U &avg, U &stddev) const
Definition: multidim_array.h:1836
correlation_matrix
void correlation_matrix(const MultidimArray< double > &m1, const MultidimArray< double > &m2, MultidimArray< double > &R, CorrelationAux &aux, bool center)
Definition: xmipp_fftw.cpp:869
MultidimArray::maxIndex
void maxIndex(int &jmax) const
Definition: multidim_array.h:1557
icdf_gauss
double icdf_gauss(double p)
Definition: xmipp_funcs.cpp:212
MULTIDIM_ARRAY
#define MULTIDIM_ARRAY(v)
Definition: multidim_array_base.h:126
MultidimArray::rangeAdjust
void rangeAdjust(T minF, T maxF)
Definition: multidim_array.h:1982
fastBestRotationAroundZ
double fastBestRotationAroundZ(const MultidimArray< double > &IrefCyl, const MultidimArray< double > &I, CorrelationAux &aux, VolumeAlignmentAux &aux2)
Definition: filters.cpp:2276
TWOPI
#define TWOPI
Definition: xmipp_macros.h:111
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
w
doublereal * w
Definition: numerical_recipes.cpp:2477
MultidimArray::initConstant
void initConstant(T val)
Definition: multidim_array.h:2705
BadPixelFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3563
fastCorrentropy
double fastCorrentropy(const MultidimArray< double > &x, const MultidimArray< double > &y, double sigma, const GaussianInterpolator &G, const MultidimArray< int > &mask)
Definition: filters.cpp:1244
Histogram1D::percentil
double percentil(double percent_mass)
Definition: histogram.cpp:160
CHECK_POINT_DIST
#define CHECK_POINT_DIST(i, j, proposedDistance)
FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:322
DiffusionFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3750
DenoiseTVFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3648
ERR_MULTIDIM_SIZE
Incorrect MultidimArray size.
Definition: xmipp_error.h:174
FINISHINGZ
#define FINISHINGZ(v)
Definition: multidim_array_base.h:71
estimateGaussian2D
void estimateGaussian2D(const MultidimArray< double > &I, double &a, double &b, Matrix1D< double > &mu, Matrix2D< double > &sigma, bool estimateMu, int iterations)
Definition: filters.cpp:2390
MetaData::ids
virtual IdIteratorProxy< false > ids()
Definition: metadata_base.h:754
eps
cmache_1 eps
Definition: numerical_recipes.cpp:2273
MedianFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3704
fourierBesselDecomposition
void fourierBesselDecomposition(const MultidimArray< double > &img_in, double r2, int k1, int k2)
Definition: filters.cpp:2469
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
matrixOperation_AtA
void matrixOperation_AtA(const Matrix2D< double > &A, Matrix2D< double > &B)
Definition: matrix2d.cpp:436
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
bestNonwrappingShift
void bestNonwrappingShift(const MultidimArray< double > &I1, const MultidimArray< double > &I2, double &shiftX, double &shiftY, CorrelationAux &aux)
Definition: filters.cpp:1895
MetaData::getValue
virtual bool getValue(MDObject &mdValueOut, size_t id) const =0
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
selfScaleToSize
void selfScaleToSize(int SplineDegree, MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim)
Definition: transformations.cpp:736
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
rotate
#define rotate(a, i, j, k, l)
computeStats_within_binary_mask
void computeStats_within_binary_mask(const MultidimArray< T1 > &mask, const MultidimArray< T > &m, double &min_val, double &max_val, double &avg, double &stddev)
Definition: mask.h:799
d
doublereal * d
Definition: numerical_recipes.cpp:2230
FOR_ALL_ELEMENTS_IN_ARRAY2D
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
Definition: multidim_array_base.h:390
theta
double theta
Definition: numerical_recipes.cpp:6214
CHECK_POINT_3D
#define CHECK_POINT_3D(k, i, j)
rotation3DMatrix
void rotation3DMatrix(double ang, char axis, Matrix2D< double > &result, bool homogeneous)
Definition: transformations.cpp:426
mutualInformation
double mutualInformation(const MultidimArray< T > &x, const MultidimArray< T > &y, int nx, int ny, const MultidimArray< int > *mask)
Definition: filters.cpp:4291
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
centerImage
Matrix2D< double > centerImage(MultidimArray< double > &I, CorrelationAux &aux, RotationalCorrelationAux &aux2, int Niter, bool limitShift)
Definition: filters.cpp:3277
detectBackground
void detectBackground(const MultidimArray< double > &vol, MultidimArray< double > &mask, double alpha, double &final_mean)
Definition: filters.cpp:197
RetinexFilter::laplacian
void laplacian(const MultidimArray< double > &img, MultidimArray< std::complex< double > > &fimg, bool direct)
Definition: filters.cpp:3857
CLIP
#define CLIP(x, x0, xF)
Definition: xmipp_macros.h:260
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
threshold
void threshold(double *phi, unsigned long nvox, double limit)
Definition: lib_vwk.cpp:524
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
Update_edge_Shah
double Update_edge_Shah(MultidimArray< double > &img, MultidimArray< double > &surface_strength, MultidimArray< double > &edge_strength, double K, const Matrix1D< double > &W)
Definition: filters.cpp:2657
ROTATE_THRESHOLD
constexpr float ROTATE_THRESHOLD
Definition: filters.cpp:2043
varianceFilter
void varianceFilter(MultidimArray< double > &I, int kernelSize, bool relative)
Definition: filters.cpp:775
solve
void solve(const Matrix2D< T > &A, const Matrix2D< T > &b, Matrix2D< T > &result)
Definition: numerical_tools.cpp:2218
DIRECT_A1D_ELEM
#define DIRECT_A1D_ELEM(v, i)
Definition: multidim_array_base.h:526
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
best_rotation
double best_rotation(const Polar< std::complex< double > > &I1, const Polar< std::complex< double > > &I2, RotationalCorrelationAux &aux)
Definition: polar.cpp:212
b
doublereal * b
Definition: numerical_recipes.cpp:2230
closing3D
void closing3D(const MultidimArray< double > &in, MultidimArray< double > &out, int neig, int count, int size)
Definition: morphology.cpp:422
ntot
int ntot
Definition: image_find_center.cpp:57
Histogram1D::index2val
void index2val(double i, double &v) const
Definition: histogram.h:295
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
log
void log(Image< double > &op)
Definition: image_operate.cpp:228
MultidimArray::computeMin
T computeMin() const
Definition: multidim_array.h:1568
MultidimArray::window
void window(MultidimArray< T1 > &result, int n0, int z0, int y0, int x0, int nF, int zF, int yF, int xF, T1 init_value=0) const
Definition: multidim_array.h:564
computeAlignmentTransforms
void computeAlignmentTransforms(const MultidimArray< double > &I, AlignmentTransforms &ITransforms, AlignmentAux &aux, CorrelationAux &aux2)
Definition: filters.cpp:2035
y0
#define y0
axis
char axis
Definition: image_operate.cpp:257
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
DAUB12
#define DAUB12
Definition: wavelet.h:39
compute_hist
void compute_hist(const MultidimArrayGeneric &array, Histogram1D &hist, int no_steps)
Definition: histogram.cpp:572
lambda
double * lambda
Definition: numerical_recipes.cpp:4414
x0
#define x0
XX
#define XX(v)
Definition: matrix1d.h:85
VolumeAlignmentAux::I12
MultidimArray< double > I12
Definition: filters.h:574
FourierTransformer::setReal
void setReal(MultidimArray< double > &img)
Definition: xmipp_fftw.cpp:129
unnormalizedGaussian2D
double unnormalizedGaussian2D(const Matrix1D< double > &r, const Matrix1D< double > &mu, const Matrix2D< double > &sigmainv)
Definition: filters.cpp:2379
distanceTransform
void distanceTransform(const MultidimArray< int > &in, MultidimArray< int > &out, bool wrap)
Definition: filters.cpp:584
xmipp_program.h
MultidimArray::binarize
void binarize(double val=0, double accuracy=XMIPP_EQUAL_ACCURACY, MultidimArray< int > *mask=NULL)
Definition: multidim_array.h:3252
index
viol index
Definition: numerical_recipes.cpp:5197
type
viol type
Definition: numerical_recipes.cpp:5198
CEIL
#define CEIL(x)
Definition: xmipp_macros.h:225
LogFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3617
in
int in
Definition: image_find_center.cpp:56
DWT
void DWT(const MultidimArray< T > &v, MultidimArray< double > &result, int isign=1)
Definition: wavelet.h:83
correlationWeighted
double correlationWeighted(MultidimArray< double > &I1, MultidimArray< double > &I2)
Definition: filters.cpp:1457
matrixOperation_AtB
void matrixOperation_AtB(const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
Definition: matrix2d.cpp:475
mu1
int mu1
Definition: image_find_center.cpp:57
fourier_filter.h
EntropyOtsuSegmentation
double EntropyOtsuSegmentation(MultidimArray< double > &V, double percentil, bool binarizeVolume)
Definition: filters.cpp:1145
computeAvgStddev
void computeAvgStddev(const std::vector< T > &V, double &avg, double &stddev)
Definition: xmipp_funcs.h:358
XMIPP_EQUAL_ACCURACY
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
MDL_FLIP
Flip the image? (bool)
Definition: metadata_label.h:236
MultidimArrayBase::sameShape
bool sameShape(const MultidimArrayBase &op) const
Definition: multidim_array_base.h:847
VolumeAlignmentAux::corr
MultidimArray< double > corr
Definition: filters.h:572
INITIAL_SHIFT_THRESHOLD
constexpr double INITIAL_SHIFT_THRESHOLD
Definition: filters.cpp:2044
produceSplineCoefficients
void produceSplineCoefficients(int SplineDegree, MultidimArray< double > &coeffs, const MultidimArray< T > &V1)
Definition: transformations.cpp:41
Matrix1D::resize
void resize(size_t Xdim, bool copy=true)
Definition: matrix1d.h:410
RetinexFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3902
imedDistance
double imedDistance(const MultidimArray< double > &I1, const MultidimArray< double > &I2)
Definition: filters.cpp:1269
alignSetOfImages
void alignSetOfImages(MetaData &MD, MultidimArray< double > &Iavg, int Niter, bool considerMirror)
Definition: filters.cpp:2199
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
MultidimArray::write
void write(const FileName &fn) const
Definition: multidim_array.cpp:970
dx
double dx
Definition: numerical_recipes.cpp:6210
max
void max(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:128
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ERR_NUMERICAL
Error related to numerical calculation.
Definition: xmipp_error.h:179
AlignmentAux::polarFourierI
Polar< std::complex< double > > polarFourierI
Definition: filters.h:515
removeSmallComponents
void removeSmallComponents(MultidimArray< double > &I, int size, int neighbourhood)
Definition: filters.cpp:689
contrastEnhancement
void contrastEnhancement(Image< double > *I)
Definition: filters.cpp:327
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
AlignmentAux::rotationalCorr
MultidimArray< double > rotationalCorr
Definition: filters.h:513
MedianFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3717
regionGrowing2D
void regionGrowing2D(const MultidimArray< double > &I_in, MultidimArray< double > &I_out, int i, int j, float stop_colour, float filling_colour, bool less, int neighbourhood)
Definition: filters.cpp:340
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
z
double z
Definition: numerical_recipes.cpp:8490
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
GaussianInterpolator::getValue
double getValue(double x) const
Definition: numerical_tools.h:118
denoiseTVproj
void denoiseTVproj(const MultidimArray< double > &X, const MultidimArray< double > &Y, double theta, MultidimArray< double > &dold)
Definition: filters.cpp:4108
BasisFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3812
ROUND
#define ROUND(x)
Definition: xmipp_macros.h:210
matrixOperation_ABt
void matrixOperation_ABt(const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
Definition: matrix2d.cpp:462
computeEdges
void computeEdges(const MultidimArray< double > &vol, MultidimArray< double > &vol_edge)
Definition: filters.cpp:3517
MDL_MAXCC
Maximum cross-correlation for the image (double)
Definition: metadata_label.h:279
BackgroundFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3660
log10
void log10(Image< double > &op)
Definition: image_operate.cpp:237
vectorR2
Matrix1D< double > vectorR2(double x, double y)
Definition: matrix1d.cpp:884
dummy
double dummy
Definition: numerical_recipes.cpp:4635
CorrelationAux::FFT1
MultidimArray< std::complex< double > > FFT1
Definition: xmipp_fftw.h:554
xF
#define xF
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
NZYX_ELEM
#define NZYX_ELEM(v, l, k, i, j)
Definition: multidim_array_base.h:154
smoothingShah
void smoothingShah(MultidimArray< double > &img, MultidimArray< double > &surface_strength, MultidimArray< double > &edge_strength, const Matrix1D< double > &W, int OuterLoops, int InnerLoops, int RefinementLoops, bool adjust_range)
Definition: filters.cpp:2720
sort
void sort(struct DCEL_T *dcel)
Definition: sorting.cpp:18
DIRECT_A3D_ELEM
#define DIRECT_A3D_ELEM(v, k, i, j)
Definition: multidim_array_base.h:239
SPEED_UP_temps
#define SPEED_UP_temps
Definition: xmipp_macros.h:419
FourierTransformer::FourierTransform
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
FFT_IDX2DIGFREQ
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
Definition: cuda_gpu_reconstruct_fourier.cpp:382
BasisFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3788
VolumeAlignmentAux::I1
MultidimArray< double > I1
Definition: filters.h:573
Polar_fftw_plans
Definition: polar.h:52
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
MAXFLOAT
#define MAXFLOAT
Definition: data_types.h:47
svdcmp
void svdcmp(const Matrix2D< T > &a, Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v)
Definition: matrix2d.cpp:125
NEXT_POWER_OF_2
#define NEXT_POWER_OF_2(x)
Definition: xmipp_macros.h:374
rotation2DMatrix
void rotation2DMatrix(T ang, Matrix2D< T > &result, bool homogeneous)
Definition: transformations.cpp:361
j
#define j
Definition: numerical_recipes.cpp:2493
scaleToSize
void scaleToSize(int SplineDegree, MultidimArrayGeneric &V2, const MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim)
Definition: transformations.cpp:746
YY
#define YY(v)
Definition: matrix1d.h:93
pixelDesvFilter
void pixelDesvFilter(MultidimArray< T > &V, double thresFactor)
Definition: filters.h:1378
bestRotationAroundZ
double bestRotationAroundZ(const MultidimArray< double > &Iref, const MultidimArray< double > &I, CorrelationAux &aux, VolumeAlignmentAux &aux2)
Definition: filters.cpp:2356
m
int m
Definition: image_find_center.cpp:57
MultidimArrayBase::copyShape
void copyShape(const MultidimArrayBase &m)
Definition: multidim_array_base.cpp:186
polarFourierTransform< true >
template void polarFourierTransform< true >(MultidimArray< double > const &, Polar< double > &, Polar< std::complex< double > > &, bool, int, int, Polar_fftw_plans *&, int)
fastCorrelation
T fastCorrelation(const MultidimArray< T > &x, const MultidimArray< T > &y)
Definition: filters.h:328
CHECK_POINT
#define CHECK_POINT(i, j)
BackgroundFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3689
MetaData::setValue
bool setValue(const MDLabel label, const T &valueIn, size_t id)
Definition: metadata_base.h:367
DenoiseTVFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3654
regionGrowing3DEqualValue
void regionGrowing3DEqualValue(const MultidimArray< double > &V_in, MultidimArray< int > &V_out, int filling_value=0)
Definition: filters.cpp:499
matlab_filter
void matlab_filter(Matrix1D< double > &B, Matrix1D< double > &A, const MultidimArray< double > &X, MultidimArray< double > &Y, MultidimArray< double > &Z)
Definition: filters.cpp:4256
SHIFT_THRESHOLD
constexpr double SHIFT_THRESHOLD
Definition: filters.cpp:2042
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
SHAH_CONVERGENCE_THRESHOLD
constexpr double SHAH_CONVERGENCE_THRESHOLD
Definition: filters.cpp:2719
VolumeAlignmentAux::I123
MultidimArray< double > I123
Definition: filters.h:575
regionGrowing3D
void regionGrowing3D(const MultidimArray< double > &V_in, MultidimArray< double > &V_out, int k, int i, int j, float stop_colour, float filling_colour, bool less)
Definition: filters.cpp:412
Polar
Definition: polar.h:67
labelImage2D
int labelImage2D(const MultidimArray< double > &I, MultidimArray< double > &label, int neighbourhood)
Definition: filters.cpp:648
boundMedianFilter
void boundMedianFilter(MultidimArray< T > &V, const MultidimArray< char > &mask, int n=0)
Definition: filters.h:1309
fillBinaryObject
void fillBinaryObject(MultidimArray< double > &I, int neighbourhood)
Definition: filters.cpp:758
fastBestRotationAroundX
double fastBestRotationAroundX(const MultidimArray< double > &IrefCyl, const MultidimArray< double > &I, CorrelationAux &aux, VolumeAlignmentAux &aux2)
Definition: filters.cpp:2304
Update_surface_Shah
double Update_surface_Shah(MultidimArray< double > &img, MultidimArray< double > &surface_strength, MultidimArray< double > &edge_strength, const Matrix1D< double > &W)
Definition: filters.cpp:2577
translate
void translate(PDB *pdb_original, PDB *pdb_move, unsigned num_atoms, double x0, double y0, double z0)
Definition: lib_pwk.cpp:125
bestShiftRealSpace
double bestShiftRealSpace(const MultidimArray< double > &I1, MultidimArray< double > &I2, double &shiftX, double &shiftY, const MultidimArray< int > *mask, int maxShift, double shiftStep)
Definition: filters.cpp:1989
VolumeAlignmentAux::IrefCyl
MultidimArray< double > IrefCyl
Definition: filters.h:570
RotationalCorrelationAux::local_transformer
FourierTransformer local_transformer
Definition: polar.h:794
typeCast
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
MedianFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3712
logFilter
void logFilter(MultidimArray< T > &V, double a, double b, double c)
Definition: filters.h:1411
keepBiggestComponent
void keepBiggestComponent(MultidimArray< double > &I, double percentage, int neighbourhood)
Definition: filters.cpp:713
medianFilter3x3
void medianFilter3x3(MultidimArray< T > &m, MultidimArray< T > &out)
Definition: filters.h:1088
BackgroundFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3674
interpolatedElementBSplineDiffX
double interpolatedElementBSplineDiffX(MultidimArray< double > &vol, double x, double y, double z, int SplineDegree)
Definition: transformations.cpp:853
String
std::string String
Definition: xmipp_strings.h:34
alglib::psi
double psi(const double x)
Definition: specialfunctions.cpp:2492
BadPixelFilter::readParams
void readParams(XmippProgram *program)
Definition: filters.cpp:3579
svdCorrelation
double svdCorrelation(const MultidimArray< double > &I1, const MultidimArray< double > &I2, const MultidimArray< int > *mask)
Definition: filters.cpp:1535
substractBackgroundRollingBall
void substractBackgroundRollingBall(MultidimArray< double > &I, int radius)
Definition: filters.cpp:75
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
realWRAP
#define realWRAP(x, x0, xF)
Definition: xmipp_macros.h:304
mu
int mu
Definition: image_find_center.cpp:56
LOG2
constexpr double LOG2
Definition: filters.h:33
tomographicDiffusion
double tomographicDiffusion(MultidimArray< double > &V, const Matrix1D< double > &alpha, double lambda)
Definition: filters.cpp:2764
IDWT
void IDWT(const MultidimArray< double > &v, MultidimArray< double > &result)
Definition: wavelet.cpp:159
u
doublereal * u
Definition: numerical_recipes.cpp:2230
bestShift
template float bestShift(MultidimArray< float > &, float &, float &, const MultidimArray< int > *, int)
AlignmentAux::R
Matrix2D< double > R
Definition: filters.h:510
MultidimArray::copy
void copy(Matrix2D< T > &op1) const
Definition: multidim_array.cpp:56
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
forceDWTSparsity
void forceDWTSparsity(MultidimArray< double > &V, double eps)
Definition: filters.cpp:3539
fillTriangle
void fillTriangle(MultidimArray< double > &img, int *tx, int *ty, double color)
Definition: filters.cpp:3012
noisyZonesFilter
void noisyZonesFilter(MultidimArray< double > &I, int kernelSize)
Definition: filters.cpp:870
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
r2
float r2
Definition: image_find_center.cpp:54
REALGAUSSIAN
#define REALGAUSSIAN
Definition: fourier_filter.h:87
Matrix1D::sort
Matrix1D< T > sort() const
Definition: matrix1d.cpp:850
EntropySegmentation
double EntropySegmentation(MultidimArray< double > &V)
Definition: filters.cpp:1078
AlignmentAux::ARS
Matrix2D< double > ARS
Definition: filters.h:508
interpolatedElementBSplineDiffY
double interpolatedElementBSplineDiffY(MultidimArray< double > &vol, double x, double y, double z, int SplineDegree)
Definition: transformations.cpp:1008
AlignmentTransforms::polarFourierI
Polar< std::complex< double > > polarFourierI
Definition: filters.h:525
MDL_SHIFT_Y
Shift for the image in the Y axis (double)
Definition: metadata_label.h:434
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
MultidimArray::computeAvg
double computeAvg() const
Definition: multidim_array.h:1739
RetinexFilter::defineParams
static void defineParams(XmippProgram *program)
Definition: filters.cpp:3835
FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(V1, V2)
Definition: multidim_array_base.h:298
localThresholding
void localThresholding(MultidimArray< double > &img, double C, double dimLocal, MultidimArray< int > &result, MultidimArray< int > *mask)
Definition: filters.cpp:3164
labelImage3D
int labelImage3D(const MultidimArray< double > &V, MultidimArray< double > &label)
Definition: filters.cpp:669
PI
#define PI
Definition: tools.h:43
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331
metadata_extension.h
epsilon
double epsilon
Definition: numerical_recipes.cpp:5702
alignImagesConsideringMirrors
double alignImagesConsideringMirrors(const MultidimArray< double > &Iref, const AlignmentTransforms &IrefTransforms, MultidimArray< double > &I, Matrix2D< double > &M, AlignmentAux &aux, CorrelationAux &aux2, RotationalCorrelationAux &aux3, bool wrap, const MultidimArray< int > *mask)
Definition: filters.cpp:2150
fastBestRotation
double fastBestRotation(const MultidimArray< double > &IrefCyl, const MultidimArray< double > &Icyl, CorrelationAux &aux, VolumeAlignmentAux &aux2, double deltaAng)
Definition: filters.cpp:2255
yF
#define yF
substractBackgroundPlane
void substractBackgroundPlane(MultidimArray< double > &I)
Definition: filters.cpp:40
STARTINGZ
#define STARTINGZ(v)
Definition: multidim_array_base.h:67
denoiseTVgradient
void denoiseTVgradient(double mu, const MultidimArray< double > &X, const MultidimArray< double > &Y, double s, double q, double lambda, double sigmag, double g, MultidimArray< double > &gradientI)
Definition: filters.cpp:4017
n
int * n
Definition: numerical_recipes.cpp:2229
MDL_IMAGE
Name of an image (std::string)
Definition: metadata_label.h:246
ZZ
#define ZZ(v)
Definition: matrix1d.h:101
a
doublereal * a
Definition: numerical_recipes.cpp:2230
MultidimArray::sum
double sum() const
Definition: multidim_array.h:3405
denoiseTVenergy
double denoiseTVenergy(double mu, const MultidimArray< double > &X, const MultidimArray< double > &Y, double s, double q, double lambda, double sigmag, double g)
Definition: filters.cpp:3959
LOWPASS
#define LOWPASS
Definition: fourier_filter.h:77
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305
MultidimArray::indexSort
void indexSort(MultidimArray< int > &indx) const
Definition: multidim_array.cpp:882
xnew
double * xnew
Definition: numerical_recipes.cpp:7187
MultidimArray::statisticsAdjust
void statisticsAdjust(U avgF, U stddevF)
Definition: multidim_array.h:2114
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710
covarianceMatrix
void covarianceMatrix(const MultidimArray< double > &I, Matrix2D< double > &C)
Definition: filters.cpp:1582
volume_convertCartesianToCylindrical
void volume_convertCartesianToCylindrical(const MultidimArray< double > &in, MultidimArray< double > &out, double Rmin, double Rmax, double deltaR, float angMin, double angMax, float deltaAng, Matrix1D< double > &axis)
Definition: polar.cpp:305
MultidimArray::computeMax
T computeMax() const
Definition: multidim_array.h:1537
inertiaMoments
void inertiaMoments(const MultidimArray< double > &img, const MultidimArray< int > *mask, Matrix1D< double > &v_out, Matrix2D< double > &u)
Definition: filters.cpp:2970
xmipp_transformation::BSPLINE3
Definition: transformations_defines.h:32
multidim_array.h
DiffusionFilter::apply
void apply(MultidimArray< double > &img)
Definition: filters.cpp:3775
delta
double * delta
Definition: numerical_recipes.cpp:7597
Generated by   doxygen 1.8.13