Xmipp  v3.23.11-Nereus
angular_resolution_alignment.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Jose Luis Vilas (jlvilas@cnb.csic.es)
4  *
5  * Yale University, New Haven, Connecticut, United States of America
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 "angular_resolution_alignment.h"
27 #include <core/metadata_extension.h>
28 #include <data/monogenic_signal.h>
29 #include <data/fourier_filter.h>
30 #include <random>
31 #include <limits>
32 #include <CTPL/ctpl_stl.h>
33 
34 void ProgAngResAlign::defineParams()
35 {
36  addUsageLine("This algorithm estimates the existence of angular alignment errors by measuring a set of directional FSC.");
37  addUsageLine("To do that, the method only makes use of two half maps, but it does not use the set of particles.");
38  addUsageLine("The result is a plot resolution-radius. When this plot presents a slope the map has angular assignment error.");
39  addUsageLine("In contrast, if the curve is flat, the map is error free or there is shift errors in the particle alignment.");
40  addUsageLine("Reference: J.L. Vilas, et. al, XXXXX (2023)");
41  addUsageLine("+ The optimal con angle is 17 degree.", true);
42  addUsageLine(" This fact was proved in J.L Vilas & H.D. Tagare Nat. Meth 2023. Other values can be used,");
43  addUsageLine(" and this parameter does not seems to affect in a significative manner");
44  addUsageLine("+* On the helix", true);
45  addUsageLine(" If the map is a helix, the helix option should be activated. This option ensures a better estimation of the angular");
46  addUsageLine(" assignment errors. The helix is assumen that is oriented along the z-axis. This flag modifies the shape of the ");
47  addUsageLine("gaussian mask usign a cylinder.");
48  addUsageLine(" ");
49  addUsageLine(" ");
50  addSeeAlsoLine("resolution_fsc");
51 
52  addParamsLine(" --half1 <input_file> : Input Half map 1");
53  addParamsLine(" --half2 <input_file> : Input Half map 2");
54  addParamsLine(" [--directional_resolution] : (Optional) Uses direcitonal FSC instead of global FSC shell ");
55  addParamsLine(" [--limit_radius] : (Optional) Limits the maximum radius ");
56 
57  addParamsLine(" [-o <output_folder=\"\">] : Folder where the results will be stored.");
58 
59  addParamsLine(" [--sampling <Ts=1>] : (Optical) Pixel size (Angstrom). If it is not provided by default will be 1 A/px.");
60  addParamsLine(" [--mask <input_file=\"\">] : (Optional) Smooth mask to remove noise. If it is not provided, the computation will be carried out without mask.");
61 
62  addParamsLine(" [--anglecone <ang_con=17>] : (Optional) Angle Cone (angle between the axis and the generatrix) for estimating the directional FSC");
63  addParamsLine(" [--threshold <thrs=0.143>] : (Optional) Threshold for the FSC/directionalFSC estimation ");
64  addParamsLine(" [--helix] : (Optional) If the reconstruction is a helix put this flag. The axis of the helix must be along the z-axis");
65  addParamsLine(" [--threads <Nthreads=1>] : (Optional) Number of threads to be used");
66 
67  addExampleLine("Resolution of two half maps half1.mrc and half2.mrc with a sampling rate of 2 A/px", false);
68  addExampleLine("xmipp_angular_resolution_alignment --half1 half1.mrc --half2 half2.mrc --sampling 2 ");
69  addExampleLine("Resolution of two half maps half1.mrc and half2.mrc with a sampling rate of 2 A/px and a mask mask.mrc", false);
70  addExampleLine("xmipp_angular_resolution_alignment --half1 half1.mrc --half2 half2.mrc --mask mask.mrc --sampling 2 ");
71 }
72 
73 void ProgAngResAlign::readParams()
74 {
75  fnhalf1 = getParam("--half1");
76  fnhalf2 = getParam("--half2");
77  fnOut = getParam("-o");
78 
79  sampling = getDoubleParam("--sampling");
80  fnmask = getParam("--mask");
81  directionalRes = checkParam("--directional_resolution");
82  limRad = checkParam("--limit_radius");
83  isHelix = checkParam("--helix");
84  ang_con = getDoubleParam("--anglecone");
85  thrs = getDoubleParam("--threshold");
86 
87  Nthreads = getIntParam("--threads");
88 }
89 
90 void ProgAngResAlign::defineFrequenciesSimple(const MultidimArray<double> &inputVol)
91 {
92 
93  xvoldim = XSIZE(inputVol);
94  yvoldim = YSIZE(inputVol);
95  zvoldim = ZSIZE(inputVol);
96 
97  // Setting the Fouer dimensions
98  size_t xdimFourier = XSIZE(inputVol)/2 + 1;
99  size_t zdimFourier = zvoldim;
100  size_t ydimFourier = yvoldim;
101 
102  // Initializing the frequency vectors
103  freq_fourier_z.initZeros(zdimFourier);
104  freq_fourier_x.initZeros(xdimFourier);
105  freq_fourier_y.initZeros(ydimFourier);
106 
107  // u is the frequency
108  double u;
109 
110  // Defining frequency components. First element should be 0, it is set as the smallest number to avoid singularities
111 
112  VEC_ELEM(freq_fourier_z,0) = std::numeric_limits<double>::min();
113  for(size_t k=1; k<zdimFourier; ++k){
114  FFT_IDX2DIGFREQ(k,zvoldim, u);
115  VEC_ELEM(freq_fourier_z, k) = u;
116  }
117 
118  VEC_ELEM(freq_fourier_y,0) = std::numeric_limits<double>::min();
119  for(size_t k=1; k<ydimFourier; ++k){
120  FFT_IDX2DIGFREQ(k,yvoldim, u);
121  VEC_ELEM(freq_fourier_y, k) = u;
122  }
123 
124  VEC_ELEM(freq_fourier_x,0) = std::numeric_limits<double>::min();
125  for(size_t k=1; k<xdimFourier; ++k){
126  FFT_IDX2DIGFREQ(k,xvoldim, u);
127  VEC_ELEM(freq_fourier_x, k) = u;
128  }
129 
130  //Initializing map with frequencies
131 
132  freqMap.initZeros(1, zdimFourier, ydimFourier, xdimFourier); //Nyquist is 2, we take 1.9 greater than Nyquist
133 
134  freqElems.initZeros(xvoldim/2+1);
135 
136  // Directional frequencies along each direction
137  double uz, uy, ux, uz2, uz2y2;
138  long n=0;
139  int idx = 0;
140 
141  // Ncomps is the number of frequencies lesser than Nyquist
142  long Ncomps = 0;
143 
144  for(size_t k=0; k<zdimFourier; ++k)
145  {
146  uz = VEC_ELEM(freq_fourier_z, k);
147  uz2 = uz*uz;
148  for(size_t i=0; i<ydimFourier; ++i)
149  {
150  uy = VEC_ELEM(freq_fourier_y, i);
151  uz2y2 = uz2 + uy*uy;
152 
153  for(size_t j=0; j<xdimFourier; ++j)
154  {
155  ux = VEC_ELEM(freq_fourier_x, j);
156  ux = sqrt(uz2y2 + ux*ux);
157 
158  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1.9;
159  if (ux<=0.5)
160  {
161  idx = (int) round(ux * xvoldim);
162  ++Ncomps;
163  DIRECT_MULTIDIM_ELEM(freqElems, idx) += 1;
164 
165  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1/ux;
166 
167  if ((j == 0) && (uy<0))
168  {
169  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1.9;
170  DIRECT_MULTIDIM_ELEM(freqElems,idx) -= 1;
171  --Ncomps;
172  }
173 
174  if ((i == 0) && (j == 0) && (uz<0))
175  {
176  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1.9;
177  DIRECT_MULTIDIM_ELEM(freqElems,idx) -= 1;
178  --Ncomps;
179  }
180 
181  }
182  ++n;
183  }
184  }
185  }
186  real_z1z2.initZeros(Ncomps);
187 
188 }
189 
190 
191 void ProgAngResAlign::arrangeFSC_and_fscGlobal()
192  {
193  // cumpos is the the cumulative number of frequencies per shell number
194  // First shell has 0 elements
195  // second shell has the number of elements of the first shell
196  // Third shell has the number of elements of the first+sencond shells and so on
197  // freqElems in .h and set in defineFreq function
198  cumpos.resizeNoCopy(NZYXSIZE(freqElems));
199 
200  DIRECT_MULTIDIM_ELEM(cumpos,0) = 0;
201  for (long n = 1; n<NZYXSIZE(cumpos); ++n)
202  {
203  DIRECT_MULTIDIM_ELEM(cumpos,n) = DIRECT_MULTIDIM_ELEM(cumpos,n-1) + DIRECT_MULTIDIM_ELEM(freqElems,n-1);
204  }
205 
206  // real_z1z2, absz1_vec and absz2_vec will storage real(z1*conj(z2)), ||z1|| and ||z2|| using the
207  //Fourier Coefficients z1, and z2 of both half maps. The storage will be in order, first elements of
208  // the vector are the freq 0, next freq 1 and so on. Ncomps is the number of components with frequency
209  // lesser than Nyquist (defined in defineFreq)
210 
211  absz1_vec = real_z1z2;
212  absz2_vec = real_z1z2;
213 
214  // num and den are de numerator and denominator of the fsc
215  MultidimArray<long> pos;
216  MultidimArray<double> num, den1, den2, noiseAmp;
217  num.initZeros(freqElems);
218  pos.resizeNoCopy(num);
219 
220  freqidx.resizeNoCopy(real_z1z2);
221 
222  // arr2indx is the position of in Fourier space. It is defined in the .h
223  // Because we only work with the frequencies lesser than nyquist their position
224  // int he Fourier space is lost. arr2indx has the position of each frequency.
225  // It is used at the end of the algorithm to generate the 3DFSC
226  freqidx.initZeros();
227  pos.initZeros();
228  den1 = num;
229  den2 = num;
230 
231  auto ZdimFT1=(int)ZSIZE(FT1);
232  auto YdimFT1=(int)YSIZE(FT1);
233  auto XdimFT1=(int)XSIZE(FT1);
234 
235  // fx, fy, fz, defined in the .h are the frequencies of each pixel with frequency
236  // lesser than Nyquist along each axis. They are used to compute the dot product
237  // between each vector position and the direction of the cone to calculate the
238  // directional FSC
239  fx.resizeNoCopy(real_z1z2);
240  fy.resizeNoCopy(real_z1z2);
241  fz.resizeNoCopy(real_z1z2);
242 
243  long n = 0;
244  for (int k=0; k<ZdimFT1; k++)
245  {
246  double uz = VEC_ELEM(freq_fourier_z, k);
247  for (int i=0; i<YdimFT1; i++)
248  {
249  double uy = VEC_ELEM(freq_fourier_y, i);
250  for (int j=0; j<XdimFT1; j++)
251  {
252  double ux = VEC_ELEM(freq_fourier_x, j);
253 
254  double iun = DIRECT_MULTIDIM_ELEM(freqMap,n);
255  double f = 1/iun;
256  ++n;
257 
258  // Only reachable frequencies
259  // To speed up the algorithm, only are considered those voxels with frequency lesser than Nyquist, 0.5. The vector idx_count
260  // stores all frequencies lesser than Nyquist. This vector determines the frequency of each component of
261  // real_z1z2, absz1_vec, absz2_vec.
262  if (f>0.5)
263  continue;
264 
265  // Index of each frequency
266  auto idx = (int) round(f * xvoldim);
267 
268  int idx_count = DIRECT_MULTIDIM_ELEM(cumpos, idx) + DIRECT_MULTIDIM_ELEM(pos, idx);
269 
270  // Storing normalized frequencies
271  DIRECT_MULTIDIM_ELEM(fx, idx_count) = (float) ux*iun;
272  DIRECT_MULTIDIM_ELEM(fy, idx_count) = (float) uy*iun;
273  DIRECT_MULTIDIM_ELEM(fz, idx_count) = (float) uz*iun;
274 
275  // In this vector we store the index of each Fourier coefficient, thus in the
276  // directional FSC estimation the calculus of the index is avoided speeding up the calculation
277  DIRECT_MULTIDIM_ELEM(freqidx, idx_count) = idx;
278 
279  // Fourier coefficients of both halves
280  std::complex<double> &z1 = dAkij(FT1, k, i, j);
281  std::complex<double> &z2 = dAkij(FT2, k, i, j);
282 
283  double absz1 = abs(z1);
284  double absz2 = abs(z2);
285 
286  DIRECT_MULTIDIM_ELEM(real_z1z2, idx_count) = (float) real(conj(z1)*z2);
287  DIRECT_MULTIDIM_ELEM(absz1_vec, idx_count) = (float) absz1*absz1;
288  DIRECT_MULTIDIM_ELEM(absz2_vec, idx_count) = (float) absz2*absz2;
289 
290  DIRECT_MULTIDIM_ELEM(pos, idx) +=1;
291  }
292  }
293  }
294  }
295 
296 void ProgAngResAlign::fscGlobal(double &threshold, double &resol)
297  {
298  // num and den are de numerator and denominator of the fsc
299  MultidimArray<double> num, den1, den2;
300  num.initZeros(freqElems);
301  den1 = num;
302  den2 = num;
303 
304  auto ZdimFT1=(int)ZSIZE(FT1);
305  auto YdimFT1=(int)YSIZE(FT1);
306  auto XdimFT1=(int)XSIZE(FT1);
307 
308  long n = 0;
309  for (int k=0; k<ZdimFT1; k++)
310  {
311  for (int i=0; i<YdimFT1; i++)
312  {
313  for (int j=0; j<XdimFT1; j++)
314  {
315  auto iun = DIRECT_MULTIDIM_ELEM(freqMap,n);
316  double f = 1/iun;
317  ++n;
318 
319  // Only reachable frequencies
320  // To speed up the algorithm, only are considered those voxels with frequency lesser than Nyquist, 0.5. The vector idx_count
321  // stores all frequencies lesser than Nyquist. This vector determines the frequency of each component of
322  // real_z1z2, absz1_vec, absz2_vec.
323  if (f>0.5)
324  continue;
325 
326  // Index of each frequency
327  auto idx = (int) round(f * xvoldim);
328 
329  // Fourier coefficients of both halves
330  std::complex<double> &z1 = dAkij(FT1, k, i, j);
331  std::complex<double> &z2 = dAkij(FT2, k, i, j);
332 
333  double absz1 = abs(z1);
334  double absz2 = abs(z2);
335 
336  DIRECT_MULTIDIM_ELEM(num, idx) = (float) real(conj(z1)*z2);
337  DIRECT_MULTIDIM_ELEM(den1, idx) = (float) absz1*absz1;
338  DIRECT_MULTIDIM_ELEM(den2, idx) = (float) absz2*absz2;
339  }
340  }
341  }
342  MultidimArray<double> fsc;
343  fsc.resizeNoCopy(num);
344  fsc.initConstant(1.0);
345 
346  //The fsc is stored in a metadata and saved
347  bool flagRes = true;
348  FOR_ALL_ELEMENTS_IN_ARRAY1D(num)
349  {
350  double auxfsc = (dAi(num,i))/(sqrt(dAi(den1,i)*dAi(den2,i))+1e-38);
351  dAi(fsc,i) = std::max(0.0, auxfsc);
352 
353  if (flagRes && (i>2) && (dAi(fsc,i)<=threshold))
354  {
355  flagRes = false;
356  double ff = (double) i / (xvoldim * sampling); // frequency
357  resol = 1./ff;
358  }
359  }
360 
361 
362  }
363 
364 
365 void ProgAngResAlign::fscInterpolation(const MultidimArray<double> &freq, const MultidimArray< double > &frc)
366 {
367  // Here the FSC at 0.143 is obtained by interpolating
368  FOR_ALL_ELEMENTS_IN_ARRAY1D(freq)
369  {
370  auto ff = dAi(frc,i);
371  if ( (ff<=thrs) && (i>2) )
372  {
373  double y2, y1, x2, x1, slope, ny;
374  y1 = ff;
375  x1 = dAi(freq,i);
376  y2 = dAi(frc, i-1);
377  x2 = dAi(freq, i-1);
378  slope = (y2 - y1)/(x2 - x1);
379  ny = y2 - slope*x2;
380 
381  double fscResolution;
382  fscResolution = (thrs - ny)/slope;
383  std::cout << "Resolution " << 1/fscResolution << std::endl;
384  break;
385  }
386  }
387 }
388 
389 
390 
391 void ProgAngResAlign::fscDir_fast(MultidimArray<float> &fsc, double rot, double tilt,
392  double &thrs, double &resol)
393 {
394  // FSCDIR_FAST: computes the directional FSC along a direction given by the angles rot and tilt.
395  // Thus the directional resolution, resol, is estimated with the threshold, thrs.
396  // The direcional FSC is stored in a metadata mdRes and in the multidimarray fsc. Later, this multidimarray
397  // will be used to estimate the FSO.
398  // In addition, the 3dfsc and a normalizationmap (needed to estimate the 3dfsc) are created. For each direction
399  // these vectors are updated
400  size_t dim = NZYXSIZE(freqElems);
401 
402  // numerator and denominator of the fsc
403  MultidimArray<float> num, den1, den2;
404  num.initZeros(dim);
405  den1.initZeros(dim);
406  den2.initZeros(dim);
407 
408  //Parameter to determine the direction of the cone
409  float x_dir, y_dir, z_dir, cosAngle, aux;
410  x_dir = sinf(tilt)*cosf(rot);
411  y_dir = sinf(tilt)*sinf(rot);
412  z_dir = cosf(tilt);
413 
414  cosAngle = (float) cos(ang_con);
415 
416  // It is multiply by 0.5 because later the weight is
417  // cosine = sqrt(exp( -((cosine -1)*(cosine -1))*aux ));
418  // thus the computation of the weight is speeded up
419  // aux = 4.0/((cos(ang_con) -1)*(cos(ang_con) -1));
420  aux = (4.0/((cosAngle -1)*(cosAngle -1)));//*0.5;
421 
422  // Computing directional resolution
423  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(real_z1z2)
424  {
425  // frecuencies of the fourier position
426  auto ux = DIRECT_MULTIDIM_ELEM(fx, n);
427  auto uy = DIRECT_MULTIDIM_ELEM(fy, n);
428  auto uz = DIRECT_MULTIDIM_ELEM(fz, n);
429 
430  // angle between the position and the direction of the cone
431  float cosine = fabs(x_dir*ux + y_dir*uy + z_dir*uz);
432 
433  // Only points inside the cone
434  if (cosine >= cosAngle)
435  {
436  // in fact it is sqrt(exp( -((cosine -1)*(cosine -1))*aux ));
437  // For sake of performance The sqrt is avoided dividing the by 2 in aux
438  cosine = expf( -((cosine -1)*(cosine -1))*aux);
439 
440  // selecting the frequency of the shell
441  size_t idxf = DIRECT_MULTIDIM_ELEM(freqidx, n);
442 
443  // computing the numerator and denominator
444  dAi(num, idxf) += DIRECT_MULTIDIM_ELEM(real_z1z2, n) * cosine;
445  dAi(den1,idxf) += DIRECT_MULTIDIM_ELEM(absz1_vec, n) * cosine;
446  dAi(den2,idxf) += DIRECT_MULTIDIM_ELEM(absz2_vec, n) * cosine;
447  }
448  }
449 
450  fsc.resizeNoCopy(num);
451  fsc.initConstant(1.0);
452 
453  //The fsc is stored in a metadata and saved
454  bool flagRes = true;
455  FOR_ALL_ELEMENTS_IN_ARRAY1D(num)
456  {
457  double auxfsc = (dAi(num,i))/(sqrt(dAi(den1,i)*dAi(den2,i))+1e-38);
458  dAi(fsc,i) = std::max(0.0, auxfsc);
459 
460  if (flagRes && (i>2) && (dAi(fsc,i)<=thrs))
461  {
462  flagRes = false;
463  double ff = (float) i / (xvoldim * sampling); // frequency
464  resol = 1./ff;
465  }
466  }
467 }
468 
469 
470 void ProgAngResAlign::generateDirections(Matrix2D<float> &angles, bool alot)
471 {
472  if (alot == true)
473  {
474  angles.initZeros(2,321);
475  MAT_ELEM(angles,0,0)=0.0f; MAT_ELEM(angles,1,0)=0.0f;
476  MAT_ELEM(angles,0,1)=324.0f; MAT_ELEM(angles,1,1)=63.4349f;
477  MAT_ELEM(angles,0,2)=36.0f; MAT_ELEM(angles,1,2)=63.4349f;
478  MAT_ELEM(angles,0,3)=180.0f; MAT_ELEM(angles,1,3)=63.435f;
479  MAT_ELEM(angles,0,4)=252.0f; MAT_ELEM(angles,1,4)=63.435f;
480  MAT_ELEM(angles,0,5)=108.0f; MAT_ELEM(angles,1,5)=63.435f;
481  MAT_ELEM(angles,0,6)=324.0f; MAT_ELEM(angles,1,6)=31.7175f;
482  MAT_ELEM(angles,0,7)=36.0f; MAT_ELEM(angles,1,7)=31.7175f;
483  MAT_ELEM(angles,0,8)=0.0f; MAT_ELEM(angles,1,8)=58.2825f;
484  MAT_ELEM(angles,0,9)=288.0f; MAT_ELEM(angles,1,9)=58.2825f;
485  MAT_ELEM(angles,0,10)=342.0f; MAT_ELEM(angles,1,10)=90.0f;
486  MAT_ELEM(angles,0,11)=306.0f; MAT_ELEM(angles,1,11)=90.0f;
487  MAT_ELEM(angles,0,12)=72.0f; MAT_ELEM(angles,1,12)=58.2825f;
488  MAT_ELEM(angles,0,13)=18.0f; MAT_ELEM(angles,1,13)=90.0f;
489  MAT_ELEM(angles,0,14)=54.0f; MAT_ELEM(angles,1,14)=90.0f;
490  MAT_ELEM(angles,0,15)=90.0f; MAT_ELEM(angles,1,15)=90.0f;
491  MAT_ELEM(angles,0,16)=216.0f; MAT_ELEM(angles,1,16)=58.282f;
492  MAT_ELEM(angles,0,17)=144.0f; MAT_ELEM(angles,1,17)=58.282f;
493  MAT_ELEM(angles,0,18)=180.0f; MAT_ELEM(angles,1,18)=31.718f;
494  MAT_ELEM(angles,0,19)=252.0f; MAT_ELEM(angles,1,19)=31.718f;
495  MAT_ELEM(angles,0,20)=108.0f; MAT_ELEM(angles,1,20)=31.718f;
496  MAT_ELEM(angles,0,21)=346.3862f; MAT_ELEM(angles,1,21)=43.6469f;
497  MAT_ELEM(angles,0,22)=58.3862f; MAT_ELEM(angles,1,22)=43.6469f;
498  MAT_ELEM(angles,0,23)=274.3862f; MAT_ELEM(angles,1,23)=43.6469f;
499  MAT_ELEM(angles,0,24)=0.0f; MAT_ELEM(angles,1,24)=90.0f;
500  MAT_ELEM(angles,0,25)=72.0f; MAT_ELEM(angles,1,25)=90.0f;
501  MAT_ELEM(angles,0,26)=288.0f; MAT_ELEM(angles,1,26)=90.0f;
502  MAT_ELEM(angles,0,27)=225.7323f; MAT_ELEM(angles,1,27)=73.955f;
503  MAT_ELEM(angles,0,28)=153.7323f; MAT_ELEM(angles,1,28)=73.955f;
504  MAT_ELEM(angles,0,29)=216.0f; MAT_ELEM(angles,1,29)=26.565f;
505  MAT_ELEM(angles,0,30)=144.0f; MAT_ELEM(angles,1,30)=26.565f;
506  MAT_ELEM(angles,0,31)=0.0f; MAT_ELEM(angles,1,31)=26.5651f;
507  MAT_ELEM(angles,0,32)=72.0f; MAT_ELEM(angles,1,32)=26.5651f;
508  MAT_ELEM(angles,0,33)=288.0f; MAT_ELEM(angles,1,33)=26.5651f;
509  MAT_ELEM(angles,0,34)=350.2677f; MAT_ELEM(angles,1,34)=73.9549f;
510  MAT_ELEM(angles,0,35)=62.2677f; MAT_ELEM(angles,1,35)=73.9549f;
511  MAT_ELEM(angles,0,36)=278.2677f; MAT_ELEM(angles,1,36)=73.9549f;
512  MAT_ELEM(angles,0,37)=206.2677f; MAT_ELEM(angles,1,37)=73.955f;
513  MAT_ELEM(angles,0,38)=134.2677f; MAT_ELEM(angles,1,38)=73.955f;
514  MAT_ELEM(angles,0,39)=202.3862f; MAT_ELEM(angles,1,39)=43.647f;
515  MAT_ELEM(angles,0,40)=130.3862f; MAT_ELEM(angles,1,40)=43.647f;
516  MAT_ELEM(angles,0,41)=13.6138f; MAT_ELEM(angles,1,41)=43.6469f;
517  MAT_ELEM(angles,0,42)=85.6138f; MAT_ELEM(angles,1,42)=43.6469f;
518  MAT_ELEM(angles,0,43)=301.6138f; MAT_ELEM(angles,1,43)=43.6469f;
519  MAT_ELEM(angles,0,44)=9.7323f; MAT_ELEM(angles,1,44)=73.9549f;
520  MAT_ELEM(angles,0,45)=81.7323f; MAT_ELEM(angles,1,45)=73.9549f;
521  MAT_ELEM(angles,0,46)=297.7323f; MAT_ELEM(angles,1,46)=73.9549f;
522  MAT_ELEM(angles,0,47)=36.0f; MAT_ELEM(angles,1,47)=90.0f;
523  MAT_ELEM(angles,0,48)=324.0f; MAT_ELEM(angles,1,48)=90.0f;
524  MAT_ELEM(angles,0,49)=229.6138f; MAT_ELEM(angles,1,49)=43.647f;
525  MAT_ELEM(angles,0,50)=157.6138f; MAT_ELEM(angles,1,50)=43.647f;
526  MAT_ELEM(angles,0,51)=324.0f; MAT_ELEM(angles,1,51)=15.8587f;
527  MAT_ELEM(angles,0,52)=36.0f; MAT_ELEM(angles,1,52)=15.8587f;
528  MAT_ELEM(angles,0,53)=341.533f; MAT_ELEM(angles,1,53)=59.6208f;
529  MAT_ELEM(angles,0,54)=306.467f; MAT_ELEM(angles,1,54)=59.6208f;
530  MAT_ELEM(angles,0,55)=333.5057f; MAT_ELEM(angles,1,55)=76.5584f;
531  MAT_ELEM(angles,0,56)=314.4943f; MAT_ELEM(angles,1,56)=76.5584f;
532  MAT_ELEM(angles,0,57)=53.533f; MAT_ELEM(angles,1,57)=59.6208f;
533  MAT_ELEM(angles,0,58)=26.4943f; MAT_ELEM(angles,1,58)=76.5584f;
534  MAT_ELEM(angles,0,59)=45.5057f; MAT_ELEM(angles,1,59)=76.5584f;
535  MAT_ELEM(angles,0,60)=197.533f; MAT_ELEM(angles,1,60)=59.621f;
536  MAT_ELEM(angles,0,61)=162.467f; MAT_ELEM(angles,1,61)=59.621f;
537  MAT_ELEM(angles,0,62)=180.0; MAT_ELEM(angles,1,62)=47.576;
538  MAT_ELEM(angles,0,63)=269.533f; MAT_ELEM(angles,1,63)=59.621;
539  MAT_ELEM(angles,0,64)=252.0f; MAT_ELEM(angles,1,64)=47.576f;
540  MAT_ELEM(angles,0,65)=108.0f; MAT_ELEM(angles,1,65)=47.576f;
541  MAT_ELEM(angles,0,66)=324.0f; MAT_ELEM(angles,1,66)=47.5762f;
542  MAT_ELEM(angles,0,67)=36.0f; MAT_ELEM(angles,1,67)=47.5762f;
543  MAT_ELEM(angles,0,68)=18.467f; MAT_ELEM(angles,1,68)=59.6208f;
544  MAT_ELEM(angles,0,69)=170.4943f; MAT_ELEM(angles,1,69)=76.558f;
545  MAT_ELEM(angles,0,70)=117.5057f; MAT_ELEM(angles,1,70)=76.558f;
546  MAT_ELEM(angles,0,71)=189.5057f; MAT_ELEM(angles,1,71)=76.558f;
547  MAT_ELEM(angles,0,72)=242.4943f; MAT_ELEM(angles,1,72)=76.558f;
548  MAT_ELEM(angles,0,73)=261.5057f; MAT_ELEM(angles,1,73)=76.558f;
549  MAT_ELEM(angles,0,74)=98.4943f; MAT_ELEM(angles,1,74)=76.558f;
550  MAT_ELEM(angles,0,75)=234.467f; MAT_ELEM(angles,1,75)=59.621f;
551  MAT_ELEM(angles,0,76)=125.533f; MAT_ELEM(angles,1,76)=59.621f;
552  MAT_ELEM(angles,0,77)=180.0f; MAT_ELEM(angles,1,77)=15.859f;
553  MAT_ELEM(angles,0,78)=252.0f; MAT_ELEM(angles,1,78)=15.859f;
554  MAT_ELEM(angles,0,79)=90.467f; MAT_ELEM(angles,1,79)=59.621f;
555  MAT_ELEM(angles,0,80)=108.0f; MAT_ELEM(angles,1,80)=15.859f;
556  MAT_ELEM(angles,0,81)=0.0f; MAT_ELEM(angles,1,81)=42.8321f;
557  MAT_ELEM(angles,0,82)=72.0f; MAT_ELEM(angles,1,82)=42.8321f;
558  MAT_ELEM(angles,0,83)=288.0f; MAT_ELEM(angles,1,83)=42.8321f;
559  MAT_ELEM(angles,0,84)=4.7693f; MAT_ELEM(angles,1,84)=81.9488f;
560  MAT_ELEM(angles,0,85)=76.7693f; MAT_ELEM(angles,1,85)=81.9488f;
561  MAT_ELEM(angles,0,86)=292.7693f; MAT_ELEM(angles,1,86)=81.9488f;
562  MAT_ELEM(angles,0,87)=220.7693f; MAT_ELEM(angles,1,87)=81.9488f;
563  MAT_ELEM(angles,0,88)=148.7693f; MAT_ELEM(angles,1,88)=81.9488f;
564  MAT_ELEM(angles,0,89)=224.2677f; MAT_ELEM(angles,1,89)=34.924f;
565  MAT_ELEM(angles,0,90)=152.2677f; MAT_ELEM(angles,1,90)=34.924f;
566  MAT_ELEM(angles,0,91)=13.5146f; MAT_ELEM(angles,1,91)=20.3172f;
567  MAT_ELEM(angles,0,92)=85.5146f; MAT_ELEM(angles,1,92)=20.3172f;
568  MAT_ELEM(angles,0,93)=301.5146f; MAT_ELEM(angles,1,93)=20.3172f;
569  MAT_ELEM(angles,0,94)=346.1363f; MAT_ELEM(angles,1,94)=66.7276f;
570  MAT_ELEM(angles,0,95)=58.1363f; MAT_ELEM(angles,1,95)=66.7276f;
571  MAT_ELEM(angles,0,96)=274.1363f; MAT_ELEM(angles,1,96)=66.7276f;
572  MAT_ELEM(angles,0,97)=197.8362f; MAT_ELEM(angles,1,97)=75.105f;
573  MAT_ELEM(angles,0,98)=269.8362f; MAT_ELEM(angles,1,98)=75.105f;
574  MAT_ELEM(angles,0,99)=125.8362f; MAT_ELEM(angles,1,99)=75.105f;
575  MAT_ELEM(angles,0,100)=199.6899f; MAT_ELEM(angles,1,100)=51.609f;
576  MAT_ELEM(angles,0,101)=127.6899f; MAT_ELEM(angles,1,101)=51.609f;
577  MAT_ELEM(angles,0,102)=334.8124f; MAT_ELEM(angles,1,102)=45.0621f;
578  MAT_ELEM(angles,0,103)=46.8124f; MAT_ELEM(angles,1,103)=45.0621f;
579  MAT_ELEM(angles,0,104)=175.3133f; MAT_ELEM(angles,1,104)=83.2562f;
580  MAT_ELEM(angles,0,105)=247.3133f; MAT_ELEM(angles,1,105)=83.2562f;
581  MAT_ELEM(angles,0,106)=103.3133f; MAT_ELEM(angles,1,106)=83.2562f;
582  MAT_ELEM(angles,0,107)=229.8637f; MAT_ELEM(angles,1,107)=66.728f;
583  MAT_ELEM(angles,0,108)=157.8637f; MAT_ELEM(angles,1,108)=66.728f;
584  MAT_ELEM(angles,0,109)=202.4854f; MAT_ELEM(angles,1,109)=20.317f;
585  MAT_ELEM(angles,0,110)=130.4854f; MAT_ELEM(angles,1,110)=20.317f;
586  MAT_ELEM(angles,0,111)=16.3101f; MAT_ELEM(angles,1,111)=51.6091f;
587  MAT_ELEM(angles,0,112)=88.3101f; MAT_ELEM(angles,1,112)=51.6091f;
588  MAT_ELEM(angles,0,113)=304.3101f; MAT_ELEM(angles,1,113)=51.6091f;
589  MAT_ELEM(angles,0,114)=18.1638f; MAT_ELEM(angles,1,114)=75.1046f;
590  MAT_ELEM(angles,0,115)=306.1638f; MAT_ELEM(angles,1,115)=75.1046f;
591  MAT_ELEM(angles,0,116)=40.6867f; MAT_ELEM(angles,1,116)=83.2562f;
592  MAT_ELEM(angles,0,117)=328.6867f; MAT_ELEM(angles,1,117)=83.2562f;
593  MAT_ELEM(angles,0,118)=241.1876f; MAT_ELEM(angles,1,118)=45.062f;
594  MAT_ELEM(angles,0,119)=97.1876f; MAT_ELEM(angles,1,119)=45.062f;
595  MAT_ELEM(angles,0,120)=169.1876f; MAT_ELEM(angles,1,120)=45.062f;
596  MAT_ELEM(angles,0,121)=351.7323f; MAT_ELEM(angles,1,121)=34.9243f;
597  MAT_ELEM(angles,0,122)=63.7323f; MAT_ELEM(angles,1,122)=34.9243f;
598  MAT_ELEM(angles,0,123)=279.7323f; MAT_ELEM(angles,1,123)=34.9243f;
599  MAT_ELEM(angles,0,124)=355.2307f; MAT_ELEM(angles,1,124)=81.9488f;
600  MAT_ELEM(angles,0,125)=67.2307f; MAT_ELEM(angles,1,125)=81.9488f;
601  MAT_ELEM(angles,0,126)=283.2307f; MAT_ELEM(angles,1,126)=81.9488f;
602  MAT_ELEM(angles,0,127)=216.0f; MAT_ELEM(angles,1,127)=73.733f;
603  MAT_ELEM(angles,0,128)=144.0f; MAT_ELEM(angles,1,128)=73.733f;
604  MAT_ELEM(angles,0,129)=207.7323f; MAT_ELEM(angles,1,129)=34.924f;
605  MAT_ELEM(angles,0,130)=135.7323f; MAT_ELEM(angles,1,130)=34.924f;
606  MAT_ELEM(angles,0,131)=346.4854f; MAT_ELEM(angles,1,131)=20.3172f;
607  MAT_ELEM(angles,0,132)=58.4854f; MAT_ELEM(angles,1,132)=20.3172f;
608  MAT_ELEM(angles,0,133)=274.4854f; MAT_ELEM(angles,1,133)=20.3172f;
609  MAT_ELEM(angles,0,134)=341.8362f; MAT_ELEM(angles,1,134)=75.1046f;
610  MAT_ELEM(angles,0,135)=53.8362f; MAT_ELEM(angles,1,135)=75.1046f;
611  MAT_ELEM(angles,0,136)=202.1363f; MAT_ELEM(angles,1,136)=66.728f;
612  MAT_ELEM(angles,0,137)=130.1363f; MAT_ELEM(angles,1,137)=66.728f;
613  MAT_ELEM(angles,0,138)=190.8124f; MAT_ELEM(angles,1,138)=45.062f;
614  MAT_ELEM(angles,0,139)=262.8124f; MAT_ELEM(angles,1,139)=45.062f;
615  MAT_ELEM(angles,0,140)=118.8124f; MAT_ELEM(angles,1,140)=45.062f;
616  MAT_ELEM(angles,0,141)=343.6899f; MAT_ELEM(angles,1,141)=51.6091f;
617  MAT_ELEM(angles,0,142)=55.6899f; MAT_ELEM(angles,1,142)=51.6091f;
618  MAT_ELEM(angles,0,143)=271.6899f; MAT_ELEM(angles,1,143)=51.6091f;
619  MAT_ELEM(angles,0,144)=184.6867f; MAT_ELEM(angles,1,144)=83.2562f;
620  MAT_ELEM(angles,0,145)=256.6867f; MAT_ELEM(angles,1,145)=83.2562f;
621  MAT_ELEM(angles,0,146)=112.6867f; MAT_ELEM(angles,1,146)=83.2562f;
622  MAT_ELEM(angles,0,147)=234.1638f; MAT_ELEM(angles,1,147)=75.105f;
623  MAT_ELEM(angles,0,148)=90.1638f; MAT_ELEM(angles,1,148)=75.105f;
624  MAT_ELEM(angles,0,149)=162.1638f; MAT_ELEM(angles,1,149)=75.105f;
625  MAT_ELEM(angles,0,150)=229.5146f; MAT_ELEM(angles,1,150)=20.317f;
626  MAT_ELEM(angles,0,151)=157.5146f; MAT_ELEM(angles,1,151)=20.317f;
627  MAT_ELEM(angles,0,152)=25.1876f; MAT_ELEM(angles,1,152)=45.0621f;
628  MAT_ELEM(angles,0,153)=313.1876f; MAT_ELEM(angles,1,153)=45.0621f;
629  MAT_ELEM(angles,0,154)=13.8637f; MAT_ELEM(angles,1,154)=66.7276f;
630  MAT_ELEM(angles,0,155)=85.8637f; MAT_ELEM(angles,1,155)=66.7276f;
631  MAT_ELEM(angles,0,156)=301.8637f; MAT_ELEM(angles,1,156)=66.7276f;
632  MAT_ELEM(angles,0,157)=31.3133f; MAT_ELEM(angles,1,157)=83.2562f;
633  MAT_ELEM(angles,0,158)=319.3133f; MAT_ELEM(angles,1,158)=83.2562f;
634  MAT_ELEM(angles,0,159)=232.3101f; MAT_ELEM(angles,1,159)=51.609f;
635  MAT_ELEM(angles,0,160)=160.3101f; MAT_ELEM(angles,1,160)=51.609f;
636  MAT_ELEM(angles,0,161)=8.2677f; MAT_ELEM(angles,1,161)=34.9243f;
637  MAT_ELEM(angles,0,162)=80.2677f; MAT_ELEM(angles,1,162)=34.9243f;
638  MAT_ELEM(angles,0,163)=296.2677f; MAT_ELEM(angles,1,163)=34.9243f;
639  MAT_ELEM(angles,0,164)=0.0f; MAT_ELEM(angles,1,164)=73.733f;
640  MAT_ELEM(angles,0,165)=72.0f; MAT_ELEM(angles,1,165)=73.733f;
641  MAT_ELEM(angles,0,166)=288.0f; MAT_ELEM(angles,1,166)=73.733f;
642  MAT_ELEM(angles,0,167)=211.2307f; MAT_ELEM(angles,1,167)=81.9488f;
643  MAT_ELEM(angles,0,168)=139.2307f; MAT_ELEM(angles,1,168)=81.9488f;
644  MAT_ELEM(angles,0,169)=216.0f; MAT_ELEM(angles,1,169)=42.832f;
645  MAT_ELEM(angles,0,170)=144.0f; MAT_ELEM(angles,1,170)=42.832f;
646  MAT_ELEM(angles,0,171)=0.0f; MAT_ELEM(angles,1,171)=12.9432f;
647  MAT_ELEM(angles,0,172)=72.0f; MAT_ELEM(angles,1,172)=12.9432f;
648  MAT_ELEM(angles,0,173)=288.0f; MAT_ELEM(angles,1,173)=12.9432f;
649  MAT_ELEM(angles,0,174)=337.2786f; MAT_ELEM(angles,1,174)=68.041f;
650  MAT_ELEM(angles,0,175)=49.2786f; MAT_ELEM(angles,1,175)=68.041f;
651  MAT_ELEM(angles,0,176)=193.2786f; MAT_ELEM(angles,1,176)=68.041f;
652  MAT_ELEM(angles,0,177)=265.2786f; MAT_ELEM(angles,1,177)=68.041f;
653  MAT_ELEM(angles,0,178)=121.2786f; MAT_ELEM(angles,1,178)=68.041f;
654  MAT_ELEM(angles,0,179)=189.4537f; MAT_ELEM(angles,1,179)=53.278f;
655  MAT_ELEM(angles,0,180)=261.4537f; MAT_ELEM(angles,1,180)=53.278f;
656  MAT_ELEM(angles,0,181)=117.4537f; MAT_ELEM(angles,1,181)=53.278f;
657  MAT_ELEM(angles,0,182)=333.4537f; MAT_ELEM(angles,1,182)=53.2783f;
658  MAT_ELEM(angles,0,183)=45.4537f; MAT_ELEM(angles,1,183)=53.2783f;
659  MAT_ELEM(angles,0,184)=180.0f; MAT_ELEM(angles,1,184)=76.378f;
660  MAT_ELEM(angles,0,185)=252.0f; MAT_ELEM(angles,1,185)=76.378f;
661  MAT_ELEM(angles,0,186)=108.0f; MAT_ELEM(angles,1,186)=76.378f;
662  MAT_ELEM(angles,0,187)=238.7214f; MAT_ELEM(angles,1,187)=68.041f;
663  MAT_ELEM(angles,0,188)=94.7214f; MAT_ELEM(angles,1,188)=68.041f;
664  MAT_ELEM(angles,0,189)=166.7214f; MAT_ELEM(angles,1,189)=68.041f;
665  MAT_ELEM(angles,0,190)=216.0f; MAT_ELEM(angles,1,190)=12.943f;
666  MAT_ELEM(angles,0,191)=144.0f; MAT_ELEM(angles,1,191)=12.943f;
667  MAT_ELEM(angles,0,192)=26.5463f; MAT_ELEM(angles,1,192)=53.2783f;
668  MAT_ELEM(angles,0,193)=314.5463f; MAT_ELEM(angles,1,193)=53.2783f;
669  MAT_ELEM(angles,0,194)=22.7214f; MAT_ELEM(angles,1,194)=68.041f;
670  MAT_ELEM(angles,0,195)=310.7214f; MAT_ELEM(angles,1,195)=68.041f;
671  MAT_ELEM(angles,0,196)=36.0f; MAT_ELEM(angles,1,196)=76.3782f;
672  MAT_ELEM(angles,0,197)=324.0f; MAT_ELEM(angles,1,197)=76.3782f;
673  MAT_ELEM(angles,0,198)=242.5463f; MAT_ELEM(angles,1,198)=53.278f;
674  MAT_ELEM(angles,0,199)=98.5463f; MAT_ELEM(angles,1,199)=53.278f;
675  MAT_ELEM(angles,0,200)=170.5463f; MAT_ELEM(angles,1,200)=53.278f;
676  MAT_ELEM(angles,0,201)=336.7264f; MAT_ELEM(angles,1,201)=37.1611f;
677  MAT_ELEM(angles,0,202)=48.7264f; MAT_ELEM(angles,1,202)=37.1611f;
678  MAT_ELEM(angles,0,203)=351.0f; MAT_ELEM(angles,1,203)=90.0f;
679  MAT_ELEM(angles,0,204)=63.0f; MAT_ELEM(angles,1,204)=90.0f;
680  MAT_ELEM(angles,0,205)=279.0f; MAT_ELEM(angles,1,205)=90.0f;
681  MAT_ELEM(angles,0,206)=221.1634f; MAT_ELEM(angles,1,206)=66.042f;
682  MAT_ELEM(angles,0,207)=149.1634f; MAT_ELEM(angles,1,207)=66.042f;
683  MAT_ELEM(angles,0,208)=196.498f; MAT_ELEM(angles,1,208)=27.943f;
684  MAT_ELEM(angles,0,209)=268.498f; MAT_ELEM(angles,1,209)=27.943f;
685  MAT_ELEM(angles,0,210)=124.498f; MAT_ELEM(angles,1,210)=27.943f;
686  MAT_ELEM(angles,0,211)=340.498f; MAT_ELEM(angles,1,211)=27.9429f;
687  MAT_ELEM(angles,0,212)=52.498f; MAT_ELEM(angles,1,212)=27.9429f;
688  MAT_ELEM(angles,0,213)=346.0516f; MAT_ELEM(angles,1,213)=81.9568f;
689  MAT_ELEM(angles,0,214)=58.0516f; MAT_ELEM(angles,1,214)=81.9568f;
690  MAT_ELEM(angles,0,215)=274.0516f; MAT_ELEM(angles,1,215)=81.9568f;
691  MAT_ELEM(angles,0,216)=210.8366f; MAT_ELEM(angles,1,216)=66.042f;
692  MAT_ELEM(angles,0,217)=138.8366f; MAT_ELEM(angles,1,217)=66.042f;
693  MAT_ELEM(angles,0,218)=192.7264f; MAT_ELEM(angles,1,218)=37.161f;
694  MAT_ELEM(angles,0,219)=264.7264f; MAT_ELEM(angles,1,219)=37.161f;
695  MAT_ELEM(angles,0,220)=120.7264f; MAT_ELEM(angles,1,220)=37.161f;
696  MAT_ELEM(angles,0,221)=6.0948f; MAT_ELEM(angles,1,221)=50.7685f;
697  MAT_ELEM(angles,0,222)=78.0948f; MAT_ELEM(angles,1,222)=50.7685f;
698  MAT_ELEM(angles,0,223)=294.0948f; MAT_ELEM(angles,1,223)=50.7685f;
699  MAT_ELEM(angles,0,224)=13.9484f; MAT_ELEM(angles,1,224)=81.9568f;
700  MAT_ELEM(angles,0,225)=85.9484f; MAT_ELEM(angles,1,225)=81.9568f;
701  MAT_ELEM(angles,0,226)=301.9484f; MAT_ELEM(angles,1,226)=81.9568f;
702  MAT_ELEM(angles,0,227)=45.0f; MAT_ELEM(angles,1,227)=90.0f;
703  MAT_ELEM(angles,0,228)=333.0f; MAT_ELEM(angles,1,228)=90.0f;
704  MAT_ELEM(angles,0,229)=239.2736f; MAT_ELEM(angles,1,229)=37.161f;
705  MAT_ELEM(angles,0,230)=95.2736f; MAT_ELEM(angles,1,230)=37.161f;
706  MAT_ELEM(angles,0,231)=167.2736f; MAT_ELEM(angles,1,231)=37.161f;
707  MAT_ELEM(angles,0,232)=324.0f; MAT_ELEM(angles,1,232)=7.9294f;
708  MAT_ELEM(angles,0,233)=36.0f; MAT_ELEM(angles,1,233)=7.9294f;
709  MAT_ELEM(angles,0,234)=332.6069f; MAT_ELEM(angles,1,234)=61.2449f;
710  MAT_ELEM(angles,0,235)=315.3931f; MAT_ELEM(angles,1,235)=61.2449f;
711  MAT_ELEM(angles,0,236)=328.9523f; MAT_ELEM(angles,1,236)=69.9333f;
712  MAT_ELEM(angles,0,237)=319.0477f; MAT_ELEM(angles,1,237)=69.9333f;
713  MAT_ELEM(angles,0,238)=44.6069f; MAT_ELEM(angles,1,238)=61.2449f;
714  MAT_ELEM(angles,0,239)=31.0477f; MAT_ELEM(angles,1,239)=69.9333f;
715  MAT_ELEM(angles,0,240)=40.9523f; MAT_ELEM(angles,1,240)=69.9333f;
716  MAT_ELEM(angles,0,241)=188.6069f; MAT_ELEM(angles,1,241)=61.245f;
717  MAT_ELEM(angles,0,242)=171.3931f; MAT_ELEM(angles,1,242)=61.245f;
718  MAT_ELEM(angles,0,243)=180.0f; MAT_ELEM(angles,1,243)=55.506f;
719  MAT_ELEM(angles,0,244)=260.6069f; MAT_ELEM(angles,1,244)=61.245f;
720  MAT_ELEM(angles,0,245)=252.0f; MAT_ELEM(angles,1,245)=55.506f;
721  MAT_ELEM(angles,0,246)=108.0f; MAT_ELEM(angles,1,246)=55.506f;
722  MAT_ELEM(angles,0,247)=324.0f; MAT_ELEM(angles,1,247)=39.6468f;
723  MAT_ELEM(angles,0,248)=36.0f; MAT_ELEM(angles,1,248)=39.6468f;
724  MAT_ELEM(angles,0,249)=9.299f; MAT_ELEM(angles,1,249)=58.6205f;
725  MAT_ELEM(angles,0,250)=278.701f; MAT_ELEM(angles,1,250)=58.6205f;
726  MAT_ELEM(angles,0,251)=166.1881f; MAT_ELEM(angles,1,251)=83.2609f;
727  MAT_ELEM(angles,0,252)=121.8119f; MAT_ELEM(angles,1,252)=83.2609f;
728  MAT_ELEM(angles,0,253)=81.299f; MAT_ELEM(angles,1,253)=58.6205f;
729  MAT_ELEM(angles,0,254)=193.8119f; MAT_ELEM(angles,1,254)=83.2609f;
730  MAT_ELEM(angles,0,255)=238.1881f; MAT_ELEM(angles,1,255)=83.2609f;
731  MAT_ELEM(angles,0,256)=265.8119f; MAT_ELEM(angles,1,256)=83.2609f;
732  MAT_ELEM(angles,0,257)=94.1881f; MAT_ELEM(angles,1,257)=83.2609f;
733  MAT_ELEM(angles,0,258)=225.299f; MAT_ELEM(angles,1,258)=58.621f;
734  MAT_ELEM(angles,0,259)=134.701f; MAT_ELEM(angles,1,259)=58.621f;
735  MAT_ELEM(angles,0,260)=180.0f; MAT_ELEM(angles,1,260)=23.788f;
736  MAT_ELEM(angles,0,261)=252.0f; MAT_ELEM(angles,1,261)=23.788f;
737  MAT_ELEM(angles,0,262)=108.0f; MAT_ELEM(angles,1,262)=23.788f;
738  MAT_ELEM(angles,0,263)=353.9052f; MAT_ELEM(angles,1,263)=50.7685f;
739  MAT_ELEM(angles,0,264)=65.9052f; MAT_ELEM(angles,1,264)=50.7685f;
740  MAT_ELEM(angles,0,265)=281.9052f; MAT_ELEM(angles,1,265)=50.7685f;
741  MAT_ELEM(angles,0,266)=9.0f; MAT_ELEM(angles,1,266)=90.0f;
742  MAT_ELEM(angles,0,267)=81.0f; MAT_ELEM(angles,1,267)=90.0f;
743  MAT_ELEM(angles,0,268)=297.0f; MAT_ELEM(angles,1,268)=90.0f;
744  MAT_ELEM(angles,0,269)=229.9484f; MAT_ELEM(angles,1,269)=81.9568f;
745  MAT_ELEM(angles,0,270)=157.9484f; MAT_ELEM(angles,1,270)=81.9568f;
746  MAT_ELEM(angles,0,271)=235.502f; MAT_ELEM(angles,1,271)=27.943f;
747  MAT_ELEM(angles,0,272)=91.502f; MAT_ELEM(angles,1,272)=27.943f;
748  MAT_ELEM(angles,0,273)=163.502f; MAT_ELEM(angles,1,273)=27.943f;
749  MAT_ELEM(angles,0,274)=19.502f; MAT_ELEM(angles,1,274)=27.9429f;
750  MAT_ELEM(angles,0,275)=307.502f; MAT_ELEM(angles,1,275)=27.9429f;
751  MAT_ELEM(angles,0,276)=354.8366f; MAT_ELEM(angles,1,276)=66.0423f;
752  MAT_ELEM(angles,0,277)=66.8366f; MAT_ELEM(angles,1,277)=66.0423f;
753  MAT_ELEM(angles,0,278)=282.8366f; MAT_ELEM(angles,1,278)=66.0423f;
754  MAT_ELEM(angles,0,279)=202.0516f; MAT_ELEM(angles,1,279)=81.9568f;
755  MAT_ELEM(angles,0,280)=130.0516f; MAT_ELEM(angles,1,280)=81.9568f;
756  MAT_ELEM(angles,0,281)=209.9052f; MAT_ELEM(angles,1,281)=50.768f;
757  MAT_ELEM(angles,0,282)=137.9052f; MAT_ELEM(angles,1,282)=50.768f;
758  MAT_ELEM(angles,0,283)=23.2736f; MAT_ELEM(angles,1,283)=37.1611f;
759  MAT_ELEM(angles,0,284)=311.2736f; MAT_ELEM(angles,1,284)=37.1611f;
760  MAT_ELEM(angles,0,285)=5.1634f; MAT_ELEM(angles,1,285)=66.0423f;
761  MAT_ELEM(angles,0,286)=77.1634f; MAT_ELEM(angles,1,286)=66.0423f;
762  MAT_ELEM(angles,0,287)=293.1634f; MAT_ELEM(angles,1,287)=66.0423f;
763  MAT_ELEM(angles,0,288)=27.0f; MAT_ELEM(angles,1,288)=90.0f;
764  MAT_ELEM(angles,0,289)=315.0f; MAT_ELEM(angles,1,289)=90.0f;
765  MAT_ELEM(angles,0,290)=222.0948f; MAT_ELEM(angles,1,290)=50.768f;
766  MAT_ELEM(angles,0,291)=150.0948f; MAT_ELEM(angles,1,291)=50.768f;
767  MAT_ELEM(angles,0,292)=324.0f; MAT_ELEM(angles,1,292)=23.7881f;
768  MAT_ELEM(angles,0,293)=36.0f; MAT_ELEM(angles,1,293)=23.7881f;
769  MAT_ELEM(angles,0,294)=350.701f; MAT_ELEM(angles,1,294)=58.6205f;
770  MAT_ELEM(angles,0,295)=297.299f; MAT_ELEM(angles,1,295)=58.6205f;
771  MAT_ELEM(angles,0,296)=337.8119f; MAT_ELEM(angles,1,296)=83.2609f;
772  MAT_ELEM(angles,0,297)=310.1881f; MAT_ELEM(angles,1,297)=83.2609f;
773  MAT_ELEM(angles,0,298)=62.701f; MAT_ELEM(angles,1,298)=58.6205f;
774  MAT_ELEM(angles,0,299)=22.1881f; MAT_ELEM(angles,1,299)=83.2609f;
775  MAT_ELEM(angles,0,300)=49.8119f; MAT_ELEM(angles,1,300)=83.2609f;
776  MAT_ELEM(angles,0,301)=206.701f; MAT_ELEM(angles,1,301)=58.621f;
777  MAT_ELEM(angles,0,302)=153.299f; MAT_ELEM(angles,1,302)=58.621f;
778  MAT_ELEM(angles,0,303)=180.0f; MAT_ELEM(angles,1,303)=39.647f;
779  MAT_ELEM(angles,0,304)=252.0f; MAT_ELEM(angles,1,304)=39.647f;
780  MAT_ELEM(angles,0,305)=108.0f; MAT_ELEM(angles,1,305)=39.647f;
781  MAT_ELEM(angles,0,306)=324.0f; MAT_ELEM(angles,1,306)=55.5056f;
782  MAT_ELEM(angles,0,307)=36.0f; MAT_ELEM(angles,1,307)=55.5056f;
783  MAT_ELEM(angles,0,308)=27.3931f; MAT_ELEM(angles,1,308)=61.2449f;
784  MAT_ELEM(angles,0,309)=175.0477f; MAT_ELEM(angles,1,309)=69.933f;
785  MAT_ELEM(angles,0,310)=112.9523f; MAT_ELEM(angles,1,310)=69.933f;
786  MAT_ELEM(angles,0,311)=184.9523f; MAT_ELEM(angles,1,311)=69.933f;
787  MAT_ELEM(angles,0,312)=247.0477f; MAT_ELEM(angles,1,312)=69.933f;
788  MAT_ELEM(angles,0,313)=256.9523f; MAT_ELEM(angles,1,313)=69.933f;
789  MAT_ELEM(angles,0,314)=103.0477f; MAT_ELEM(angles,1,314)=69.933f;
790  MAT_ELEM(angles,0,315)=243.3931f; MAT_ELEM(angles,1,315)=61.245f;
791  MAT_ELEM(angles,0,316)=116.6069f; MAT_ELEM(angles,1,316)=61.245f;
792  MAT_ELEM(angles,0,317)=180.0f; MAT_ELEM(angles,1,317)=7.929f;
793  MAT_ELEM(angles,0,318)=252.0f; MAT_ELEM(angles,1,318)=7.929f;
794  MAT_ELEM(angles,0,319)=99.3931f; MAT_ELEM(angles,1,319)=61.245f;
795  MAT_ELEM(angles,0,320)=108.0f; MAT_ELEM(angles,1,320)=7.929f;
796  }
797  else
798  {
799  angles.initZeros(2,81);
800  MAT_ELEM(angles, 0, 0) = 0.000000f; MAT_ELEM(angles, 1, 0) = 0.000000f;
801  MAT_ELEM(angles, 0, 1) = 36.000000f; MAT_ELEM(angles, 1, 1) = 15.858741f;
802  MAT_ELEM(angles, 0, 2) = 36.000000f; MAT_ELEM(angles, 1, 2) = 31.717482f;
803  MAT_ELEM(angles, 0, 3) = 36.000000f; MAT_ELEM(angles, 1, 3) = 47.576224f;
804  MAT_ELEM(angles, 0, 4) = 36.000000f; MAT_ELEM(angles, 1, 4) = 63.434965f;
805  MAT_ELEM(angles, 0, 5) = 62.494295f; MAT_ELEM(angles, 1, 5) = -76.558393f;
806  MAT_ELEM(angles, 0, 6) = 54.000000f; MAT_ELEM(angles, 1, 6) = 90.000000f;
807  MAT_ELEM(angles, 0, 7) = 45.505705f; MAT_ELEM(angles, 1, 7) = 76.558393f;
808  MAT_ELEM(angles, 0, 8) = 108.000000f; MAT_ELEM(angles, 1, 8) = 15.858741f;
809  MAT_ELEM(angles, 0, 9) = 108.000000f; MAT_ELEM(angles, 1, 9) = 31.717482f;
810  MAT_ELEM(angles, 0, 10) = 108.000000f; MAT_ELEM(angles, 1, 10) = 47.576224f;
811  MAT_ELEM(angles, 0, 11) = 108.000000f; MAT_ELEM(angles, 1, 11) = 63.434965f;
812  MAT_ELEM(angles, 0, 12) = 134.494295f; MAT_ELEM(angles, 1, 12) = -76.558393f;
813  MAT_ELEM(angles, 0, 13) = 126.000000f; MAT_ELEM(angles, 1, 13) = 90.000000f;
814  MAT_ELEM(angles, 0, 14) = 117.505705f; MAT_ELEM(angles, 1, 14) = 76.558393f;
815  MAT_ELEM(angles, 0, 15) = 144.000000f; MAT_ELEM(angles, 1, 15) = -15.858741f;
816  MAT_ELEM(angles, 0, 16) = 144.000000f; MAT_ELEM(angles, 1, 16) = -31.717482f;
817  MAT_ELEM(angles, 0, 17) = 144.000000f; MAT_ELEM(angles, 1, 17) = -47.576224f;
818  MAT_ELEM(angles, 0, 18) = 144.000000f; MAT_ELEM(angles, 1, 18) = -63.434965f;
819  MAT_ELEM(angles, 0, 19) = 170.494295f; MAT_ELEM(angles, 1, 19) = 76.558393f;
820  MAT_ELEM(angles, 0, 20) = 162.000000f; MAT_ELEM(angles, 1, 20) = 90.000000f;
821  MAT_ELEM(angles, 0, 21) = 153.505705f; MAT_ELEM(angles, 1, 21) = -76.558393f;
822  MAT_ELEM(angles, 0, 22) = 72.000000f; MAT_ELEM(angles, 1, 22) = -15.858741f;
823  MAT_ELEM(angles, 0, 23) = 72.000000f; MAT_ELEM(angles, 1, 23) = -31.717482f;
824  MAT_ELEM(angles, 0, 24) = 72.000000f; MAT_ELEM(angles, 1, 24) = -47.576224f;
825  MAT_ELEM(angles, 0, 25) = 72.000000f; MAT_ELEM(angles, 1, 25) = -63.434965f;
826  MAT_ELEM(angles, 0, 26) = 98.494295f; MAT_ELEM(angles, 1, 26) = 76.558393f;
827  MAT_ELEM(angles, 0, 27) = 90.000000f; MAT_ELEM(angles, 1, 27) = 90.000000f;
828  MAT_ELEM(angles, 0, 28) = 81.505705f; MAT_ELEM(angles, 1, 28) = -76.558393f;
829  MAT_ELEM(angles, 0, 29) = 0.000000f; MAT_ELEM(angles, 1, 29) = -15.858741f;
830  MAT_ELEM(angles, 0, 30) = 0.000000f; MAT_ELEM(angles, 1, 30) = -31.717482f;
831  MAT_ELEM(angles, 0, 31) = 0.000000f; MAT_ELEM(angles, 1, 31) = -47.576224f;
832  MAT_ELEM(angles, 0, 32) = 0.000000f; MAT_ELEM(angles, 1, 32) = -63.434965f;
833  MAT_ELEM(angles, 0, 33) = 26.494295f; MAT_ELEM(angles, 1, 33) = 76.558393f;
834  MAT_ELEM(angles, 0, 34) = 18.000000f; MAT_ELEM(angles, 1, 34) = 90.000000f;
835  MAT_ELEM(angles, 0, 35) = 9.505705f; MAT_ELEM(angles, 1, 35) = -76.558393f;
836  MAT_ELEM(angles, 0, 36) = 12.811021f; MAT_ELEM(angles, 1, 36) = 42.234673f;
837  MAT_ELEM(angles, 0, 37) = 18.466996f; MAT_ELEM(angles, 1, 37) = 59.620797f;
838  MAT_ELEM(angles, 0, 38) = 0.000000f; MAT_ELEM(angles, 1, 38) = 90.000000f;
839  MAT_ELEM(angles, 0, 39) = 8.867209f; MAT_ELEM(angles, 1, 39) = 75.219088f;
840  MAT_ELEM(angles, 0, 40) = 72.000000f; MAT_ELEM(angles, 1, 40) = 26.565058f;
841  MAT_ELEM(angles, 0, 41) = 59.188979f; MAT_ELEM(angles, 1, 41) = 42.234673f;
842  MAT_ELEM(angles, 0, 42) = 84.811021f; MAT_ELEM(angles, 1, 42) = 42.234673f;
843  MAT_ELEM(angles, 0, 43) = 53.533003f; MAT_ELEM(angles, 1, 43) = 59.620797f;
844  MAT_ELEM(angles, 0, 44) = 72.000000f; MAT_ELEM(angles, 1, 44) = 58.282544f;
845  MAT_ELEM(angles, 0, 45) = 90.466996f; MAT_ELEM(angles, 1, 45) = 59.620797f;
846  MAT_ELEM(angles, 0, 46) = 72.000000f; MAT_ELEM(angles, 1, 46) = 90.000000f;
847  MAT_ELEM(angles, 0, 47) = 63.132791f; MAT_ELEM(angles, 1, 47) = 75.219088f;
848  MAT_ELEM(angles, 0, 48) = 80.867209f; MAT_ELEM(angles, 1, 48) = 75.219088f;
849  MAT_ELEM(angles, 0, 49) = 144.000000f; MAT_ELEM(angles, 1, 49) = 26.565058f;
850  MAT_ELEM(angles, 0, 50) = 131.188979f; MAT_ELEM(angles, 1, 50) = 42.234673f;
851  MAT_ELEM(angles, 0, 51) = 156.811021f; MAT_ELEM(angles, 1, 51) = 42.234673f;
852  MAT_ELEM(angles, 0, 52) = 125.533003f; MAT_ELEM(angles, 1, 52) = 59.620797f;
853  MAT_ELEM(angles, 0, 53) = 144.000000f; MAT_ELEM(angles, 1, 53) = 58.282544f;
854  MAT_ELEM(angles, 0, 54) = 162.466996f; MAT_ELEM(angles, 1, 54) = 59.620797f;
855  MAT_ELEM(angles, 0, 55) = 144.000000f; MAT_ELEM(angles, 1, 55) = 90.000000f;
856  MAT_ELEM(angles, 0, 56) = 135.132791f; MAT_ELEM(angles, 1, 56) = 75.219088f;
857  MAT_ELEM(angles, 0, 57) = 152.867209f; MAT_ELEM(angles, 1, 57) = 75.219088f;
858  MAT_ELEM(angles, 0, 58) = 180.000000f; MAT_ELEM(angles, 1, 58) = -26.565058f;
859  MAT_ELEM(angles, 0, 59) = 167.188979f; MAT_ELEM(angles, 1, 59) = -42.234673f;
860  MAT_ELEM(angles, 0, 60) = 180.000000f; MAT_ELEM(angles, 1, 60) = -58.282544f;
861  MAT_ELEM(angles, 0, 61) = 161.533003f; MAT_ELEM(angles, 1, 61) = -59.620797f;
862  MAT_ELEM(angles, 0, 62) = 171.132791f; MAT_ELEM(angles, 1, 62) = -75.219088f;
863  MAT_ELEM(angles, 0, 63) = 108.000000f; MAT_ELEM(angles, 1, 63) = -26.565058f;
864  MAT_ELEM(angles, 0, 64) = 120.811021f; MAT_ELEM(angles, 1, 64) = -42.234673f;
865  MAT_ELEM(angles, 0, 65) = 95.188979f; MAT_ELEM(angles, 1, 65) = -42.234673f;
866  MAT_ELEM(angles, 0, 66) = 126.466996f; MAT_ELEM(angles, 1, 66) = -59.620797f;
867  MAT_ELEM(angles, 0, 67) = 108.000000f; MAT_ELEM(angles, 1, 67) = -58.282544f;
868  MAT_ELEM(angles, 0, 68) = 89.533003f; MAT_ELEM(angles, 1, 68) = -59.620797f;
869  MAT_ELEM(angles, 0, 69) = 108.000000f; MAT_ELEM(angles, 1, 69) = 90.000000f;
870  MAT_ELEM(angles, 0, 70) = 116.867209f; MAT_ELEM(angles, 1, 70) = -75.219088f;
871  MAT_ELEM(angles, 0, 71) = 99.132791f; MAT_ELEM(angles, 1, 71) = -75.219088f;
872  MAT_ELEM(angles, 0, 72) = 36.000000f; MAT_ELEM(angles, 1, 72) = -26.565058f;
873  MAT_ELEM(angles, 0, 73) = 48.811021f; MAT_ELEM(angles, 1, 73) = -42.234673f;
874  MAT_ELEM(angles, 0, 74) = 23.188979f; MAT_ELEM(angles, 1, 74) = -42.234673f;
875  MAT_ELEM(angles, 0, 75) = 54.466996f; MAT_ELEM(angles, 1, 75) = -59.620797f;
876  MAT_ELEM(angles, 0, 76) = 36.000000f; MAT_ELEM(angles, 1, 76) = -58.282544f;
877  MAT_ELEM(angles, 0, 77) = 17.533003f; MAT_ELEM(angles, 1, 77) = -59.620797f;
878  MAT_ELEM(angles, 0, 78) = 36.000000f; MAT_ELEM(angles, 1, 78) = 90.000000f;
879  MAT_ELEM(angles, 0, 79) = 44.867209f; MAT_ELEM(angles, 1, 79) = -75.219088f;
880  MAT_ELEM(angles, 0, 80) = 27.132791f; MAT_ELEM(angles, 1, 80) = -75.219088f;
881  }
882 
883  angles *= PI/180.0;
884 }
885 
886 
887 void ProgAngResAlign::readData(MultidimArray<double> &half1, MultidimArray<double> &half2)
888 {
889  std::cout << "Reading data..." << std::endl;
890  Image<double> imgHalf1, imgHalf2, inMask;
891  imgHalf1.read(fnhalf1);
892  imgHalf2.read(fnhalf2);
893 
894  half1 = imgHalf1();
895  half2 = imgHalf2();
896 
897  maxRadius = XSIZE(half1)*0.5;
898 
899  if (fnmask!= "")
900  {
901  Image<double> inMask;
902  inMask.read(fnmask);
903  auto mapMask = inMask();
904 
905  // Check maximum radius
906  size_t xdim = XSIZE(mapMask)*0.5;
907  size_t ydim = YSIZE(mapMask)*0.5;
908  size_t zdim = ZSIZE(mapMask)*0.5;
909 
910  size_t Nelems;
911  if (isHelix)
912  {
913  Nelems = round(sqrt(xdim*xdim + ydim*ydim)+1);
914  }
915  else
916  {
917  Nelems = round(sqrt(xdim*xdim + ydim*ydim + zdim*zdim)+1);
918  }
919 
920  MultidimArray<double> radiusElems, rr, radiusTheoretical;
921  radiusElems.initZeros(Nelems);
922  rr.initZeros(Nelems);
923  radiusTheoretical.initZeros(Nelems);
924 
925 
926  long n=0;
927  if (isHelix)
928  {
929  for(size_t k=0; k<ZSIZE(mapMask); ++k)
930  {
931  for(size_t i=0; i<YSIZE(mapMask); ++i)
932  {
933  size_t yy = (i-ydim);
934  yy *= yy;
935 
936  for(size_t j=0; j<XSIZE(mapMask); ++j)
937  {
938  size_t xx = (j-xdim);
939  xx *= xx;
940  int rad = round(sqrt(xx + yy));
941  auto idx = (int) round(rad);
942  dAi(radiusTheoretical, idx) += 1;
943 
944  if (DIRECT_MULTIDIM_ELEM(mapMask, n)>0)
945  {
946 
947  dAi(radiusElems, idx) += 1;
948 
949  dAi(rr, idx) = rad;
950  }
951  ++n;
952  }
953  }
954  }
955  }
956  else
957  {
958  for(size_t k=0; k<ZSIZE(mapMask); ++k)
959  {
960  size_t zz = (k-zdim);
961  zz *= zz;
962 
963  for(size_t i=0; i<YSIZE(mapMask); ++i)
964  {
965  size_t yy = (i-ydim);
966  yy *= yy;
967 
968  for(size_t j=0; j<XSIZE(mapMask); ++j)
969  {
970  size_t xx = (j-xdim);
971  xx *= xx;
972  int rad = round(sqrt(xx + yy + zz));
973  auto idx = (int) round(rad);
974  dAi(radiusTheoretical, idx) += 1;
975 
976  if (DIRECT_MULTIDIM_ELEM(mapMask, n)>0)
977  {
978 
979  dAi(radiusElems, idx) += 1;
980 
981  dAi(rr, idx) = rad;
982  }
983  ++n;
984  }
985  }
986  }
987  }
988 
989 
990  if (limRad)
991  {
992  double radiusThreshold;
993 
994  radiusThreshold = 0.75;
995 
996 
997  for (size_t k=0; k<radiusElems.nzyxdim; k++)
998  {
999  if (dAi(radiusElems, k) > radiusThreshold * dAi(radiusTheoretical, k))
1000  maxRadius = k;
1001  }
1002 
1003  std::cout << "maxRadius = -> " << maxRadius << std::endl;
1004  }
1005 
1006 
1007  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mapMask)
1008  {
1009 
1010  double valmask = DIRECT_MULTIDIM_ELEM(mapMask, n);
1011  DIRECT_MULTIDIM_ELEM(half1, n) = DIRECT_MULTIDIM_ELEM(half1, n) * valmask;
1012  DIRECT_MULTIDIM_ELEM(half2, n) = DIRECT_MULTIDIM_ELEM(half2, n) * valmask;
1013  }
1014  }
1015 }
1016 
1017 
1018 void ProgAngResAlign::generateShellMask(MultidimArray<double> &shellMask, size_t shellNum, bool ishelix)
1019 {
1020  double sigma2 = 5*5;
1021  double inv2sigma2;
1022  inv2sigma2 = 1/(2*sigma2);
1023  if (ishelix)
1024  {
1025  FOR_ALL_ELEMENTS_IN_ARRAY3D(shellMask)
1026  {
1027  double radius = sqrt(i*i + j*j);
1028  radius = radius - shellNum;
1029  A3D_ELEM(shellMask, k, i, j) = exp(-(radius*radius)*inv2sigma2);
1030  }
1031  }else
1032  {
1033  FOR_ALL_ELEMENTS_IN_ARRAY3D(shellMask)
1034  {
1035  double radius = sqrt(k*k + i*i + j*j);
1036  radius = radius - shellNum;
1037  A3D_ELEM(shellMask, k, i, j) = exp(-(radius*radius)*inv2sigma2);
1038  }
1039  }
1040 }
1041 
1042 
1043 void ProgAngResAlign::applyShellMask(const MultidimArray<double> &half1, const MultidimArray<double> &half2,
1044  const MultidimArray<double> &shellMask,
1045  MultidimArray<double> &half1_aux,
1046  MultidimArray<double> &half2_aux)
1047 {
1048  half1_aux = half1;
1049  half2_aux = half2;
1050 
1051  // Applying the mask
1052  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(shellMask)
1053  {
1054  double valmask = DIRECT_MULTIDIM_ELEM(shellMask, n);
1055  DIRECT_MULTIDIM_ELEM(half1_aux, n) = DIRECT_MULTIDIM_ELEM(half1, n) * valmask;
1056  DIRECT_MULTIDIM_ELEM(half2_aux, n) = DIRECT_MULTIDIM_ELEM(half2, n) * valmask;
1057  }
1058 
1059  half1_aux.setXmippOrigin();
1060  half2_aux.setXmippOrigin();
1061 
1062  // fftshift must by applied before computing the fft. This will be computed later
1063  CenterFFT(half1_aux, true);
1064  CenterFFT(half2_aux, true);
1065 }
1066 
1067 
1068 void ProgAngResAlign::run()
1069  {
1070  std::cout << "Starting ... " << std::endl << std::endl;
1071 
1072  //This read the data and applies a fftshift to in the next step compute the fft
1073  // The maps half1_aux, and half2_aux are used to recover the initial maps after masking in each iteration
1074  MultidimArray<double> half1, half2;
1075  MultidimArray<double> &phalf1 = half1, &phalf2 = half2;
1076  readData(phalf1, phalf2);
1077 
1078  // Defining frequencies freq_fourier_x,y,z and freqMap
1079  // The number of frequencies in each shell freqElem is determined
1080  defineFrequenciesSimple(phalf1);
1081 
1082  //Define shell Mask
1083  MultidimArray<double> shellMask;
1084  shellMask.resizeNoCopy(phalf1);
1085 
1086  // Generating the set of directions to be analyzed
1087  // And converting the cone angle to radians
1088  generateDirections(angles, false);
1089  ang_con = ang_con*PI/180.0;
1090 
1091  MultidimArray<double> half1_aux;
1092  MultidimArray<double> half2_aux;
1093 
1094  //Generate shell mask and computing radial resolution
1095 
1096  MetaDataVec mdOut;
1097 
1098  //std::cout << "max radius = " << maxRadius << std::endl;
1099 
1100  for (size_t shellRadius = 0; shellRadius<maxRadius; shellRadius++)
1101  {
1102  shellMask.initZeros();
1103  shellMask.setXmippOrigin();
1104  generateShellMask(shellMask, shellRadius, isHelix);
1105 
1106  // Image<double> img;
1107  // img() = shellMask;
1108  // FileName fn;
1109  // fn = formatString("mk_%i.mrc", shellRadius);
1110  // img.write(fn);
1111 
1112  applyShellMask(half1, half2, shellMask, half1_aux, half2_aux);
1113 
1114  //Computing the FFT
1115  FourierTransformer transformer2(FFTW_BACKWARD), transformer1(FFTW_BACKWARD);
1116  //transformer1.setThreadsNumber(Nthreads);
1117  //transformer2.setThreadsNumber(Nthreads);
1118 
1119  transformer1.FourierTransform(half1_aux, FT1, false);
1120  transformer2.FourierTransform(half2_aux, FT2, false);
1121 
1122  // Storing the shell of both maps as vectors global
1123  // The global FSC is also computed
1124 
1125  arrangeFSC_and_fscGlobal();
1126 
1127 
1128  double res;
1129  if (directionalRes)
1130  {
1131  std::vector<double> radialResolution(angles.mdimx);
1132 
1133  for (size_t k = 0; k<angles.mdimx; k++)
1134  {
1135  float rot = MAT_ELEM(angles, 0, k);
1136  float tilt = MAT_ELEM(angles, 1, k);
1137 
1138  double resInterp = -1;
1139  MultidimArray<float> fsc;
1140 
1141  // Estimating the direction FSC along the direction given by rot and tilt
1142  fscDir_fast(fsc, rot, tilt, thrs, resInterp);
1143 
1144  radialResolution[k] = resInterp;
1145  }
1146 
1147  std::sort(radialResolution.begin(),radialResolution.end());
1148 
1149  res = radialResolution[round(0.5*angles.mdimx)];
1150  radialResolution.clear();
1151  }
1152  else
1153  {
1154  fscGlobal(thrs, res);
1155  }
1156 
1157  std::cout << "Radius "<< shellRadius << "/" << maxRadius << " " << res << "A" << std::endl;
1158 
1159  MDRowVec row;
1160  row.setValue(MDL_RESOLUTION_FRC, res);
1161  row.setValue(MDL_IDX, shellRadius);
1162  mdOut.addRow(row);
1163 
1164  }
1165 
1166  mdOut.write(fnOut);
1167 }
dAi
#define dAi(v, i)
Definition: multidim_array_base.h:530
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
min
void min(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:118
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
alglib::conj
alglib::complex conj(const alglib::complex &z)
Definition: ap.cpp:4988
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
MDRowVec::setValue
void setValue(const MDObject &object) override
Definition: metadata_row_vec.cpp:199
MultidimArray::initConstant
void initConstant(T val)
Definition: multidim_array.h:2705
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
angular_resolution_alignment.h
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
MetaDataVec::addRow
size_t addRow(const MDRow &row) override
Definition: metadata_vec.cpp:225
XmippProgram::addSeeAlsoLine
void addSeeAlsoLine(const char *seeAlso)
Definition: xmipp_program.cpp:290
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
CenterFFT
void CenterFFT(MultidimArray< T > &v, bool forward)
Definition: xmipp_fft.h:291
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
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
ctpl_stl.h
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
f
double * f
Definition: numerical_recipes.cpp:4414
fourier_filter.h
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
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
dAkij
#define dAkij(V, k, i, j)
Definition: multidim_array_base.h:244
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
XmippProgram::addExampleLine
void addExampleLine(const char *example, bool verbatim=true)
Definition: xmipp_program.cpp:286
monogenic_signal.h
NZYXSIZE
#define NZYXSIZE(v)
Definition: multidim_array_base.h:119
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
sort
void sort(struct DCEL_T *dcel)
Definition: sorting.cpp:18
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
j
#define j
Definition: numerical_recipes.cpp:2493
FOR_ALL_ELEMENTS_IN_ARRAY1D
#define FOR_ALL_ELEMENTS_IN_ARRAY1D(v)
Definition: multidim_array_base.h:554
alglib::round
int round(double x)
Definition: ap.cpp:7245
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
MDL_IDX
Index within a list (size_t)
Definition: metadata_label.h:245
u
doublereal * u
Definition: numerical_recipes.cpp:2230
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
ImageBase::read
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
Definition: xmipp_image_base.cpp:119
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
XmippProgram::addUsageLine
void addUsageLine(const char *line, bool verbatim=false)
Definition: xmipp_program.cpp:282
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
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
n
int * n
Definition: numerical_recipes.cpp:2229
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305
MDL_RESOLUTION_FRC
Fourier shell correlation (double)
Definition: metadata_label.h:386
Generated by   doxygen 1.8.13