Xmipp  v3.23.11-Nereus
symmetries.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 <stdio.h>
27 #include <limits>
28 #include <chrono>
29 #include <random>
30 #include "symmetries.h"
31 
32 // Symmetrize_crystal_vectors==========================================
33 //A note: Please realize that we are not repeating code here.
34 //The class SymList deals with symmetries when expressed in
35 //Cartesian space, that is the basis is orthonormal. Here
36 //we describe symmetries in the crystallographic way
37 //that is, the basis and the crystal vectors are the same.
38 //For same symmetries both representations are almost the same
39 //but in general they are rather different.
40 
41 //IMPORTANT: matrix orden should match the one used in "readSymmetryFile"
42 //if not the wrong angles are assigned to the different matrices
43 
44 void symmetrizeCrystalVectors(Matrix1D<double> &aint,
45  Matrix1D<double> &bint,
46  Matrix1D<double> &shift,
47  int space_group,
48  int sym_no,
49  const Matrix1D<double> &eprm_aint,
50  const Matrix1D<double> &eprm_bint)
51 {
52  //Notice that we should use R.inv and not R to relate eprm.aint and aint
53  shift.initZeros();//I think this init is OK even the vector dim=0
54  switch (space_group)
55  {
56  case(sym_undefined):
57  case sym_P1:
58  XX(aint) = XX(eprm_aint);
59  YY(aint) = YY(eprm_aint);
60  XX(bint) = XX(eprm_bint);
61  YY(bint) = YY(eprm_bint);
62  break;
63  case sym_P2: std::cerr << "\n Group P2 not implemented\n";
64  exit(1);
65  break;
66  case sym_P2_1: std::cerr << "\n Group P2_1 not implemented\n";
67  exit(1);
68  break;
69  case sym_C2: std::cerr << "\n Group C2 not implemented\n";
70  exit(1);
71  break;
72  case sym_P222: std::cerr << "\n Group P222 not implemented\n";
73  exit(1);
74  break;
75  case sym_P2_122:
76  switch (sym_no)
77  {
78  case -1:
79  XX(aint) = XX(eprm_aint);
80  YY(aint) = YY(eprm_aint);
81  XX(bint) = XX(eprm_bint);
82  YY(bint) = YY(eprm_bint);
83  break;
84  case 0:
85  XX(aint) = -XX(eprm_aint);
86  YY(aint) = -YY(eprm_aint);
87  XX(bint) = -XX(eprm_bint);
88  YY(bint) = -YY(eprm_bint);
89  break;
90  case 1:
91  XX(aint) = -XX(eprm_aint);
92  YY(aint) = +YY(eprm_aint);
93  XX(bint) = -XX(eprm_bint);
94  YY(bint) = +YY(eprm_bint);
95  VECTOR_R3(shift, 0.5, 0.0, 0.0);
96  break;
97  case 2:
98  XX(aint) = +XX(eprm_aint);
99  YY(aint) = -YY(eprm_aint);
100  XX(bint) = +XX(eprm_bint);
101  YY(bint) = -YY(eprm_bint);
102  VECTOR_R3(shift, 0.5, 0.0, 0.0);
103  break;
104  }//switch P2_122 end
105  break;
106  case sym_P22_12:
107  switch (sym_no)
108  {
109  case -1:
110  XX(aint) = XX(eprm_aint);
111  YY(aint) = YY(eprm_aint);
112  XX(bint) = XX(eprm_bint);
113  YY(bint) = YY(eprm_bint);
114  break;
115  case 0:
116  XX(aint) = -XX(eprm_aint);
117  YY(aint) = -YY(eprm_aint);
118  XX(bint) = -XX(eprm_bint);
119  YY(bint) = -YY(eprm_bint);
120  break;
121  case 1:
122  XX(aint) = XX(eprm_aint);
123  YY(aint) = -YY(eprm_aint);
124  XX(bint) = XX(eprm_bint);
125  YY(bint) = -YY(eprm_bint);
126  VECTOR_R3(shift, 0.0, 0.5, 0.0);
127  break;
128  case 2:
129  XX(aint) = -XX(eprm_aint);
130  YY(aint) = YY(eprm_aint);
131  XX(bint) = -XX(eprm_bint);
132  YY(bint) = YY(eprm_bint);
133  VECTOR_R3(shift, 0.0, 0.5, 0.0);
134  break;
135  }//switch P22_12 end
136  break;
137 
138  case sym_P22_12_1: std::cerr << "\n Group P22_12_1 not implemented\n";
139  exit(1);
140  break;
141  case sym_P4:
142  switch (sym_no)
143  {
144  case -1:
145  XX(aint) = XX(eprm_aint);
146  YY(aint) = YY(eprm_aint);
147  XX(bint) = XX(eprm_bint);
148  YY(bint) = YY(eprm_bint);
149  break;
150  case 0:
151  XX(aint) = -YY(eprm_aint);
152  YY(aint) = XX(eprm_aint);
153  XX(bint) = -YY(eprm_bint);
154  YY(bint) = XX(eprm_bint);
155  break;
156  case 1:
157  XX(aint) = -XX(eprm_aint);
158  YY(aint) = -YY(eprm_aint);
159  XX(bint) = -XX(eprm_bint);
160  YY(bint) = -YY(eprm_bint);
161  break;
162  case 2:
163  XX(aint) = YY(eprm_aint);
164  YY(aint) = -XX(eprm_aint);
165  XX(bint) = YY(eprm_bint);
166  YY(bint) = -XX(eprm_bint);
167  break;
168  }//switch P4 end
169  break;
170  case sym_P422: REPORT_ERROR(ERR_NOT_IMPLEMENTED, "Group P422 not implemented");
171  break;
172  case sym_P42_12:
173  switch (sym_no)
174  {
175  case -1:
176  XX(aint) = XX(eprm_aint);
177  YY(aint) = YY(eprm_aint);
178  XX(bint) = XX(eprm_bint);
179  YY(bint) = YY(eprm_bint);
180  break;
181  case 0:
182  XX(aint) = -XX(eprm_aint);
183  YY(aint) = -YY(eprm_aint);
184  XX(bint) = -XX(eprm_bint);
185  YY(bint) = -YY(eprm_bint);
186  break;
187  case 1:
188  XX(aint) = +YY(eprm_aint);
189  YY(aint) = +XX(eprm_aint);
190  XX(bint) = +YY(eprm_bint);
191  YY(bint) = +XX(eprm_bint);
192  break;
193  case 2:
194  XX(aint) = -YY(eprm_aint);
195  YY(aint) = -XX(eprm_aint);
196  XX(bint) = -YY(eprm_bint);
197  YY(bint) = -XX(eprm_bint);
198  break;
199  case 3:
200  XX(aint) = +YY(eprm_aint);
201  YY(aint) = -XX(eprm_aint);
202  XX(bint) = +YY(eprm_bint);
203  YY(bint) = -XX(eprm_bint);
204  VECTOR_R3(shift, 0.5, 0.5, 0);
205  break;
206  case 4:
207  XX(aint) = -YY(eprm_aint);
208  YY(aint) = +XX(eprm_aint);
209  XX(bint) = -YY(eprm_bint);
210  YY(bint) = +XX(eprm_bint);
211  VECTOR_R3(shift, 0.5, 0.5, 0);
212  break;
213  case 5:
214  XX(aint) = -XX(eprm_aint);
215  YY(aint) = +YY(eprm_aint);
216  XX(bint) = -XX(eprm_bint);
217  YY(bint) = +YY(eprm_bint);
218  VECTOR_R3(shift, 0.5, 0.5, 0);
219  break;
220  case 6:
221  XX(aint) = +XX(eprm_aint);
222  YY(aint) = -YY(eprm_aint);
223  XX(bint) = +XX(eprm_bint);
224  YY(bint) = -YY(eprm_bint);
225  VECTOR_R3(shift, 0.5, 0.5, 0);
226  break;
227  default:
228  std::cout << "\n Wrong symmetry number "
229  "in symmetrize_crystal_vectors, bye"
230  << std::endl;
231  exit(1);
232  break;
233 
234 
235  }//switch P4212 end
236  break;
237  case sym_P3: std::cerr << "\n Group P3 not implemented\n";
238  exit(1);
239  break;
240  case sym_P312: std::cerr << "\n Group P312 not implemented\n";
241  exit(1);
242  break;
243  case sym_P6:
244  switch (sym_no)
245  {
246  case -1:
247  XX(aint) = XX(eprm_aint);
248  YY(aint) = YY(eprm_aint);
249  XX(bint) = XX(eprm_bint);
250  YY(bint) = YY(eprm_bint);
251  break;
252  case 0:
253  XX(aint) = XX(eprm_aint) - YY(eprm_aint);
254  YY(aint) = XX(eprm_aint);
255  XX(bint) = XX(eprm_bint) - YY(eprm_bint);
256  YY(bint) = XX(eprm_bint);
257  break;
258  case 1:
259  XX(aint) = -YY(eprm_aint);
260  YY(aint) = XX(eprm_aint) - YY(eprm_aint);
261  XX(bint) = -YY(eprm_bint);
262  YY(bint) = XX(eprm_bint) - YY(eprm_bint);
263  break;
264  case 2:
265  XX(aint) = -XX(eprm_aint);
266  YY(aint) = -YY(eprm_aint);
267  XX(bint) = -XX(eprm_bint);
268  YY(bint) = -YY(eprm_bint);
269  break;
270  case 3:
271  XX(aint) = -XX(eprm_aint) + YY(eprm_aint);
272  YY(aint) = -XX(eprm_aint);
273  XX(bint) = -XX(eprm_bint) + YY(eprm_bint);
274  YY(bint) = -XX(eprm_bint);
275  break;
276  case 4:
277  XX(aint) = +YY(eprm_aint);
278  YY(aint) = -XX(eprm_aint) + YY(eprm_aint);
279  XX(bint) = +YY(eprm_bint);
280  YY(bint) = -XX(eprm_bint) + YY(eprm_bint);
281  break;
282  }//switch P6 end
283  break;
284 
285  case sym_P622: std::cerr << "\n Group P622 not implemented\n";
286  exit(1);
287  break;
288  }
289 
290 }//symmetrizeCrystalVectors end
291 
292 #define Symmetrize_Vol(X) {\
293  for (size_t i=0; i<vol_in.VolumesNo(); i++)\
294  X(vol_in(i),vol_in.grid(i),eprm_aint,eprm_bint,mask,i, \
295  grid_type);\
296  }
297 
298 // Symmetrize_crystal_volume==========================================
299 //IMPORTANT: matrix orden should match the one used in "readSymmetryFile"
300 //if not the wrong angles are assigned to the different matrices
301 void symmetrizeCrystalVolume(GridVolume &vol_in,
302  const Matrix1D<double> &eprm_aint,
303  const Matrix1D<double> &eprm_bint,
304  int eprm_space_group,
305  const MultidimArray<int> &mask, int grid_type)
306 {
307  //SO FAR ONLY THE GRID CENTERED IN 0,0,0 IS SYMMETRIZED, THE OTHER
308  //ONE SINCE REQUIRE INTERPOLATION IS IGNORED
309  switch (eprm_space_group)
310  {
311  case sym_undefined:
312  case sym_P1:
313  break;
314  case sym_P2: std::cerr << "\n Group P2 not implemented\n";
315  exit(1);
316  break;
317  case sym_P2_1: std::cerr << "\n Group P2_1 not implemented\n";
318  exit(1);
319  break;
320  case sym_C2: std::cerr << "\n Group C2 not implemented\n";
321  exit(1);
322  break;
323  case sym_P222: std::cerr << "\n Group P222 not implemented\n";
324  exit(1);
325  break;
326  case sym_P2_122:
327  Symmetrize_Vol(symmetry_P2_122)//already has ;
328  break;
329  case sym_P22_12:
330  Symmetrize_Vol(symmetry_P22_12)//already has ;
331  break;
332  case sym_P22_12_1: std::cerr << "\n Group P22_12_1 not implemented\n";
333  exit(1);
334  break;
335  case sym_P4:
336  Symmetrize_Vol(symmetry_P4)//already has ;
337  break;
338  case sym_P422: std::cerr << "\n Group P422 not implemented\n";
339  exit(1);
340  break;
341  case sym_P42_12:
342  Symmetrize_Vol(symmetry_P42_12)//already has ;
343  break;
344  case sym_P3: std::cerr << "\n Group P3 not implemented\n";
345  exit(1);
346  break;
347  case sym_P312: std::cerr << "\n Group P312 not implemented\n";
348  exit(1);
349  break;
350  case sym_P6:
351  Symmetrize_Vol(symmetry_P6)//already has ;
352  break;
353  case sym_P622: std::cerr << "\n Group P622 not implemented\n";
354  exit(1);
355  break;
356  }
357 
358 
359 }//symmetrizeCrystalVectors end
360 #define put_inside(j,j_min,j_max,jint) \
361  if( (j) < (j_min) ) { (j) = (j) + (jint);}\
362  else if( (j) > (j_max) ) { (j) = (j) - (jint);};
363 
364 /* Symmetrizes a simple grid with P2_122 symmetry--------------------------*/
365 
366 void symmetry_P2_122(Image<double> &vol, const SimpleGrid &grid,
367  const Matrix1D<double> &eprm_aint,
368  const Matrix1D<double> &eprm_bint,
369  const MultidimArray<int> &mask, int volume_no,
370  int grid_type)
371 {
372 
373  auto ZZ_lowest = (int) ZZ(grid.lowest);
374  int YY_lowest = STARTINGY(mask);
375  int XX_lowest = STARTINGX(mask);
376  auto ZZ_highest = (int) ZZ(grid.highest);
377  int YY_highest = FINISHINGY(mask);
378  int XX_highest = FINISHINGX(mask);
379 
380  //if there is an extra slice in the z direction there is no way
381  //to calculate the -z slice
382  if (ABS(ZZ_lowest) > ABS(ZZ_highest))
383  ZZ_lowest = -ZZ_highest;
384  else
385  ZZ_highest = ABS(ZZ_lowest);
386 
387  while (1)
388  {
389  if (mask(0, XX_lowest) == 0)
390  XX_lowest++;
391  else
392  break;
393  if (XX_lowest == XX_highest)
394  {
395  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
396  exit(0);
397  }
398  }
399  while (1)
400  {
401  if (mask(0, XX_highest) == 0)
402  XX_highest--;
403  else
404  break;
405  if (XX_lowest == XX_highest)
406  {
407  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
408  exit(0);
409  }
410  }
411  while (1)
412  {
413  if (mask(YY_lowest, 0) == 0)
414  YY_lowest++;
415  else
416  break;
417  if (YY_lowest == YY_highest)
418  {
419  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
420  exit(0);
421  }
422  }
423  while (1)
424  {
425  if (mask(YY_highest, 0) == 0)
426  YY_highest--;
427  else
428  break;
429  if (YY_lowest == YY_highest)
430  {
431  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
432  exit(0);
433  }
434 
435  }
436 
437 
438  int maxZ, maxY, maxX, minZ, minY, minX;
439  int x, y, z;
440 
441  int x0, y0, z0;
442  int x1, y1, z1;
443  int x2, y2, z2;
444 
445  int XXaint, YYbint;
446  XXaint = (int) XX(eprm_aint);
447  YYbint = (int)YY(eprm_bint);
448 
449  int XXaint_2;
450  XXaint_2 = XXaint / 2;
451  int xx, yy, zz;
452 
453  if (ABS(XX_lowest) > ABS(XX_highest))
454  {
455  minX = XX_lowest;
456  maxX = 0;
457  }
458  else
459  {
460  minX = 0;
461  maxX = XX_highest;
462  }
463  if (ABS(YY_lowest) > ABS(YY_highest))
464  {
465  minY = YY_lowest;
466  maxY = 0;
467  }
468  else
469  {
470  minY = 0;
471  maxY = YY_highest;
472  }
473  if (ABS(ZZ_lowest) > ABS(ZZ_highest))
474  {
475  minZ = ZZ_lowest;
476  maxZ = 0;
477  }
478  else
479  {
480  minZ = 0;
481  maxZ = ZZ_highest;
482  }
483 
484  //FCC non supported yet
485  if (volume_no == 1 && grid_type == FCC)
486  {
487  std::cerr << "\nSimetries using FCC not implemented\n";
488  exit(1);
489  }
490 
491  for (z = minZ;z <= maxZ;z++)
492  for (y = minY;y <= maxY;y++)
493  for (x = minX;x <= maxX;x++)
494  {
495  //sym=-1---------------------------------------------------------
496  if (!A2D_ELEM(mask, y, x) || z < ZZ_lowest || z > ZZ_highest)
497  continue;
498 
499  //sym=0 ---------------------------------------------------------
500  xx = -x;
501  yy = -y;
502  zz = z;
503  // xx = x; yy=-y; zz=-z;
504  if (volume_no == 1)
505  {
506  xx--;
507  yy--;
508  }
509  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
510  {
511  put_inside(xx, XX_lowest, XX_highest, XXaint)
512  put_inside(yy, YY_lowest, YY_highest, YYbint)
513  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
514  std::cerr << "ERROR in symmetry_P function"
515  << "after correction spot is still"
516  << "outside mask\a" << std::endl;
517  }
518  x0 = xx;
519  y0 = yy;
520  z0 = zz;
521 
522  //sym=1----------------------------------------------------------
523  xx = XXaint_2 - x;
524  yy = y;
525  zz = -z;
526  // xx = -x; yy= -y+YYbint_2; zz= -z;
527  if (volume_no == 1)
528  {
529  xx--;
530  zz--;
531  }
532  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
533  {
534  put_inside(xx, XX_lowest, XX_highest, XXaint)
535  put_inside(yy, YY_lowest, YY_highest, YYbint)
536  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
537  std::cerr << "ERROR in symmetry_P function"
538  << "after correction spot is still"
539  << "outside mask\a" << std::endl;
540  }//sym=1 end
541  x1 = xx;
542  y1 = yy;
543  z1 = zz;
544 
545  //sym=2----------------------------------------------------------
546  xx = XXaint_2 + x;
547  yy = -y;
548  zz = -z;
549  // xx = -x; yy= y+YYbint_2; zz= -z;
550  if (volume_no == 1)
551  {
552  yy--;
553  zz--;
554  }
555  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
556  {
557  put_inside(xx, XX_lowest, XX_highest, XXaint)
558  put_inside(yy, YY_lowest, YY_highest, YYbint)
559  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
560  std::cerr << "ERROR in symmetry_P function"
561  << "after correction spot is still"
562  << "outside mask\a" << std::endl;
563  }//sym=2 end
564  x2 = xx;
565  y2 = yy;
566  z2 = zz;
567 
568  //only the first simple grid center and the origen is the same
569  //point.
570  // if(volume_no==1)
571  // {
572  // switch (grid_type) {
573  // case FCC: std::cerr<< "\nSimetries using FCC not implemented\n";break;
574  // case BCC: x0--;y0--; z1--;
575  // x2--;y2--;z2--; y3--;
576  // x4--; x5--; z5--;
577  // y6--;z6--;
578  // break;
579  // case CC: break;
580  // }
581  // }
582 
583  VOLVOXEL(vol, z , y , x) = VOLVOXEL(vol, z0, y0, x0) =
584  VOLVOXEL(vol, z1, y1, x1) = VOLVOXEL(vol, z2, y2, x2) =
585  (VOLVOXEL(vol, z , y , x) + VOLVOXEL(vol, z0, y0, x0) +
586  VOLVOXEL(vol, z1, y1, x1) + VOLVOXEL(vol, z2, y2, x2)) / 4.0;
587  }//for end
588 }//symmetryP2_122 end
589 /* Symmetrizes a simple grid with P2_122 symmetry--------------------------*/
590 
591 void symmetry_P22_12(Image<double> &vol, const SimpleGrid &grid,
592  const Matrix1D<double> &eprm_aint,
593  const Matrix1D<double> &eprm_bint,
594  const MultidimArray<int> &mask, int volume_no,
595  int grid_type)
596 {
597 
598  auto ZZ_lowest = (int) ZZ(grid.lowest);
599  int YY_lowest = STARTINGY(mask);
600  int XX_lowest = STARTINGX(mask);
601  auto ZZ_highest = (int) ZZ(grid.highest);
602  int YY_highest = FINISHINGY(mask);
603  int XX_highest = FINISHINGX(mask);
604 
605  //if there is an extra slice in the z direction there is no way
606  //to calculate the -z slice
607  if (ABS(ZZ_lowest) > ABS(ZZ_highest))
608  ZZ_lowest = -ZZ_highest;
609  else
610  ZZ_highest = ABS(ZZ_lowest);
611 
612  while (1)
613  {
614  if (mask(0, XX_lowest) == 0)
615  XX_lowest++;
616  else
617  break;
618  if (XX_lowest == XX_highest)
619  {
620  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
621  exit(0);
622  }
623  }
624  while (1)
625  {
626  if (mask(0, XX_highest) == 0)
627  XX_highest--;
628  else
629  break;
630  if (XX_lowest == XX_highest)
631  {
632  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
633  exit(0);
634  }
635  }
636  while (1)
637  {
638  if (mask(YY_lowest, 0) == 0)
639  YY_lowest++;
640  else
641  break;
642  if (YY_lowest == YY_highest)
643  {
644  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
645  exit(0);
646  }
647  }
648  while (1)
649  {
650  if (mask(YY_highest, 0) == 0)
651  YY_highest--;
652  else
653  break;
654  if (YY_lowest == YY_highest)
655  {
656  std::cerr << "Error in symmetry_P2_122, while(1)" << std::endl;
657  exit(0);
658  }
659 
660  }
661 
662 
663  int maxZ, maxY, maxX, minZ, minY, minX;
664  int x, y, z;
665 
666  int x0, y0, z0;
667  int x1, y1, z1;
668  int x2, y2, z2;
669 
670  int XXaint, YYbint;
671  XXaint = (int) XX(eprm_aint);
672  YYbint = (int)YY(eprm_bint);
673 
674  int YYbint_2;
675  YYbint_2 = YYbint / 2;
676  int xx, yy, zz;
677 
678  if (ABS(XX_lowest) > ABS(XX_highest))
679  {
680  minX = XX_lowest;
681  maxX = 0;
682  }
683  else
684  {
685  minX = 0;
686  maxX = XX_highest;
687  }
688  if (ABS(YY_lowest) > ABS(YY_highest))
689  {
690  minY = YY_lowest;
691  maxY = 0;
692  }
693  else
694  {
695  minY = 0;
696  maxY = YY_highest;
697  }
698  if (ABS(ZZ_lowest) > ABS(ZZ_highest))
699  {
700  minZ = ZZ_lowest;
701  maxZ = 0;
702  }
703  else
704  {
705  minZ = 0;
706  maxZ = ZZ_highest;
707  }
708 
709  //FCC non supported yet
710  if (volume_no == 1 && grid_type == FCC)
711  {
712  std::cerr << "\nSimetries using FCC not implemented\n";
713  exit(1);
714  }
715 
716  for (z = minZ;z <= maxZ;z++)
717  for (y = minY;y <= maxY;y++)
718  for (x = minX;x <= maxX;x++)
719  {
720  //sym=-1---------------------------------------------------------
721  if (!A2D_ELEM(mask, y, x) || z < ZZ_lowest || z > ZZ_highest)
722  continue;
723 
724  //sym=0 ---------------------------------------------------------
725  xx = -x;
726  yy = -y;
727  zz = z;
728  // xx = x; yy=-y; zz=-z;
729  if (volume_no == 1)
730  {
731  xx--;
732  yy--;
733  }
734  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
735  {
736  put_inside(xx, XX_lowest, XX_highest, XXaint)
737  put_inside(yy, YY_lowest, YY_highest, YYbint)
738  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
739  std::cerr << "ERROR in symmetry_P function"
740  << "after correction spot is still"
741  << "outside mask\a" << std::endl;
742  }
743  x0 = xx;
744  y0 = yy;
745  z0 = zz;
746 
747  //sym=1----------------------------------------------------------
748  // xx = XXaint_2-x; yy= y; zz= -z;
749  xx = x;
750  yy = -y + YYbint_2;
751  zz = -z;
752  if (volume_no == 1)
753  {
754  yy--;
755  zz--;
756  }
757  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
758  {
759  put_inside(xx, XX_lowest, XX_highest, XXaint)
760  put_inside(yy, YY_lowest, YY_highest, YYbint)
761  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
762  std::cerr << "ERROR in symmetry_P function"
763  << "after correction spot is still"
764  << "outside mask\a" << std::endl;
765  }//sym=1 end
766  x1 = xx;
767  y1 = yy;
768  z1 = zz;
769 
770  //sym=2----------------------------------------------------------
771  // xx = XXaint_2+x; yy= -y; zz= -z;
772  xx = -x;
773  yy = y + YYbint_2;
774  zz = -z;
775  if (volume_no == 1)
776  {
777  xx--;
778  zz--;
779  }
780  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
781  {
782  put_inside(xx, XX_lowest, XX_highest, XXaint)
783  put_inside(yy, YY_lowest, YY_highest, YYbint)
784  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
785  std::cerr << "ERROR in symmetry_P function"
786  << "after correction spot is still"
787  << "outside mask\a" << std::endl;
788  }//sym=2 end
789  x2 = xx;
790  y2 = yy;
791  z2 = zz;
792 
793  //only the first simple grid center and the origen is the same
794  //point.
795  // if(volume_no==1)
796  // {
797  // switch (grid_type) {
798  // case FCC: std::cerr<< "\nSimetries using FCC not implemented\n";break;
799  // case BCC: x0--;y0--; z1--;
800  // x2--;y2--;z2--; y3--;
801  // x4--; x5--; z5--;
802  // y6--;z6--;
803  // break;
804  // case CC: break;
805  // }
806  // }
807 
808  VOLVOXEL(vol, z , y , x) = VOLVOXEL(vol, z0, y0, x0) =
809  VOLVOXEL(vol, z1, y1, x1) = VOLVOXEL(vol, z2, y2, x2) =
810  (VOLVOXEL(vol, z , y , x) + VOLVOXEL(vol, z0, y0, x0) +
811  VOLVOXEL(vol, z1, y1, x1) + VOLVOXEL(vol, z2, y2, x2)) / 4.0;
812  }//for end
813 }//symmetryP2_122 end
814 /* Symmetrizes a simple grid with P4 symmetry --------------------------*/
815 void symmetry_P4(Image<double> &vol, const SimpleGrid &grid,
816  const Matrix1D<double> &eprm_aint,
817  const Matrix1D<double> &eprm_bint,
818  const MultidimArray<int> &mask, int volume_no, int grid_type)
819 {
820  auto ZZ_lowest = (int) ZZ(grid.lowest);
821  int YY_lowest = STARTINGY(mask);
822  int XX_lowest = STARTINGX(mask);
823  auto ZZ_highest = (int) ZZ(grid.highest);
824  int YY_highest = FINISHINGY(mask);
825  int XX_highest = FINISHINGX(mask);
826 
827  while (1)
828  {
829  if (mask(0, XX_lowest) == 0)
830  XX_lowest++;
831  else
832  break;
833  if (XX_lowest == XX_highest)
834  {
835  std::cerr << "Error in symmetry_P4, while(1)" << std::endl;
836  exit(0);
837  }
838  }
839  while (1)
840  {
841  if (mask(0, XX_highest) == 0)
842  XX_highest--;
843  else
844  break;
845  if (XX_lowest == XX_highest)
846  {
847  std::cerr << "Error in symmetry_P4, while(1)" << std::endl;
848  exit(0);
849  }
850  }
851  while (1)
852  {
853  if (mask(YY_lowest, 0) == 0)
854  YY_lowest++;
855  else
856  break;
857  if (YY_lowest == YY_highest)
858  {
859  std::cerr << "Error in symmetry_P4, while(1)" << std::endl;
860  exit(0);
861  }
862  }
863  while (1)
864  {
865  if (mask(YY_highest, 0) == 0)
866  YY_highest--;
867  else
868  break;
869  if (YY_lowest == YY_highest)
870  {
871  std::cerr << "Error in symmetry_P4, while(1)" << std::endl;
872  exit(0);
873  }
874 
875  }
876 
877  // int ZZ_lowest =(int) ZZ(grid.lowest);
878  // int YY_lowest =MAX((int) YY(grid.lowest),STARTINGY(mask));
879  // int XX_lowest =MAX((int) XX(grid.lowest),STARTINGX(mask));
880  // int ZZ_highest=(int) ZZ(grid.highest);
881  // int YY_highest=MIN((int) YY(grid.highest),FINISHINGY(mask));
882  // int XX_highest=MIN((int) XX(grid.highest),FINISHINGX(mask));
883 
884  int maxZ, maxY, maxX, minZ, minY, minX;
885  int x, y, z;
886 
887  int x0, y0, z0;
888  int x1, y1, z1;
889  int x2, y2, z2;
890 
891  int XXaint, YYbint;
892  XXaint = (int) XX(eprm_aint);
893  YYbint = (int)YY(eprm_bint);
894 
895  int xx, yy, zz;
896 
897  if (ABS(XX_lowest) > ABS(XX_highest))
898  {
899  minX = XX_lowest;
900  maxX = 0;
901  }
902  else
903  {
904  minX = 0;
905  maxX = XX_highest;
906  }
907  if (ABS(YY_lowest) > ABS(YY_highest))
908  {
909  minY = YY_lowest;
910  maxY = 0;
911  }
912  else
913  {
914  minY = 0;
915  maxY = YY_highest;
916  }
917 
918  minZ = ZZ_lowest;
919  maxZ = ZZ_highest;
920  //FCC non supported yet
921  if (volume_no == 1 && grid_type == FCC)
922  {
923  std::cerr << "\nSimetries using FCC not implemented\n";
924  exit(1);
925  }
926  for (z = minZ;z <= maxZ;z++)
927  for (y = minY;y <= maxY;y++)
928  for (x = minX;x <= maxX;x++)
929  {
930  //sym=-1---------------------------------------------------------
931  if (!A2D_ELEM(mask, y, x) || z < ZZ_lowest || z > ZZ_highest)
932  continue;
933 
934  //sym=0 ---------------------------------------------------------
935  xx = -y;
936  yy = x;
937  zz = z;
938  //only the first simple grid center and the origen is the same
939  //point.
940  if (volume_no == 1)
941  {
942  xx--;
943  }
944  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
945  {
946  put_inside(xx, XX_lowest, XX_highest, XXaint)
947  put_inside(yy, YY_lowest, YY_highest, YYbint)
948  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
949  std::cerr << "ERROR in symmetry_P function"
950  << "after correction spot is still"
951  << "outside mask\a" << std::endl;
952  }
953  x0 = xx;
954  y0 = yy;
955  z0 = zz;
956 
957  //sym=1----------------------------------------------------------
958  xx = -x;
959  yy = -y;
960  zz = z;
961  if (volume_no == 1)
962  {
963  xx--;
964  yy--;
965  }
966  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
967  {
968  put_inside(xx, XX_lowest, XX_highest, XXaint)
969  put_inside(yy, YY_lowest, YY_highest, YYbint)
970  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
971  std::cerr << "ERROR in symmetry_P function"
972  << "after correction spot is still"
973  << "outside mask\a" << std::endl;
974  }//sym=1 end
975  x1 = xx;
976  y1 = yy;
977  z1 = zz;
978 
979  //sym=2----------------------------------------------------------
980  xx = y;
981  yy = -x;
982  zz = z;
983  if (volume_no == 1)
984  {
985  yy--;
986  }
987  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
988  {
989  put_inside(xx, XX_lowest, XX_highest, XXaint)
990  put_inside(yy, YY_lowest, YY_highest, YYbint)
991  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
992  std::cerr << "ERROR in symmetry_P function"
993  << "after correction spot is still"
994  << "outside mask\a" << std::endl;
995  }//sym=2 end
996  x2 = xx;
997  y2 = yy;
998  z2 = zz;
999 
1000  VOLVOXEL(vol, z , y , x) = VOLVOXEL(vol, z0, y0, x0) =
1001  VOLVOXEL(vol, z1, y1, x1) = VOLVOXEL(vol, z2, y2, x2) =
1002  (VOLVOXEL(vol, z , y , x) + VOLVOXEL(vol, z0, y0, x0) +
1003  VOLVOXEL(vol, z1, y1, x1) + VOLVOXEL(vol, z2, y2, x2)) / 4.0;
1004  }//for end
1005 
1006 }
1007 /* Symmetrizes a simple grid with P4212 symmetry--------------------------*/
1008 
1009 void symmetry_P42_12(Image<double> &vol, const SimpleGrid &grid,
1010  const Matrix1D<double> &eprm_aint,
1011  const Matrix1D<double> &eprm_bint,
1012  const MultidimArray<int> &mask, int volume_no,
1013  int grid_type)
1014 {
1015 
1016  auto ZZ_lowest = (int) ZZ(grid.lowest);
1017  int YY_lowest = STARTINGY(mask);
1018  int XX_lowest = STARTINGX(mask);
1019  auto ZZ_highest = (int) ZZ(grid.highest);
1020  int YY_highest = FINISHINGY(mask);
1021  int XX_highest = FINISHINGX(mask);
1022 
1023  //if there is an extra slice in the z direction there is no way
1024  //to calculate the -z slice
1025  if (ABS(ZZ_lowest) > ABS(ZZ_highest))
1026  ZZ_lowest = -ZZ_highest;
1027  else
1028  ZZ_highest = ABS(ZZ_lowest);
1029 
1030  while (1)
1031  {
1032  if (mask(0, XX_lowest) == 0)
1033  XX_lowest++;
1034  else
1035  break;
1036  if (XX_lowest == XX_highest)
1037  {
1038  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1039  exit(0);
1040  }
1041  }
1042  while (1)
1043  {
1044  if (mask(0, XX_highest) == 0)
1045  XX_highest--;
1046  else
1047  break;
1048  if (XX_lowest == XX_highest)
1049  {
1050  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1051  exit(0);
1052  }
1053  }
1054  while (1)
1055  {
1056  if (mask(YY_lowest, 0) == 0)
1057  YY_lowest++;
1058  else
1059  break;
1060  if (YY_lowest == YY_highest)
1061  {
1062  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1063  exit(0);
1064  }
1065  }
1066  while (1)
1067  {
1068  if (mask(YY_highest, 0) == 0)
1069  YY_highest--;
1070  else
1071  break;
1072  if (YY_lowest == YY_highest)
1073  {
1074  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1075  exit(0);
1076  }
1077 
1078  }
1079 
1080  // int ZZ_lowest =(int) ZZ(grid.lowest);
1081  // int YY_lowest =MAX((int) YY(grid.lowest),STARTINGY(mask));
1082  // int XX_lowest =MAX((int) XX(grid.lowest),STARTINGX(mask));
1083  // int ZZ_highest=(int) ZZ(grid.highest);
1084  // int YY_highest=MIN((int) YY(grid.highest),FINISHINGY(mask));
1085  // int XX_highest=MIN((int) XX(grid.highest),FINISHINGX(mask));
1086 
1087  int maxZ, maxY, maxX, minZ, minY, minX;
1088  int x, y, z;
1089 
1090  int x0, y0, z0;
1091  int x1, y1, z1;
1092  int x2, y2, z2;
1093  int x3, y3, z3;
1094  int x4, y4, z4;
1095  int x5, y5, z5;
1096  int x6, y6, z6;
1097 
1098  int XXaint, YYbint;
1099  XXaint = (int) XX(eprm_aint);
1100  YYbint = (int)YY(eprm_bint);
1101 
1102  int XXaint_2, YYbint_2;
1103  XXaint_2 = XXaint / 2;
1104  YYbint_2 = YYbint / 2;
1105  int xx, yy, zz;
1106 
1107  if (ABS(XX_lowest) > ABS(XX_highest))
1108  {
1109  minX = XX_lowest;
1110  maxX = 0;
1111  }
1112  else
1113  {
1114  minX = 0;
1115  maxX = XX_highest;
1116  }
1117  if (ABS(YY_lowest) > ABS(YY_highest))
1118  {
1119  minY = YY_lowest;
1120  maxY = 0;
1121  }
1122  else
1123  {
1124  minY = 0;
1125  maxY = YY_highest;
1126  }
1127  if (ABS(ZZ_lowest) > ABS(ZZ_highest))
1128  {
1129  minZ = ZZ_lowest;
1130  maxZ = 0;
1131  }
1132  else
1133  {
1134  minZ = 0;
1135  maxZ = ZZ_highest;
1136  }
1137 
1138  //FCC non supported yet
1139  if (volume_no == 1 && grid_type == FCC)
1140  {
1141  std::cerr << "\nSimetries using FCC not implemented\n";
1142  exit(1);
1143  }
1144 
1145  for (z = minZ;z <= maxZ;z++)
1146  for (y = minY;y <= maxY;y++)
1147  for (x = minX;x <= maxX;x++)
1148  {
1149  //sym=-1---------------------------------------------------------
1150  if (!A2D_ELEM(mask, y, x) || z < ZZ_lowest || z > ZZ_highest)
1151  continue;
1152 
1153  //sym=0 ---------------------------------------------------------
1154  xx = -x;
1155  yy = -y;
1156  zz = z;
1157  if (volume_no == 1)
1158  {
1159  xx--;
1160  yy--;
1161  }
1162  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1163  {
1164  put_inside(xx, XX_lowest, XX_highest, XXaint)
1165  put_inside(yy, YY_lowest, YY_highest, YYbint)
1166  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1167  std::cerr << "ERROR in symmetry_P function"
1168  << "after correction spot is still"
1169  << "outside mask\a" << std::endl;
1170  }
1171  x0 = xx;
1172  y0 = yy;
1173  z0 = zz;
1174 
1175  //sym=1----------------------------------------------------------
1176  xx = y;
1177  yy = x;
1178  zz = -z;//I think z-- is always inside the grid
1179  //we do not need to check
1180  if (volume_no == 1)
1181  {
1182  zz--;
1183  }
1184  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1185  {
1186  put_inside(xx, XX_lowest, XX_highest, XXaint)
1187  put_inside(yy, YY_lowest, YY_highest, YYbint)
1188  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1189  std::cerr << "ERROR in symmetry_P function"
1190  << "after correction spot is still"
1191  << "outside mask\a" << std::endl;
1192  }//sym=1 end
1193  x1 = xx;
1194  y1 = yy;
1195  z1 = zz;
1196 
1197  //sym=2----------------------------------------------------------
1198  xx = -y;
1199  yy = -x;
1200  zz = -z;
1201  if (volume_no == 1)
1202  {
1203  xx--;
1204  yy--;
1205  zz--;
1206  }
1207  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1208  {
1209  put_inside(xx, XX_lowest, XX_highest, XXaint)
1210  put_inside(yy, YY_lowest, YY_highest, YYbint)
1211  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1212  std::cerr << "ERROR in symmetry_P function"
1213  << "after correction spot is still"
1214  << "outside mask\a" << std::endl;
1215  }//sym=2 end
1216  x2 = xx;
1217  y2 = yy;
1218  z2 = zz;
1219 
1220  //sym=3----------------------------------------------------------
1221  xx = y + XXaint_2;
1222  yy = -x + YYbint_2;
1223  zz = z;
1224  if (volume_no == 1)
1225  {
1226  yy--;
1227  }
1228  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1229  {
1230  put_inside(xx, XX_lowest, XX_highest, XXaint)
1231  put_inside(yy, YY_lowest, YY_highest, YYbint)
1232  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1233  std::cerr << "ERROR in symmetry_P function"
1234  << "after correction spot is still"
1235  << "outside mask\a" << std::endl;
1236  }//sym=3 end
1237  x3 = xx;
1238  y3 = yy;
1239  z3 = zz;
1240 
1241  //sym=4----------------------------------------------------------
1242  xx = -y + XXaint_2;
1243  yy = + x + YYbint_2;
1244  zz = z;
1245  if (volume_no == 1)
1246  {
1247  xx--;
1248  }
1249  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1250  {
1251  put_inside(xx, XX_lowest, XX_highest, XXaint)
1252  put_inside(yy, YY_lowest, YY_highest, YYbint)
1253  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1254  std::cerr << "ERROR in symmetry_P function"
1255  << "after correction spot is still"
1256  << "outside mask\a" << std::endl;
1257  }//sym=4 end
1258  x4 = xx;
1259  y4 = yy;
1260  z4 = zz;
1261  //sym=5----------------------------------------------------------
1262  xx = -x + XXaint_2;
1263  yy = + y + YYbint_2;
1264  zz = -z;
1265  if (volume_no == 1)
1266  {
1267  xx--;
1268  zz--;
1269  }
1270  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1271  {
1272  put_inside(xx, XX_lowest, XX_highest, XXaint)
1273  put_inside(yy, YY_lowest, YY_highest, YYbint)
1274  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1275  std::cerr << "ERROR in symmetry_P function"
1276  << "after correction spot is still"
1277  << "outside mask\a" << std::endl;
1278  }//sym=5 end
1279  x5 = xx;
1280  y5 = yy;
1281  z5 = zz;
1282 
1283  //sym=6----------------------------------------------------------
1284  xx = + x + XXaint_2;
1285  yy = -y + YYbint_2;
1286  zz = -z;
1287  if (volume_no == 1)
1288  {
1289  yy--;
1290  zz--;
1291  }
1292  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1293  {
1294  put_inside(xx, XX_lowest, XX_highest, XXaint)
1295  put_inside(yy, YY_lowest, YY_highest, YYbint)
1296  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1297  std::cerr << "ERROR in symmetry_P function"
1298  << "after correction spot is still"
1299  << "outside mask\a" << std::endl;
1300  }//sym=6 end
1301  x6 = xx;
1302  y6 = yy;
1303  z6 = zz;
1304 
1305  //only the first simple grid center and the origen is the same
1306  //point.
1307  // if(volume_no==1)
1308  // {
1309  // switch (grid_type) {
1310  // case FCC: std::cerr<< "\nSimetries using FCC not implemented\n";break;
1311  // case BCC: x0--;y0--; z1--;
1312  // x2--;y2--;z2--; y3--;
1313  // x4--; x5--; z5--;
1314  // y6--;z6--;
1315  // break;
1316  // case CC: break;
1317  // }
1318  // }
1319 
1320  VOLVOXEL(vol, z , y , x) = VOLVOXEL(vol, z0, y0, x0) =
1321  VOLVOXEL(vol, z1, y1, x1) = VOLVOXEL(vol, z2, y2, x2) =
1322  VOLVOXEL(vol, z3, y3, x3) = VOLVOXEL(vol, z4, y4, x4) =
1323  VOLVOXEL(vol, z5, y5, x5) = VOLVOXEL(vol, z6, y6, x6) =
1324  (VOLVOXEL(vol, z , y , x) + VOLVOXEL(vol, z0, y0, x0) +
1325  VOLVOXEL(vol, z1, y1, x1) + VOLVOXEL(vol, z2, y2, x2) +
1326  VOLVOXEL(vol, z3, y3, x3) + VOLVOXEL(vol, z4, y4, x4) +
1327  VOLVOXEL(vol, z5, y5, x5) + VOLVOXEL(vol, z6, y6, x6)) / 8.0;
1328  }//for end
1329 }//symmetryP42_12 end
1330 /* Symmetrizes a simple grid with P6 symmetry-----------------------------*/
1331 void symmetry_P6(Image<double> &vol, const SimpleGrid &grid,
1332  const Matrix1D<double> &eprm_aint,
1333  const Matrix1D<double> &eprm_bint,
1334  const MultidimArray<int> &mask, int volume_no,
1335  int grid_type)
1336 {
1337 
1338  auto ZZ_lowest = (int) ZZ(grid.lowest);
1339  int YY_lowest = STARTINGY(mask);
1340  int XX_lowest = STARTINGX(mask);
1341  auto ZZ_highest = (int) ZZ(grid.highest);
1342  int YY_highest = FINISHINGY(mask);
1343  int XX_highest = FINISHINGX(mask);
1344 
1345 
1346  while (1)
1347  {
1348  if (mask(0, XX_lowest) == 0)
1349  XX_lowest++;
1350  else
1351  break;
1352  if (XX_lowest == XX_highest)
1353  {
1354  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1355  exit(0);
1356  }
1357  }
1358  while (1)
1359  {
1360  if (mask(0, XX_highest) == 0)
1361  XX_highest--;
1362  else
1363  break;
1364  if (XX_lowest == XX_highest)
1365  {
1366  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1367  exit(0);
1368  }
1369  }
1370  while (1)
1371  {
1372  if (mask(YY_lowest, 0) == 0)
1373  YY_lowest++;
1374  else
1375  break;
1376  if (YY_lowest == YY_highest)
1377  {
1378  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1379  exit(0);
1380  }
1381  }
1382  while (1)
1383  {
1384  if (mask(YY_highest, 0) == 0)
1385  YY_highest--;
1386  else
1387  break;
1388  if (YY_lowest == YY_highest)
1389  {
1390  std::cerr << "Error in symmetry_P42_12, while(1)" << std::endl;
1391  exit(0);
1392  }
1393 
1394  }
1395 
1396  // int ZZ_lowest =(int) ZZ(grid.lowest);
1397  // int YY_lowest =MAX((int) YY(grid.lowest),STARTINGY(mask));
1398  // int XX_lowest =MAX((int) XX(grid.lowest),STARTINGX(mask));
1399  // int ZZ_highest=(int) ZZ(grid.highest);
1400  // int YY_highest=MIN((int) YY(grid.highest),FINISHINGY(mask));
1401  // int XX_highest=MIN((int) XX(grid.highest),FINISHINGX(mask));
1402 
1403  int maxZ, maxY, maxX, minZ, minY, minX;
1404  int x, y, z;
1405 
1406  int x0, y0, z0;
1407  int x1, y1, z1;
1408  int x2, y2, z2;
1409  int x3, y3, z3;
1410  int x4, y4, z4;
1411 
1412  int XXaint, YYbint;
1413  XXaint = (int) XX(eprm_aint);
1414  YYbint = (int)YY(eprm_bint);
1415 
1416  int xx, yy, zz;
1417 
1418  if (ABS(XX_lowest) > ABS(XX_highest))
1419  {
1420  minX = XX_lowest;
1421  maxX = 0;
1422  }
1423  else
1424  {
1425  minX = 0;
1426  maxX = XX_highest;
1427  }
1428  //P6 is tricky. I have decide to apply it to half the volume
1429  //instead of to 1 sizth. I think the amount of ifs that I save
1430  //are worth this larger loop
1431 
1432  minY = YY_lowest;
1433  maxY = YY_highest;
1434 
1435  minZ = ZZ_lowest;
1436  maxZ = ZZ_highest;
1437 
1438  //FCC non supported yet
1439  if (volume_no == 1 && grid_type == FCC)
1440  {
1441  std::cerr << "\nSimetries using FCC not implemented\n";
1442  exit(1);
1443  }
1444 
1445  for (z = minZ;z <= maxZ;z++)
1446  for (y = minY;y <= maxY;y++)
1447  for (x = minX;x <= maxX;x++)
1448  {
1449  //sym=-1---------------------------------------------------------
1450  if (!A2D_ELEM(mask, y, x) || z < ZZ_lowest || z > ZZ_highest)
1451  continue;
1452 
1453  //sym=0 ---------------------------------------------------------
1454  xx = x - y;
1455  yy = x;
1456  zz = z;
1457  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1458  {
1459  put_inside(xx, XX_lowest, XX_highest, XXaint)
1460  put_inside(yy, YY_lowest, YY_highest, YYbint)
1461  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1462  std::cerr << "ERROR in symmetry_P function"
1463  << "after correction spot is still"
1464  << "outside mask\a" << std::endl;
1465  }
1466  x0 = xx;
1467  y0 = yy;
1468  z0 = zz;
1469 
1470  //sym=1----------------------------------------------------------
1471  xx = -y;
1472  yy = x - y;
1473  zz = z;
1474  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1475  {
1476  put_inside(xx, XX_lowest, XX_highest, XXaint)
1477  put_inside(yy, YY_lowest, YY_highest, YYbint)
1478  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1479  std::cerr << "ERROR in symmetry_P function"
1480  << "after correction spot is still"
1481  << "outside mask\a" << std::endl;
1482  }//sym=1 end
1483  x1 = xx;
1484  y1 = yy;
1485  z1 = zz;
1486 
1487  //sym=2----------------------------------------------------------
1488  xx = -x;
1489  yy = -y;
1490  zz = z;
1491  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1492  {
1493  put_inside(xx, XX_lowest, XX_highest, XXaint)
1494  put_inside(yy, YY_lowest, YY_highest, YYbint)
1495  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1496  std::cerr << "ERROR in symmetry_P function"
1497  << "after correction spot is still"
1498  << "outside mask\a" << std::endl;
1499  }//sym=2 end
1500  x2 = xx;
1501  y2 = yy;
1502  z2 = zz;
1503 
1504  //sym=3----------------------------------------------------------
1505  xx = -x + y;
1506  yy = -x;
1507  zz = z;
1508  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1509  {
1510  put_inside(xx, XX_lowest, XX_highest, XXaint)
1511  put_inside(yy, YY_lowest, YY_highest, YYbint)
1512  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1513  std::cerr << "ERROR in symmetry_P function"
1514  << "after correction spot is still"
1515  << "outside mask\a" << std::endl;
1516  }//sym=3 end
1517  x3 = xx;
1518  y3 = yy;
1519  z3 = zz;
1520 
1521  //sym=4----------------------------------------------------------
1522  xx = + y;
1523  yy = -x + y;
1524  zz = z;
1525  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1526  {
1527  put_inside(xx, XX_lowest, XX_highest, XXaint)
1528  put_inside(yy, YY_lowest, YY_highest, YYbint)
1529  if (!A2D_ELEM(mask, yy, xx) || mask.outside(yy, xx))
1530  std::cerr << "ERROR in symmetry_P function"
1531  << "after correction spot is still"
1532  << "outside mask\a" << std::endl;
1533  }//sym=4 end
1534  x4 = xx;
1535  y4 = yy;
1536  z4 = zz;
1537 
1538  if (volume_no == 1)
1539  {
1540  switch (grid_type)
1541  {
1542  case FCC:
1543  std::cerr << "\nSimetries using FCC not implemented\n";
1544  break;
1545  //there is no way to reinforce P6 in the second grid without
1546  //interpolation. This is the best we can do.
1547  case BCC:
1548  x1 = x0 = x;
1549  y1 = y0 = y;
1550  x2 = x3 = x4 = -x - 1;
1551  y2 = y3 = y4 = -y - 1;
1552  break;
1553  case CC:
1554  break;
1555  }
1556  }
1557 
1558 
1559  VOLVOXEL(vol, z , y , x) = VOLVOXEL(vol, z0, y0, x0) =
1560  VOLVOXEL(vol, z1, y1, x1) = VOLVOXEL(vol, z2, y2, x2) =
1561  VOLVOXEL(vol, z3, y3, x3) = VOLVOXEL(vol, z4, y4, x4) =
1562  (VOLVOXEL(vol, z , y , x) + VOLVOXEL(vol, z0, y0, x0) +
1563  VOLVOXEL(vol, z1, y1, x1) + VOLVOXEL(vol, z2, y2, x2) +
1564  VOLVOXEL(vol, z3, y3, x3) + VOLVOXEL(vol, z4, y4, x4)) / 6.0;
1565 
1566  }//for end
1567 
1568 }
1569 #undef wrap_as_Crystal
1570 #undef DEBUG
1571 
1572 
1573 // Forward declaration
1574 double interpolatedElement3DHelical(const MultidimArray<double> &Vin, double x, double y, double z, double zHelical,
1575  double sinRotHelical, double cosRotHelical);
1576 
1577 double interpolatedElement3DHelicalInt(const MultidimArray<double> &Vin, int x, int y, int z, double zHelical,
1578  double sinRotHelical, double cosRotHelical)
1579 {
1580  if (x<STARTINGX(Vin) || x>FINISHINGX(Vin) || y<STARTINGY(Vin) || y>FINISHINGY(Vin))
1581  return 0.0;
1582  if (z>=STARTINGZ(Vin) && z<=FINISHINGZ(Vin))
1583  return A3D_ELEM(Vin,z,y,x);
1584  else if (z<STARTINGZ(Vin))
1585  {
1586  double newx=cosRotHelical*x-sinRotHelical*y;
1587  double newy=sinRotHelical*x+cosRotHelical*y;
1588  return interpolatedElement3DHelical(Vin,newx,newy,z+zHelical,zHelical, sinRotHelical, cosRotHelical);
1589  }
1590  else
1591  {
1592  double newx=cosRotHelical*x+sinRotHelical*y;
1593  double newy=-sinRotHelical*x+cosRotHelical*y;
1594  return interpolatedElement3DHelical(Vin,newx,newy,z-zHelical,zHelical, sinRotHelical, cosRotHelical);
1595  }
1596 }
1597 
1598 double interpolatedElement3DHelical(const MultidimArray<double> &Vin, double x, double y, double z, double zHelical,
1599  double sinRotHelical, double cosRotHelical)
1600 {
1601  int x0 = floor(x);
1602  double fx = x - x0;
1603  int x1 = x0 + 1;
1604 
1605  int y0 = floor(y);
1606  double fy = y - y0;
1607  int y1 = y0 + 1;
1608 
1609  int z0 = floor(z);
1610  double fz = z - z0;
1611  int z1 = z0 + 1;
1612 
1613  double d000 = interpolatedElement3DHelicalInt(Vin, x0, y0, z0, zHelical, sinRotHelical, cosRotHelical);
1614  double d001 = interpolatedElement3DHelicalInt(Vin, x1, y0, z0, zHelical, sinRotHelical, cosRotHelical);
1615  double d010 = interpolatedElement3DHelicalInt(Vin, x0, y1, z0, zHelical, sinRotHelical, cosRotHelical);
1616  double d011 = interpolatedElement3DHelicalInt(Vin, x1, y1, z0, zHelical, sinRotHelical, cosRotHelical);
1617  double d100 = interpolatedElement3DHelicalInt(Vin, x0, y0, z1, zHelical, sinRotHelical, cosRotHelical);
1618  double d101 = interpolatedElement3DHelicalInt(Vin, x1, y0, z1, zHelical, sinRotHelical, cosRotHelical);
1619  double d110 = interpolatedElement3DHelicalInt(Vin, x0, y1, z1, zHelical, sinRotHelical, cosRotHelical);
1620  double d111 = interpolatedElement3DHelicalInt(Vin, x1, y1, z1, zHelical, sinRotHelical, cosRotHelical);
1621 
1622  double dx00 = LIN_INTERP(fx, d000, d001);
1623  double dx01 = LIN_INTERP(fx, d100, d101);
1624  double dx10 = LIN_INTERP(fx, d010, d011);
1625  double dx11 = LIN_INTERP(fx, d110, d111);
1626  double dxy0 = LIN_INTERP(fy, dx00, dx10);
1627  double dxy1 = LIN_INTERP(fy, dx01, dx11);
1628 
1629  return LIN_INTERP(fz, dxy0, dxy1);
1630 }
1631 
1632 void symmetry_Helical(MultidimArray<double> &Vout, const MultidimArray<double> &Vin, double zHelical, double rotHelical,
1633  double rot0, MultidimArray<int> *mask, bool dihedral, double heightFraction, int Cn)
1634 {
1635  int zFirst=FIRST_XMIPP_INDEX(round(heightFraction*ZSIZE(Vin)));
1636  int zLast=LAST_XMIPP_INDEX(round(heightFraction*ZSIZE(Vin)));
1637 
1638  Vout.initZeros(Vin);
1639  double izHelical=1.0/zHelical;
1640  int zHelical2=(int)std::floor(0.5*zHelical);
1641  double sinRotHelical = sin(rotHelical);
1642  double cosRotHelical = cos(rotHelical);
1643 
1644  int Llength=ceil(ZSIZE(Vin)*izHelical);
1645 
1646  Matrix1D<double> sinCn, cosCn;
1647  sinCn.initZeros(Cn);
1648  cosCn.initZeros(Cn);
1649  for (int n=0; n<Cn; ++n)
1650  {
1651  VEC_ELEM(sinCn,n)=sin(n*TWOPI/Cn);
1652  VEC_ELEM(cosCn,n)=cos(n*TWOPI/Cn);
1653  }
1654 
1655  FOR_ALL_ELEMENTS_IN_ARRAY3D(Vin)
1656  {
1657  if (mask!=nullptr && !A3D_ELEM(*mask,k,i,j))
1658  continue;
1659  double rot=atan2((double)i,(double)j)+rot0;
1660  double rho=sqrt((double)i*i+(double)j*j);
1661  int l0=(int)ceil((STARTINGZ(Vin)-k)*izHelical);
1662  int lF=l0+Llength;
1663  double finalValue=0;
1664  double L=0;
1665  for (int il=l0; il<=lF; ++il)
1666  {
1667  double l=il;
1668  double kp=k+l*zHelical;
1669  if (kp>=zFirst && kp<=zLast)
1670  {
1671  double rotp=rot+l*rotHelical;
1672  double ip = rho*sin(rotp);
1673  double jp = rho*cos(rotp);
1674  double weight = 1.0;
1675  if (kp-zFirst<=zHelical2)
1676  weight=(kp-zFirst+1)/(zHelical2+1);
1677  else if (zLast-kp<=zHelical2)
1678  weight=(zLast+1-kp)/(zHelical2+1);
1679  finalValue+=weight*interpolatedElement3DHelical(Vin,jp,ip,kp,zHelical,sinRotHelical,cosRotHelical);
1680  L+=weight;
1681 
1682  for (int n=1; n<Cn; ++n)
1683  {
1684  double jpp=VEC_ELEM(cosCn,n)*jp-VEC_ELEM(sinCn,n)*ip;
1685  double ipp=VEC_ELEM(sinCn,n)*jp+VEC_ELEM(cosCn,n)*ip;
1686  finalValue+=weight*interpolatedElement3DHelical(Vin,jpp,ipp,kp,zHelical,sinRotHelical,cosRotHelical);
1687  L+=weight;
1688  }
1689  if (dihedral)
1690  {
1691  finalValue+=weight*interpolatedElement3DHelical(Vin,jp,-ip,-kp,zHelical,sinRotHelical,cosRotHelical);
1692  L+=weight;
1693  for (int n=1; n<Cn; ++n)
1694  {
1695  double jpp=VEC_ELEM(cosCn,n)*jp-VEC_ELEM(sinCn,n)*(-ip);
1696  double ipp=VEC_ELEM(sinCn,n)*jp+VEC_ELEM(cosCn,n)*(-ip);
1697  finalValue+=weight*interpolatedElement3DHelical(Vin,jpp,ipp,-kp,zHelical,sinRotHelical,cosRotHelical);
1698  L+=weight;
1699  }
1700  }
1701  }
1702  }
1703  A3D_ELEM(Vout,k,i,j)=finalValue/L;
1704  }
1705 }
1706 
1707 void symmetry_HelicalLowRes(MultidimArray<double> &Vout, const MultidimArray<double> &Vin, double zHelical, double rotHelical,
1708  double rot0, MultidimArray<int> *mask, double heightFraction)
1709 {
1710  MultidimArray<double> Vaux;
1711  Vout.initZeros(Vin);
1712  double helicalStep=rotHelical/zHelical;
1713  Matrix2D<double> A;
1714 
1715  for (int k=0; k<(int)ZSIZE(Vin); ++k)
1716  {
1717  double angle=RAD2DEG(helicalStep*k)+rot0;
1718  rotation3DMatrix(angle,'Z',A,true);
1719  MAT_ELEM(A,2,3)=-k;
1720  applyGeometry(xmipp_transformation::LINEAR,Vaux,Vin,A,xmipp_transformation::IS_NOT_INV,false,0.0);
1721  Vout+=Vaux;
1722 
1723  rotation3DMatrix(-angle,'Z',A,true);
1724  MAT_ELEM(A,2,3)=k;
1725  applyGeometry(xmipp_transformation::LINEAR,Vaux,Vin,A,xmipp_transformation::IS_NOT_INV,false,0.0);
1726  Vout+=Vaux;
1727  }
1728  Vout/=2*ZSIZE(Vin);
1729  if (mask!=nullptr)
1730  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vout)
1731  if (!DIRECT_MULTIDIM_ELEM(*mask,n))
1732  DIRECT_MULTIDIM_ELEM(Vout,n)=0.0;
1733 }
1734 
1735 void symmetry_Dihedral(MultidimArray<double> &Vout, const MultidimArray<double> &Vin,
1736  double rotStep, double zmin, double zmax, double zStep, MultidimArray<int> *mask)
1737 {
1738  // Find the best rotation
1739  MultidimArray<double> V180;
1740  Matrix2D<double> AZ, AX;
1741  rotation3DMatrix(180,'X',AX,true);
1742  applyGeometry(xmipp_transformation::LINEAR,V180,Vin,AX,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
1743  double bestCorr, bestRot, bestZ;
1744  bestCorr = bestRot = bestZ = std::numeric_limits<double>::min();
1745  double rot=-180.0;
1746  while (rot<180.0)
1747  {
1748  rotation3DMatrix(rot,'Z',AZ,true);
1749  double z=zmin;
1750  while (z<=zmax)
1751  {
1752  MAT_ELEM(AZ,2,3)=z;
1753  applyGeometry(xmipp_transformation::LINEAR,Vout,Vin,AZ,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
1754  double corr=correlationIndex(Vout,V180,mask);
1755  if (corr>bestCorr)
1756  {
1757  bestCorr=corr;
1758  bestRot=rot;
1759  bestZ=z;
1760  }
1761  z+=zStep;
1762  }
1763  rot+=rotStep;
1764  }
1765 
1766  rotation3DMatrix(-bestRot/2,'Z',AZ,true);
1767  MAT_ELEM(AZ,2,3)=-bestZ/2;
1768  applyGeometry(xmipp_transformation::BSPLINE3,V180,Vin,AZ*AX,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
1769  rotation3DMatrix(bestRot/2,'Z',AZ,true);
1770  MAT_ELEM(AZ,2,3)=bestZ/2;
1771  applyGeometry(xmipp_transformation::BSPLINE3,Vout,Vin,AZ,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
1772  Vout+=V180;
1773  Vout*=0.5;
1774 }
Symmetrize_Vol
#define Symmetrize_Vol(X)
Definition: symmetries.cpp:292
sym_P22_12_1
#define sym_P22_12_1
Definition: symmetries.h:62
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
VOLVOXEL
#define VOLVOXEL(V, k, i, j)
Definition: xmipp_image_base.h:214
SimpleGrid::highest
Matrix1D< double > highest
Definition: grids.h:161
ERR_NOT_IMPLEMENTED
Case or algorithm not implemented yet.
Definition: xmipp_error.h:177
sym_P22_12
#define sym_P22_12
Definition: symmetries.h:60
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
interpolatedElement3DHelical
double interpolatedElement3DHelical(const MultidimArray< double > &Vin, double x, double y, double z, double zHelical, double sinRotHelical, double cosRotHelical)
Definition: symmetries.cpp:1598
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
rho
double rho
Definition: numerical_recipes.cpp:6214
ipp
ipp
Definition: numerical_recipes.cpp:4714
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
y
static double * y
Definition: numerical_recipes.cpp:8487
xmipp_transformation::LINEAR
Definition: transformations_defines.h:32
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
sym_P4
#define sym_P4
Definition: symmetries.h:63
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
sym_P6
#define sym_P6
Definition: symmetries.h:68
TWOPI
#define TWOPI
Definition: xmipp_macros.h:111
symmetry_Dihedral
void symmetry_Dihedral(MultidimArray< double > &Vout, const MultidimArray< double > &Vin, double rotStep, double zmin, double zmax, double zStep, MultidimArray< int > *mask)
Definition: symmetries.cpp:1735
z0
#define z0
sym_P622
#define sym_P622
Definition: symmetries.h:69
FINISHINGZ
#define FINISHINGZ(v)
Definition: multidim_array_base.h:71
CC
#define CC
CC identifier.
Definition: grids.h:585
symmetry_HelicalLowRes
void symmetry_HelicalLowRes(MultidimArray< double > &Vout, const MultidimArray< double > &Vin, double zHelical, double rotHelical, double rot0, MultidimArray< int > *mask, double heightFraction)
Definition: symmetries.cpp:1707
symmetry_P2_122
void symmetry_P2_122(Image< double > &vol, const SimpleGrid &grid, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint, const MultidimArray< int > &mask, int volume_no, int grid_type)
Definition: symmetries.cpp:366
SimpleGrid::lowest
Matrix1D< double > lowest
Definition: grids.h:158
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
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
sym_P2_1
#define sym_P2_1
Definition: symmetries.h:56
rotation3DMatrix
void rotation3DMatrix(double ang, char axis, Matrix2D< double > &result, bool homogeneous)
Definition: transformations.cpp:426
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
symmetry_P42_12
void symmetry_P42_12(Image< double > &vol, const SimpleGrid &grid, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint, const MultidimArray< int > &mask, int volume_no, int grid_type)
Definition: symmetries.cpp:1009
BCC
#define BCC
BCC identifier.
Definition: grids.h:589
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
sym_P2
#define sym_P2
Definition: symmetries.h:55
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
sym_P2_122
#define sym_P2_122
Definition: symmetries.h:59
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
y0
#define y0
x0
#define x0
XX
#define XX(v)
Definition: matrix1d.h:85
sym_P312
#define sym_P312
Definition: symmetries.h:67
symmetry_P4
void symmetry_P4(Image< double > &vol, const SimpleGrid &grid, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint, const MultidimArray< int > &mask, int volume_no, int grid_type)
Definition: symmetries.cpp:815
symmetry_P6
void symmetry_P6(Image< double > &vol, const SimpleGrid &grid, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint, const MultidimArray< int > &mask, int volume_no, int grid_type)
Definition: symmetries.cpp:1331
FCC
#define FCC
FCC identifier.
Definition: grids.h:587
symmetry_Helical
void symmetry_Helical(MultidimArray< double > &Vout, const MultidimArray< double > &Vin, double zHelical, double rotHelical, double rot0, MultidimArray< int > *mask, bool dihedral, double heightFraction, int Cn)
Definition: symmetries.cpp:1632
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
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
interpolatedElement3DHelicalInt
double interpolatedElement3DHelicalInt(const MultidimArray< double > &Vin, int x, int y, int z, double zHelical, double sinRotHelical, double cosRotHelical)
Definition: symmetries.cpp:1577
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
sym_P42_12
#define sym_P42_12
Definition: symmetries.h:65
symmetry_P22_12
void symmetry_P22_12(Image< double > &vol, const SimpleGrid &grid, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint, const MultidimArray< int > &mask, int volume_no, int grid_type)
Definition: symmetries.cpp:591
j
#define j
Definition: numerical_recipes.cpp:2493
YY
#define YY(v)
Definition: matrix1d.h:93
sym_C2
#define sym_C2
Definition: symmetries.h:57
symmetrizeCrystalVectors
void symmetrizeCrystalVectors(Matrix1D< double > &aint, Matrix1D< double > &bint, Matrix1D< double > &shift, int space_group, int sym_no, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint)
Definition: symmetries.cpp:44
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
sym_P3
#define sym_P3
Definition: symmetries.h:66
alglib::round
int round(double x)
Definition: ap.cpp:7245
put_inside
#define put_inside(j, j_min, j_max, jint)
Definition: symmetries.cpp:360
sym_P422
#define sym_P422
Definition: symmetries.h:64
MultidimArrayBase::outside
bool outside(int k, int i, int j) const
Definition: multidim_array_base.h:956
RAD2DEG
#define RAD2DEG(r)
Definition: xmipp_macros.h:320
FIRST_XMIPP_INDEX
#define FIRST_XMIPP_INDEX(size)
Definition: xmipp_macros.h:439
sym_P1
#define sym_P1
Definition: symmetries.h:54
VECTOR_R3
#define VECTOR_R3(v, x, y, z)
Definition: matrix1d.h:124
sym_P222
#define sym_P222
Definition: symmetries.h:58
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
LAST_XMIPP_INDEX
#define LAST_XMIPP_INDEX(size)
Definition: xmipp_macros.h:448
LIN_INTERP
#define LIN_INTERP(a, l, h)
Definition: xmipp_macros.h:381
STARTINGZ
#define STARTINGZ(v)
Definition: multidim_array_base.h:67
n
int * n
Definition: numerical_recipes.cpp:2229
sym_undefined
#define sym_undefined
Definition: symmetries.h:53
symmetrizeCrystalVolume
void symmetrizeCrystalVolume(GridVolume &vol_in, const Matrix1D< double > &eprm_aint, const Matrix1D< double > &eprm_bint, int eprm_space_group, const MultidimArray< int > &mask, int grid_type)
Definition: symmetries.cpp:301
ZZ
#define ZZ(v)
Definition: matrix1d.h:101
xmipp_transformation::BSPLINE3
Definition: transformations_defines.h:32
Generated by   doxygen 1.8.13