Xmipp  v3.23.11-Nereus
resolution_fso.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Jose Luis Vilas (joseluis.vilas-prieto@yale.edu)
4  * or (jlvilas@cnb.csic.es)
5  * Hemant. D. Tagare (hemant.tagare@yale.edu)
6  *
7  * Yale University, New Haven, Connecticut, United States of America
8  *
9  * This program is free software; you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation; either version 2 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program; if not, write to the Free Software
21  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
22  * 02111-1307 USA
23  *
24  * All comments concerning this program package may be sent to the
25  * e-mail address 'xmipp@cnb.csic.es'
26  ***************************************************************************/
27 
28 #include "resolution_fso.h"
29 #include <core/metadata_extension.h>
30 #include <data/monogenic_signal.h>
31 #include <data/fourier_filter.h>
32 #include <random>
33 #include <limits>
34 #include <CTPL/ctpl_stl.h>
35 
36 void ProgFSO::defineParams()
37 {
38  addUsageLine("Calculate Fourier Shell Occupancy - FSO curve - via directional FSC measurements.");
39  addUsageLine("Following outputs are generated:");
40  addUsageLine(" 1) FSO curve");
41  addUsageLine(" 2) Global resolution from FSO and FSC");
42  addUsageLine(" 3) 3DFSC");
43  addUsageLine(" 4) Anisotropic filter");
44  addUsageLine("Reference: J.L. Vilas, H.D. Tagare, XXXXX (2021)");
45  addUsageLine("+* Fourier Shell Occupancy (FSO)", true);
46  addUsageLine("+ The Fourier Shell Occupancy Curve can be obtained from a set of directional FSC (see below).");
47  addUsageLine("+ To do that, two half maps are used to determine the Global FSC at threshold 0.143. Then, the ratio between the number");
48  addUsageLine("+ of directions with resolution higher (better) than the Global resolution and the total number of measured directions is");
49  addUsageLine("+ calculated at different frequencies (resolutions). Note that this ratio is between 0 (resolution of all directions is worse than the global FSC)");
50  addUsageLine("+ resolution than the global FSC) and 1 (all directions present better resolution than the FSC) at a given resolution.");
51  addUsageLine("+ In the particular case, FSO curve takes the value of 0.5 (FSO=0.5), then half of the directions are better, and.");
52  addUsageLine("+ the other half are worse than the FSC, this situation occurs at the resoltuion of the map. It means the FSO = 0.5 at the ");
53  addUsageLine("+ FSC resolution. A map is isotropic if all directional resolution are similar, and anisotropic is there are significant resolution values along");
54  addUsageLine("+ different directions. Thus, when the FSO presents a sharp cliff, a step-like function, the map will be isotropic.");
55  addUsageLine("+ In contrast, when the OFSC shows a slope the map will be anisotropic. The lesser slope the higher resolution isotropy.");
56  addUsageLine("+ ");
57  addUsageLine("+* Directional Fourier Shell Correlation (dFSC)", true);
58  addUsageLine("+ This program estimates the directional FSC between two half maps along all posible directions on the projection sphere.");
59  addUsageLine("+ The directionality is measured by means of conical-like filters in Fourier Space. To avoid possible Gibbs effects ");
60  addUsageLine("+ the filters are gaussian functions with their respective maxima along the filtering direction. A set of 321 directions ");
61  addUsageLine("+ is used to cover the projection sphere, computing for each direction the directional FSC at 0.143 between the two half maps.");
62  addUsageLine("+ The result is a set of 321 FSC curves (321 is the number of analyzed directions, densely covering the projection sphere).");
63  addUsageLine("+ The 3DFSC is then obtained from all curves by interpolation. Note that as well as it occurs with global FSC, the directional FSC is mask dependent.");
64  addUsageLine(" ");
65  addUsageLine("+* Resolution Distribution and 3DFSC", true);
66  addUsageLine("+ The directional-FSC, dFSC is estimated along 321 directions on the projection sphere. For each direction the corresponding");
67  addUsageLine("+ resolution is determined. Thus, it is possible to determine the resolution distribution on the projection sphere.");
68  addUsageLine("+ This distribution is saved in the output metadata named resolutionDistribution.xmd. Also by means of all dFSC, the 3DFSC");
69  addUsageLine("+ is calculated and saved as 3dFSC.mrc, which gives an idea about the information distributionin Fourier Space.");
70  addUsageLine(" ");
71  addUsageLine(" ");
72  addSeeAlsoLine("resolution_fsc");
73 
74  addParamsLine(" --half1 <input_file> : Input Half map 1");
75  addParamsLine(" --half2 <input_file> : Input Half map 2");
76 
77  addParamsLine(" [-o <output_folder=\"\">] : Folder where the results will be stored.");
78 
79  addParamsLine(" [--sampling <Ts=1>] : (Optical) Pixel size (Angstrom). If it is not provided by default will be 1 A/px.");
80  addParamsLine(" [--mask <input_file=\"\">] : (Optional) Smooth mask to remove noise. If it is not provided, the computation will be carried out without mask.");
81 
82  addParamsLine(" [--anglecone <ang_con=17>] : (Optional) Angle Cone (angle between the axis and the generatrix) for estimating the directional FSC");
83  addParamsLine(" [--threshold <thrs=0.143>] : (Optional) Threshold for the FSC/directionalFSC estimation ");
84 
85  addParamsLine(" [--threedfsc_filter] : (Optional) Put this flag to estimate the 3DFSC, and apply it as low pass filter to obtain a directionally filtered map. It mean to apply an anisotropic filter.");
86 
87  addParamsLine(" [--threads <Nthreads=1>] : (Optional) Number of threads to be used");
88 
89  addExampleLine("Resolution of two half maps half1.mrc and half2.mrc with a sampling rate of 2 A/px", false);
90  addExampleLine("xmipp_resolution_fso --half1 half1.mrc --half2 half2.mrc --sampling_rate 2 ");
91  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);
92  addExampleLine("xmipp_resolution_fso --half1 half1.mrc --half2 half2.mrc --mask mask.mrc --sampling_rate 2 ");
93 }
94 
95 void ProgFSO::readParams()
96 {
97  fnhalf1 = getParam("--half1");
98  fnhalf2 = getParam("--half2");
99  fnOut = getParam("-o");
100 
101  sampling = getDoubleParam("--sampling");
102  fnmask = getParam("--mask");
103  ang_con = getDoubleParam("--anglecone");
104  thrs = getDoubleParam("--threshold");
105  do_3dfsc_filter = checkParam("--threedfsc_filter");
106 
107  Nthreads = getIntParam("--threads");
108 }
109 
110 
111 void ProgFSO::defineFrequencies(const MultidimArray< std::complex<double> > &mapfftV,
112  const MultidimArray<double> &inputVol)
113 {
114  // Initializing the frequency vectors
115  freq_fourier_z.initZeros(ZSIZE(mapfftV));
116  freq_fourier_x.initZeros(XSIZE(mapfftV));
117  freq_fourier_y.initZeros(YSIZE(mapfftV));
118 
119  // u is the frequency
120  double u;
121 
122  // Defining frequency components. First element should be 0, it is set as the smallest number to avoid singularities
123 
124  VEC_ELEM(freq_fourier_z,0) = std::numeric_limits<double>::min();
125  for(size_t k=1; k<ZSIZE(mapfftV); ++k){
126  FFT_IDX2DIGFREQ(k,ZSIZE(inputVol), u);
127  VEC_ELEM(freq_fourier_z, k) = u;
128  }
129 
130  VEC_ELEM(freq_fourier_y,0) = std::numeric_limits<double>::min();
131  for(size_t k=1; k<YSIZE(mapfftV); ++k){
132  FFT_IDX2DIGFREQ(k,YSIZE(inputVol), u);
133  VEC_ELEM(freq_fourier_y, k) = u;
134  }
135 
136  VEC_ELEM(freq_fourier_x,0) = std::numeric_limits<double>::min();
137  for(size_t k=1; k<XSIZE(mapfftV); ++k){
138  FFT_IDX2DIGFREQ(k,XSIZE(inputVol), u);
139  VEC_ELEM(freq_fourier_x, k) = u;
140  }
141 
142  //Initializing map with frequencies
143  freqMap.resizeNoCopy(mapfftV);
144  freqMap.initConstant(1.9); //Nyquist is 2, we take 1.9 greater than Nyquist
145 
146  xvoldim = XSIZE(inputVol);
147  yvoldim = YSIZE(inputVol);
148  zvoldim = ZSIZE(inputVol);
149  freqElems.initZeros(xvoldim/2+1);
150 
151  // Directional frequencies along each direction
152  double uz, uy, ux, uz2, uz2y2;
153  long n=0;
154  int idx = 0;
155 
156  // Ncomps is the number of frequencies lesser than Nyquist
157  long Ncomps = 0;
158 
159  for(size_t k=0; k<ZSIZE(mapfftV); ++k)
160  {
161  uz = VEC_ELEM(freq_fourier_z, k);
162  uz2 = uz*uz;
163  for(size_t i=0; i<YSIZE(mapfftV); ++i)
164  {
165  uy = VEC_ELEM(freq_fourier_y, i);
166  uz2y2 = uz2 + uy*uy;
167 
168  for(size_t j=0; j<XSIZE(mapfftV); ++j)
169  {
170  ux = VEC_ELEM(freq_fourier_x, j);
171  ux = sqrt(uz2y2 + ux*ux);
172 
173  if (ux<=0.5)
174  {
175  idx = (int) round(ux * xvoldim);
176  ++Ncomps;
177  DIRECT_MULTIDIM_ELEM(freqElems, idx) += 1;
178 
179  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1/ux;
180 
181  if ((j == 0) && (uy<0))
182  {
183  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1.9;
184  DIRECT_MULTIDIM_ELEM(freqElems,idx) -= 1;
185  --Ncomps;
186  }
187 
188  if ((i == 0) && (j == 0) && (uz<0))
189  {
190  DIRECT_MULTIDIM_ELEM(freqMap,n) = 1.9;
191  DIRECT_MULTIDIM_ELEM(freqElems,idx) -= 1;
192  --Ncomps;
193  }
194 
195  }
196  ++n;
197  }
198  }
199  }
200 
201  // Image<double> img;
202  // img() = freqMap;
203  // img.write("freqMap.mrc");
204  real_z1z2.initZeros(Ncomps);
205 }
206 
207 
208 
209 void ProgFSO::arrangeFSC_and_fscGlobal(double sampling_rate,
210  double &thrs, MultidimArray<double> &freq)
211  {
212  // cumpos is the the cumulative number of frequencies per shell number
213  // First shell has 0 elements
214  // second shell has the number of elements of the first shell
215  // Third shell has the number of elements of the first+sencond shells and so on
216  // freqElems in .h and set in defineFreq function
217  cumpos.resizeNoCopy(NZYXSIZE(freqElems));
218 
219  DIRECT_MULTIDIM_ELEM(cumpos,0) = 0;
220  for (long n = 1; n<NZYXSIZE(cumpos); ++n)
221  {
222  DIRECT_MULTIDIM_ELEM(cumpos,n) = DIRECT_MULTIDIM_ELEM(cumpos,n-1) + DIRECT_MULTIDIM_ELEM(freqElems,n-1);
223  }
224 
225  // real_z1z2, absz1_vec and absz2_vec will storage real(z1*conj(z2)), ||z1|| and ||z2|| using the
226  //Fourier Coefficients z1, and z2 of both half maps. The storage will be in order, first elements of
227  // the vector are the freq 0, next freq 1 and so on. Ncomps is the number of components with frequency
228  // lesser than Nyquist (defined in defineFreq)
229 
230  absz1_vec = real_z1z2;
231  absz2_vec = real_z1z2;
232 
233  // num and den are de numerator and denominator of the fsc
234  MultidimArray<long> pos;
235  MultidimArray<double> num, den1, den2, noiseAmp;
236  num.initZeros(freqElems);
237  pos.resizeNoCopy(num);
238 
239  freqidx.resizeNoCopy(real_z1z2);
240 
241  // arr2indx is the position of in Fourier space. It is defined in the .h
242  // Because we only work with the frequencies lesser than nyquist their position
243  // int he Fourier space is lost. arr2indx has the position of each frequency.
244  // It is used at the end of the algorithm to generate the 3DFSC
245  arr2indx.resizeNoCopy(real_z1z2);
246  arr2indx.initZeros();
247  freqidx.initZeros();
248  pos.initZeros();
249  den1 = num;
250  den2 = num;
251 
252 
253  auto ZdimFT1=(int)ZSIZE(FT1);
254  auto YdimFT1=(int)YSIZE(FT1);
255  auto XdimFT1=(int)XSIZE(FT1);
256 
257  // fx, fy, fz, defined in the .h are the frequencies of each pixel with frequency
258  // lesser than Nyquist along each axis. They are used to compute the dot product
259  // between each vector position and the direction of the cone to calculate the
260  // directional FSC
261  fx.resizeNoCopy(real_z1z2);
262  fy.resizeNoCopy(real_z1z2);
263  fz.resizeNoCopy(real_z1z2);
264 
265  long n = 0;
266  for (int k=0; k<ZdimFT1; k++)
267  {
268  double uz = VEC_ELEM(freq_fourier_z, k);
269  for (int i=0; i<YdimFT1; i++)
270  {
271  double uy = VEC_ELEM(freq_fourier_y, i);
272  for (int j=0; j<XdimFT1; j++)
273  {
274  double ux = VEC_ELEM(freq_fourier_x, j);
275 
276  double iun = DIRECT_MULTIDIM_ELEM(freqMap,n);
277  double f = 1/iun;
278  ++n;
279 
280  // Only reachable frequencies
281  // To speed up the algorithm, only are considered those voxels with frequency lesser than Nyquist, 0.5. The vector idx_count
282  // stores all frequencies lesser than Nyquist. This vector determines the frequency of each component of
283  // real_z1z2, absz1_vec, absz2_vec.
284  if (f>0.5)
285  continue;
286 
287  // Index of each frequency
288  auto idx = (int) round(f * xvoldim);
289 
290  int idx_count = DIRECT_MULTIDIM_ELEM(cumpos, idx) + DIRECT_MULTIDIM_ELEM(pos, idx);
291 
292  DIRECT_MULTIDIM_ELEM(arr2indx, idx_count) = n - 1;
293 
294  // Storing normalized frequencies
295  DIRECT_MULTIDIM_ELEM(fx, idx_count) = (float) ux*iun;
296  DIRECT_MULTIDIM_ELEM(fy, idx_count) = (float) uy*iun;
297  DIRECT_MULTIDIM_ELEM(fz, idx_count) = (float) uz*iun;
298 
299  // In this vector we store the index of each Fourier coefficient, thus in the
300  // directional FSC estimation the calculus of the index is avoided speeding up the calculation
301  DIRECT_MULTIDIM_ELEM(freqidx, idx_count) = idx;
302 
303  // Fourier coefficients of both halves
304  std::complex<double> &z1 = dAkij(FT1, k, i, j);
305  std::complex<double> &z2 = dAkij(FT2, k, i, j);
306 
307  double absz1 = abs(z1);
308  double absz2 = abs(z2);
309 
310  DIRECT_MULTIDIM_ELEM(real_z1z2, idx_count) = (float) real(conj(z1)*z2);
311  DIRECT_MULTIDIM_ELEM(absz1_vec, idx_count) = (float) absz1*absz1;
312  DIRECT_MULTIDIM_ELEM(absz2_vec, idx_count) = (float) absz2*absz2;
313 
314  dAi(num,idx) += real(conj(z1) * z2);
315  dAi(den1,idx) += absz1*absz1;
316  dAi(den2,idx) += absz2*absz2;
317 
318  DIRECT_MULTIDIM_ELEM(pos, idx) +=1;
319  }
320  }
321  }
322 
323  // The global FSC is stored as a metadata
324  size_t id;
325  MetaDataVec mdRes;
326  MultidimArray< double > frc;
327  freq.initZeros(freqElems);
328  frc.initZeros(freqElems);
329  FOR_ALL_ELEMENTS_IN_ARRAY1D(freq)
330  {
331  dAi(frc,i) = dAi(num,i)/sqrt( (dAi(den1,i)) * (dAi(den2,i)) );
332  dAi(freq,i) = (float) i / (xvoldim * sampling_rate);
333 
334  if (i>0)
335  {
336  id=mdRes.addObject();
337 
338  mdRes.setValue(MDL_RESOLUTION_FREQ,dAi(freq, i),id);
339  mdRes.setValue(MDL_RESOLUTION_FRC,dAi(frc, i),id);
340  mdRes.setValue(MDL_RESOLUTION_FREQREAL, 1./dAi(freq, i), id);
341  }
342  }
343  mdRes.write(fnOut+"/GlobalFSC.xmd");
344 
345  // Here the FSC at 0.143 is obtained by interpolating
346  fscInterpolation(freq, frc);
347 
348  aniFilter.initZeros(freqidx);
349  }
350 
351 
352 void ProgFSO::fscInterpolation(const MultidimArray<double> &freq, const MultidimArray< double > &frc)
353 {
354  // Here the FSC at 0.143 is obtained by interpolating
355  FOR_ALL_ELEMENTS_IN_ARRAY1D(freq)
356  {
357  double ff = dAi(frc,i);
358  if ( (ff<=thrs) && (i>2) )
359  {
360  double y2, y1, x2, x1, slope, ny;
361  y1 = ff;
362  x1 = dAi(freq,i);
363  y2 = dAi(frc, i-1);
364  x2 = dAi(freq, i-1);
365  slope = (y2 - y1)/(x2 - x1);
366  ny = y2 - slope*x2;
367 
368  double fscResolution;
369  fscResolution = (thrs - ny)/slope;
370  std::cout << "Resolution " << 1/fscResolution << std::endl;
371  break;
372  }
373  }
374 }
375 
376 
377 
378 void ProgFSO::fscDir_fast(MultidimArray<float> &fsc, double rot, double tilt,
379  MultidimArray<float> &threeD_FSC, MultidimArray<float> &normalizationMap,
380  double &thrs, double &resol, std::vector<Matrix2D<double>> &freqMat, size_t dirnum)
381 {
382  // FSCDIR_FAST: computes the directional FSC along a direction given by the angles rot and tilt.
383  // Thus the directional resolution, resol, is estimated with the threshold, thrs.
384  // The direcional FSC is stored in a metadata mdRes and in the multidimarray fsc. Later, this multidimarray
385  // will be used to estimate the FSO.
386  // In addition, the 3dfsc and a normalizationmap (needed to estimate the 3dfsc) are created. For each direction
387  // these vectors are updated
388  size_t dim = NZYXSIZE(freqElems);
389 
390  // numerator and denominator of the fsc
391  MultidimArray<float> num, den1, den2;
392  num.initZeros(dim);
393  den1.initZeros(dim);
394  den2.initZeros(dim);
395 
396  // in vecidx the position n of the vector real_z1z2 is stored.
397  // this allow to have a quick access to the frequencies and the
398  // positions of the threeD_FSC without sweeping all positions
399  std::vector<long> vecidx;
400  vecidx.reserve(real_z1z2.nzyxdim);
401  // the used weight of the directional filter
402  std::vector<double> weightFSC3D;
403  weightFSC3D.reserve(real_z1z2.nzyxdim);
404 
405 
406  //Parameter to determine the direction of the cone
407  float x_dir, y_dir, z_dir, cosAngle, aux;
408  x_dir = sinf(tilt)*cosf(rot);
409  y_dir = sinf(tilt)*sinf(rot);
410  z_dir = cosf(tilt);
411 
412  // For the Bingham Test
413  Matrix2D<double> T;
414  T.initZeros(3,3);
415  MAT_ELEM(T, 0, 0) = x_dir*x_dir; MAT_ELEM(T, 0, 1) = x_dir*y_dir; MAT_ELEM(T, 0, 2) = x_dir*z_dir;
416  MAT_ELEM(T, 1, 0) = y_dir*x_dir; MAT_ELEM(T, 1, 1) = y_dir*y_dir; MAT_ELEM(T, 1, 2) = y_dir*z_dir;
417  MAT_ELEM(T, 2, 0) = z_dir*x_dir; MAT_ELEM(T, 2, 1) = z_dir*y_dir; MAT_ELEM(T, 2, 2) = z_dir*z_dir;
418 
419  cosAngle = (float) cos(ang_con);
420 
421  // It is multiply by 0.5 because later the weight is
422  // cosine = sqrt(exp( -((cosine -1)*(cosine -1))*aux ));
423  // thus the computation of the weight is speeded up
424  // aux = 4.0/((cos(ang_con) -1)*(cos(ang_con) -1));
425  aux = (4.0/((cosAngle -1)*(cosAngle -1)));//*0.5;
426 
427  // Computing directional resolution
428  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(real_z1z2)
429  {
430  // frecuencies of the fourier position
431  float ux = DIRECT_MULTIDIM_ELEM(fx, n);
432  float uy = DIRECT_MULTIDIM_ELEM(fy, n);
433  float uz = DIRECT_MULTIDIM_ELEM(fz, n);
434 
435  // angle between the position and the direction of the cone
436  float cosine = fabs(x_dir*ux + y_dir*uy + z_dir*uz);
437 
438  // Only points inside the cone
439  if (cosine >= cosAngle)
440  {
441  // in fact it is sqrt(exp( -((cosine -1)*(cosine -1))*aux ));
442  // For sake of performance The sqrt is avoided dividing the by 2 in aux
443  cosine = expf( -((cosine -1)*(cosine -1))*aux);
444 
445  vecidx.push_back(n);
446  // cosine *= cosine; Commented because is equivalent to remove the root square in aux
447  weightFSC3D.push_back(cosine);
448 
449  // selecting the frequency of the shell
450  size_t idxf = DIRECT_MULTIDIM_ELEM(freqidx, n);
451 
452  // computing the numerator and denominator
453  dAi(num, idxf) += DIRECT_MULTIDIM_ELEM(real_z1z2, n) * cosine;
454  dAi(den1,idxf) += DIRECT_MULTIDIM_ELEM(absz1_vec, n) * cosine;
455  dAi(den2,idxf) += DIRECT_MULTIDIM_ELEM(absz2_vec, n) * cosine;
456  }
457  }
458 
459  fsc.resizeNoCopy(num);
460  fsc.initConstant(1.0);
461 
462  MetaDataVec mdOut;
463 
464  //The fsc is stored in a metadata and saved
465  bool flagRes = true;
466  FOR_ALL_ELEMENTS_IN_ARRAY1D(num)
467  {
468  double auxfsc = (dAi(num,i))/(sqrt(dAi(den1,i)*dAi(den2,i))+1e-38);
469  dAi(fsc,i) = std::max(0.0, auxfsc);
470 
471  if (dAi(fsc,i)>=thrs)
472  {
473  freqMat[i] = freqMat[i] + T;
474  }
475 
476  // MDRowVec row;
477  // row.setValue(MDL_RESOLUTION_FRC, (double) dAi(fsc,i));
478  // mdOut.addRow(row);
479 
480  if (flagRes && (i>2) && (dAi(fsc,i)<=thrs))
481  {
482  flagRes = false;
483  double ff = (float) i / (xvoldim * sampling); // frequency
484  resol = 1./ff;
485  }
486  }
487  // FileName auxfn;
488  // auxfn = formatString("fscDirection_%i.xmd", dirnum);
489  // mdOut.write(auxfn);
490 
491  // the 3dfsc is computed and updated
492  if (do_3dfsc_filter)
493  {
494  // sizevec is the number of points inside the cone
495  size_t sizevec = vecidx.size();
496  for (size_t kk = 0; kk< sizevec; ++kk)
497  {
498  double w = weightFSC3D[kk]; // the weights
499  long n = vecidx[kk]; // position in fourier space
500  size_t ind = DIRECT_MULTIDIM_ELEM(freqidx, n); // position in the fsc array
501 
502  // 3dfsc considering the weight
503  DIRECT_MULTIDIM_ELEM(threeD_FSC, n) += w*DIRECT_MULTIDIM_ELEM(fsc, ind);
504  // normalization map, we will divide by it once all direction are computed
505  DIRECT_MULTIDIM_ELEM(normalizationMap, n) += w;
506  }
507  }
508 }
509 
510 double ProgFSO::incompleteGammaFunction(double &x)
511 {
512  size_t idx;
513  idx = round(2*x);
514  if (idx>40)
515  idx = 40;
516  if (idx<0)
517  idx = 0;
518 
519  float val;
520 
521  // Table with the values of the incomplete lower gamma function. The set of values of the table can be
522  // obtained in matlab with the function gammainc(x,5). The implementation of this funcitonis not easy
523  // for that reason, a numerical table was put here.
524  Matrix1D<double> incompgamma;
525  incompgamma.initZeros(41);
526  VEC_ELEM(incompgamma,0) =0.0f;
527  VEC_ELEM(incompgamma,1) =0.00017212f;
528  VEC_ELEM(incompgamma,2) =0.0036598f;
529  VEC_ELEM(incompgamma,3) =0.018576f;
530  VEC_ELEM(incompgamma,4) =0.052653f;
531  VEC_ELEM(incompgamma,5) =0.10882f;
532  VEC_ELEM(incompgamma,6) =0.18474f;
533  VEC_ELEM(incompgamma,7) =0.27456f;
534  VEC_ELEM(incompgamma,8) =0.37116f;
535  VEC_ELEM(incompgamma,9) =0.4679f;
536  VEC_ELEM(incompgamma,10) =0.55951f;
537  VEC_ELEM(incompgamma,11) =0.64248f;
538  VEC_ELEM(incompgamma,12) =0.71494f;
539  VEC_ELEM(incompgamma,13) =0.77633f;
540  VEC_ELEM(incompgamma,14) =0.82701f;
541  VEC_ELEM(incompgamma,15) =0.86794f;
542  VEC_ELEM(incompgamma,16) =0.90037f;
543  VEC_ELEM(incompgamma,17) =0.92564f;
544  VEC_ELEM(incompgamma,18) =0.94504f;
545  VEC_ELEM(incompgamma,19) =0.95974f;
546  VEC_ELEM(incompgamma,20) =0.97075f;
547  VEC_ELEM(incompgamma,21) =0.97891f;
548  VEC_ELEM(incompgamma,22) =0.9849f;
549  VEC_ELEM(incompgamma,23) =0.98925f;
550  VEC_ELEM(incompgamma,24) =0.9924f;
551  VEC_ELEM(incompgamma,25) =0.99465f;
552  VEC_ELEM(incompgamma,26) =0.99626f;
553  VEC_ELEM(incompgamma,27) =0.9974f;
554  VEC_ELEM(incompgamma,28) =0.99819f;
555  VEC_ELEM(incompgamma,29) =0.99875f;
556  VEC_ELEM(incompgamma,30) =0.99914f;
557  VEC_ELEM(incompgamma,31) =0.99941f;
558  VEC_ELEM(incompgamma,32) =0.9996f;
559  VEC_ELEM(incompgamma,33) =0.99973f;
560  VEC_ELEM(incompgamma,34) =0.99982f;
561  VEC_ELEM(incompgamma,35) =0.99988f;
562  VEC_ELEM(incompgamma,36) =0.99992f;
563  VEC_ELEM(incompgamma,37) =0.99994f;
564  VEC_ELEM(incompgamma,38) =0.99996f;
565  VEC_ELEM(incompgamma,39) =0.99997f;
566  VEC_ELEM(incompgamma,40) =0.99998f;
567  val = VEC_ELEM(incompgamma,idx);
568  return val;
569 }
570 
571 
572 void ProgFSO::generateDirections(Matrix2D<float> &angles, bool alot)
573 {
574  if (alot == true)
575  {
576  angles.initZeros(2,321);
577  MAT_ELEM(angles,0,0)=0.0f; MAT_ELEM(angles,1,0)=0.0f;
578  MAT_ELEM(angles,0,1)=324.0f; MAT_ELEM(angles,1,1)=63.4349f;
579  MAT_ELEM(angles,0,2)=36.0f; MAT_ELEM(angles,1,2)=63.4349f;
580  MAT_ELEM(angles,0,3)=180.0f; MAT_ELEM(angles,1,3)=63.435f;
581  MAT_ELEM(angles,0,4)=252.0f; MAT_ELEM(angles,1,4)=63.435f;
582  MAT_ELEM(angles,0,5)=108.0f; MAT_ELEM(angles,1,5)=63.435f;
583  MAT_ELEM(angles,0,6)=324.0f; MAT_ELEM(angles,1,6)=31.7175f;
584  MAT_ELEM(angles,0,7)=36.0f; MAT_ELEM(angles,1,7)=31.7175f;
585  MAT_ELEM(angles,0,8)=0.0f; MAT_ELEM(angles,1,8)=58.2825f;
586  MAT_ELEM(angles,0,9)=288.0f; MAT_ELEM(angles,1,9)=58.2825f;
587  MAT_ELEM(angles,0,10)=342.0f; MAT_ELEM(angles,1,10)=90.0f;
588  MAT_ELEM(angles,0,11)=306.0f; MAT_ELEM(angles,1,11)=90.0f;
589  MAT_ELEM(angles,0,12)=72.0f; MAT_ELEM(angles,1,12)=58.2825f;
590  MAT_ELEM(angles,0,13)=18.0f; MAT_ELEM(angles,1,13)=90.0f;
591  MAT_ELEM(angles,0,14)=54.0f; MAT_ELEM(angles,1,14)=90.0f;
592  MAT_ELEM(angles,0,15)=90.0f; MAT_ELEM(angles,1,15)=90.0f;
593  MAT_ELEM(angles,0,16)=216.0f; MAT_ELEM(angles,1,16)=58.282f;
594  MAT_ELEM(angles,0,17)=144.0f; MAT_ELEM(angles,1,17)=58.282f;
595  MAT_ELEM(angles,0,18)=180.0f; MAT_ELEM(angles,1,18)=31.718f;
596  MAT_ELEM(angles,0,19)=252.0f; MAT_ELEM(angles,1,19)=31.718f;
597  MAT_ELEM(angles,0,20)=108.0f; MAT_ELEM(angles,1,20)=31.718f;
598  MAT_ELEM(angles,0,21)=346.3862f; MAT_ELEM(angles,1,21)=43.6469f;
599  MAT_ELEM(angles,0,22)=58.3862f; MAT_ELEM(angles,1,22)=43.6469f;
600  MAT_ELEM(angles,0,23)=274.3862f; MAT_ELEM(angles,1,23)=43.6469f;
601  MAT_ELEM(angles,0,24)=0.0f; MAT_ELEM(angles,1,24)=90.0f;
602  MAT_ELEM(angles,0,25)=72.0f; MAT_ELEM(angles,1,25)=90.0f;
603  MAT_ELEM(angles,0,26)=288.0f; MAT_ELEM(angles,1,26)=90.0f;
604  MAT_ELEM(angles,0,27)=225.7323f; MAT_ELEM(angles,1,27)=73.955f;
605  MAT_ELEM(angles,0,28)=153.7323f; MAT_ELEM(angles,1,28)=73.955f;
606  MAT_ELEM(angles,0,29)=216.0f; MAT_ELEM(angles,1,29)=26.565f;
607  MAT_ELEM(angles,0,30)=144.0f; MAT_ELEM(angles,1,30)=26.565f;
608  MAT_ELEM(angles,0,31)=0.0f; MAT_ELEM(angles,1,31)=26.5651f;
609  MAT_ELEM(angles,0,32)=72.0f; MAT_ELEM(angles,1,32)=26.5651f;
610  MAT_ELEM(angles,0,33)=288.0f; MAT_ELEM(angles,1,33)=26.5651f;
611  MAT_ELEM(angles,0,34)=350.2677f; MAT_ELEM(angles,1,34)=73.9549f;
612  MAT_ELEM(angles,0,35)=62.2677f; MAT_ELEM(angles,1,35)=73.9549f;
613  MAT_ELEM(angles,0,36)=278.2677f; MAT_ELEM(angles,1,36)=73.9549f;
614  MAT_ELEM(angles,0,37)=206.2677f; MAT_ELEM(angles,1,37)=73.955f;
615  MAT_ELEM(angles,0,38)=134.2677f; MAT_ELEM(angles,1,38)=73.955f;
616  MAT_ELEM(angles,0,39)=202.3862f; MAT_ELEM(angles,1,39)=43.647f;
617  MAT_ELEM(angles,0,40)=130.3862f; MAT_ELEM(angles,1,40)=43.647f;
618  MAT_ELEM(angles,0,41)=13.6138f; MAT_ELEM(angles,1,41)=43.6469f;
619  MAT_ELEM(angles,0,42)=85.6138f; MAT_ELEM(angles,1,42)=43.6469f;
620  MAT_ELEM(angles,0,43)=301.6138f; MAT_ELEM(angles,1,43)=43.6469f;
621  MAT_ELEM(angles,0,44)=9.7323f; MAT_ELEM(angles,1,44)=73.9549f;
622  MAT_ELEM(angles,0,45)=81.7323f; MAT_ELEM(angles,1,45)=73.9549f;
623  MAT_ELEM(angles,0,46)=297.7323f; MAT_ELEM(angles,1,46)=73.9549f;
624  MAT_ELEM(angles,0,47)=36.0f; MAT_ELEM(angles,1,47)=90.0f;
625  MAT_ELEM(angles,0,48)=324.0f; MAT_ELEM(angles,1,48)=90.0f;
626  MAT_ELEM(angles,0,49)=229.6138f; MAT_ELEM(angles,1,49)=43.647f;
627  MAT_ELEM(angles,0,50)=157.6138f; MAT_ELEM(angles,1,50)=43.647f;
628  MAT_ELEM(angles,0,51)=324.0f; MAT_ELEM(angles,1,51)=15.8587f;
629  MAT_ELEM(angles,0,52)=36.0f; MAT_ELEM(angles,1,52)=15.8587f;
630  MAT_ELEM(angles,0,53)=341.533f; MAT_ELEM(angles,1,53)=59.6208f;
631  MAT_ELEM(angles,0,54)=306.467f; MAT_ELEM(angles,1,54)=59.6208f;
632  MAT_ELEM(angles,0,55)=333.5057f; MAT_ELEM(angles,1,55)=76.5584f;
633  MAT_ELEM(angles,0,56)=314.4943f; MAT_ELEM(angles,1,56)=76.5584f;
634  MAT_ELEM(angles,0,57)=53.533f; MAT_ELEM(angles,1,57)=59.6208f;
635  MAT_ELEM(angles,0,58)=26.4943f; MAT_ELEM(angles,1,58)=76.5584f;
636  MAT_ELEM(angles,0,59)=45.5057f; MAT_ELEM(angles,1,59)=76.5584f;
637  MAT_ELEM(angles,0,60)=197.533f; MAT_ELEM(angles,1,60)=59.621f;
638  MAT_ELEM(angles,0,61)=162.467f; MAT_ELEM(angles,1,61)=59.621f;
639  MAT_ELEM(angles,0,62)=180.0; MAT_ELEM(angles,1,62)=47.576;
640  MAT_ELEM(angles,0,63)=269.533f; MAT_ELEM(angles,1,63)=59.621;
641  MAT_ELEM(angles,0,64)=252.0f; MAT_ELEM(angles,1,64)=47.576f;
642  MAT_ELEM(angles,0,65)=108.0f; MAT_ELEM(angles,1,65)=47.576f;
643  MAT_ELEM(angles,0,66)=324.0f; MAT_ELEM(angles,1,66)=47.5762f;
644  MAT_ELEM(angles,0,67)=36.0f; MAT_ELEM(angles,1,67)=47.5762f;
645  MAT_ELEM(angles,0,68)=18.467f; MAT_ELEM(angles,1,68)=59.6208f;
646  MAT_ELEM(angles,0,69)=170.4943f; MAT_ELEM(angles,1,69)=76.558f;
647  MAT_ELEM(angles,0,70)=117.5057f; MAT_ELEM(angles,1,70)=76.558f;
648  MAT_ELEM(angles,0,71)=189.5057f; MAT_ELEM(angles,1,71)=76.558f;
649  MAT_ELEM(angles,0,72)=242.4943f; MAT_ELEM(angles,1,72)=76.558f;
650  MAT_ELEM(angles,0,73)=261.5057f; MAT_ELEM(angles,1,73)=76.558f;
651  MAT_ELEM(angles,0,74)=98.4943f; MAT_ELEM(angles,1,74)=76.558f;
652  MAT_ELEM(angles,0,75)=234.467f; MAT_ELEM(angles,1,75)=59.621f;
653  MAT_ELEM(angles,0,76)=125.533f; MAT_ELEM(angles,1,76)=59.621f;
654  MAT_ELEM(angles,0,77)=180.0f; MAT_ELEM(angles,1,77)=15.859f;
655  MAT_ELEM(angles,0,78)=252.0f; MAT_ELEM(angles,1,78)=15.859f;
656  MAT_ELEM(angles,0,79)=90.467f; MAT_ELEM(angles,1,79)=59.621f;
657  MAT_ELEM(angles,0,80)=108.0f; MAT_ELEM(angles,1,80)=15.859f;
658  MAT_ELEM(angles,0,81)=0.0f; MAT_ELEM(angles,1,81)=42.8321f;
659  MAT_ELEM(angles,0,82)=72.0f; MAT_ELEM(angles,1,82)=42.8321f;
660  MAT_ELEM(angles,0,83)=288.0f; MAT_ELEM(angles,1,83)=42.8321f;
661  MAT_ELEM(angles,0,84)=4.7693f; MAT_ELEM(angles,1,84)=81.9488f;
662  MAT_ELEM(angles,0,85)=76.7693f; MAT_ELEM(angles,1,85)=81.9488f;
663  MAT_ELEM(angles,0,86)=292.7693f; MAT_ELEM(angles,1,86)=81.9488f;
664  MAT_ELEM(angles,0,87)=220.7693f; MAT_ELEM(angles,1,87)=81.9488f;
665  MAT_ELEM(angles,0,88)=148.7693f; MAT_ELEM(angles,1,88)=81.9488f;
666  MAT_ELEM(angles,0,89)=224.2677f; MAT_ELEM(angles,1,89)=34.924f;
667  MAT_ELEM(angles,0,90)=152.2677f; MAT_ELEM(angles,1,90)=34.924f;
668  MAT_ELEM(angles,0,91)=13.5146f; MAT_ELEM(angles,1,91)=20.3172f;
669  MAT_ELEM(angles,0,92)=85.5146f; MAT_ELEM(angles,1,92)=20.3172f;
670  MAT_ELEM(angles,0,93)=301.5146f; MAT_ELEM(angles,1,93)=20.3172f;
671  MAT_ELEM(angles,0,94)=346.1363f; MAT_ELEM(angles,1,94)=66.7276f;
672  MAT_ELEM(angles,0,95)=58.1363f; MAT_ELEM(angles,1,95)=66.7276f;
673  MAT_ELEM(angles,0,96)=274.1363f; MAT_ELEM(angles,1,96)=66.7276f;
674  MAT_ELEM(angles,0,97)=197.8362f; MAT_ELEM(angles,1,97)=75.105f;
675  MAT_ELEM(angles,0,98)=269.8362f; MAT_ELEM(angles,1,98)=75.105f;
676  MAT_ELEM(angles,0,99)=125.8362f; MAT_ELEM(angles,1,99)=75.105f;
677  MAT_ELEM(angles,0,100)=199.6899f; MAT_ELEM(angles,1,100)=51.609f;
678  MAT_ELEM(angles,0,101)=127.6899f; MAT_ELEM(angles,1,101)=51.609f;
679  MAT_ELEM(angles,0,102)=334.8124f; MAT_ELEM(angles,1,102)=45.0621f;
680  MAT_ELEM(angles,0,103)=46.8124f; MAT_ELEM(angles,1,103)=45.0621f;
681  MAT_ELEM(angles,0,104)=175.3133f; MAT_ELEM(angles,1,104)=83.2562f;
682  MAT_ELEM(angles,0,105)=247.3133f; MAT_ELEM(angles,1,105)=83.2562f;
683  MAT_ELEM(angles,0,106)=103.3133f; MAT_ELEM(angles,1,106)=83.2562f;
684  MAT_ELEM(angles,0,107)=229.8637f; MAT_ELEM(angles,1,107)=66.728f;
685  MAT_ELEM(angles,0,108)=157.8637f; MAT_ELEM(angles,1,108)=66.728f;
686  MAT_ELEM(angles,0,109)=202.4854f; MAT_ELEM(angles,1,109)=20.317f;
687  MAT_ELEM(angles,0,110)=130.4854f; MAT_ELEM(angles,1,110)=20.317f;
688  MAT_ELEM(angles,0,111)=16.3101f; MAT_ELEM(angles,1,111)=51.6091f;
689  MAT_ELEM(angles,0,112)=88.3101f; MAT_ELEM(angles,1,112)=51.6091f;
690  MAT_ELEM(angles,0,113)=304.3101f; MAT_ELEM(angles,1,113)=51.6091f;
691  MAT_ELEM(angles,0,114)=18.1638f; MAT_ELEM(angles,1,114)=75.1046f;
692  MAT_ELEM(angles,0,115)=306.1638f; MAT_ELEM(angles,1,115)=75.1046f;
693  MAT_ELEM(angles,0,116)=40.6867f; MAT_ELEM(angles,1,116)=83.2562f;
694  MAT_ELEM(angles,0,117)=328.6867f; MAT_ELEM(angles,1,117)=83.2562f;
695  MAT_ELEM(angles,0,118)=241.1876f; MAT_ELEM(angles,1,118)=45.062f;
696  MAT_ELEM(angles,0,119)=97.1876f; MAT_ELEM(angles,1,119)=45.062f;
697  MAT_ELEM(angles,0,120)=169.1876f; MAT_ELEM(angles,1,120)=45.062f;
698  MAT_ELEM(angles,0,121)=351.7323f; MAT_ELEM(angles,1,121)=34.9243f;
699  MAT_ELEM(angles,0,122)=63.7323f; MAT_ELEM(angles,1,122)=34.9243f;
700  MAT_ELEM(angles,0,123)=279.7323f; MAT_ELEM(angles,1,123)=34.9243f;
701  MAT_ELEM(angles,0,124)=355.2307f; MAT_ELEM(angles,1,124)=81.9488f;
702  MAT_ELEM(angles,0,125)=67.2307f; MAT_ELEM(angles,1,125)=81.9488f;
703  MAT_ELEM(angles,0,126)=283.2307f; MAT_ELEM(angles,1,126)=81.9488f;
704  MAT_ELEM(angles,0,127)=216.0f; MAT_ELEM(angles,1,127)=73.733f;
705  MAT_ELEM(angles,0,128)=144.0f; MAT_ELEM(angles,1,128)=73.733f;
706  MAT_ELEM(angles,0,129)=207.7323f; MAT_ELEM(angles,1,129)=34.924f;
707  MAT_ELEM(angles,0,130)=135.7323f; MAT_ELEM(angles,1,130)=34.924f;
708  MAT_ELEM(angles,0,131)=346.4854f; MAT_ELEM(angles,1,131)=20.3172f;
709  MAT_ELEM(angles,0,132)=58.4854f; MAT_ELEM(angles,1,132)=20.3172f;
710  MAT_ELEM(angles,0,133)=274.4854f; MAT_ELEM(angles,1,133)=20.3172f;
711  MAT_ELEM(angles,0,134)=341.8362f; MAT_ELEM(angles,1,134)=75.1046f;
712  MAT_ELEM(angles,0,135)=53.8362f; MAT_ELEM(angles,1,135)=75.1046f;
713  MAT_ELEM(angles,0,136)=202.1363f; MAT_ELEM(angles,1,136)=66.728f;
714  MAT_ELEM(angles,0,137)=130.1363f; MAT_ELEM(angles,1,137)=66.728f;
715  MAT_ELEM(angles,0,138)=190.8124f; MAT_ELEM(angles,1,138)=45.062f;
716  MAT_ELEM(angles,0,139)=262.8124f; MAT_ELEM(angles,1,139)=45.062f;
717  MAT_ELEM(angles,0,140)=118.8124f; MAT_ELEM(angles,1,140)=45.062f;
718  MAT_ELEM(angles,0,141)=343.6899f; MAT_ELEM(angles,1,141)=51.6091f;
719  MAT_ELEM(angles,0,142)=55.6899f; MAT_ELEM(angles,1,142)=51.6091f;
720  MAT_ELEM(angles,0,143)=271.6899f; MAT_ELEM(angles,1,143)=51.6091f;
721  MAT_ELEM(angles,0,144)=184.6867f; MAT_ELEM(angles,1,144)=83.2562f;
722  MAT_ELEM(angles,0,145)=256.6867f; MAT_ELEM(angles,1,145)=83.2562f;
723  MAT_ELEM(angles,0,146)=112.6867f; MAT_ELEM(angles,1,146)=83.2562f;
724  MAT_ELEM(angles,0,147)=234.1638f; MAT_ELEM(angles,1,147)=75.105f;
725  MAT_ELEM(angles,0,148)=90.1638f; MAT_ELEM(angles,1,148)=75.105f;
726  MAT_ELEM(angles,0,149)=162.1638f; MAT_ELEM(angles,1,149)=75.105f;
727  MAT_ELEM(angles,0,150)=229.5146f; MAT_ELEM(angles,1,150)=20.317f;
728  MAT_ELEM(angles,0,151)=157.5146f; MAT_ELEM(angles,1,151)=20.317f;
729  MAT_ELEM(angles,0,152)=25.1876f; MAT_ELEM(angles,1,152)=45.0621f;
730  MAT_ELEM(angles,0,153)=313.1876f; MAT_ELEM(angles,1,153)=45.0621f;
731  MAT_ELEM(angles,0,154)=13.8637f; MAT_ELEM(angles,1,154)=66.7276f;
732  MAT_ELEM(angles,0,155)=85.8637f; MAT_ELEM(angles,1,155)=66.7276f;
733  MAT_ELEM(angles,0,156)=301.8637f; MAT_ELEM(angles,1,156)=66.7276f;
734  MAT_ELEM(angles,0,157)=31.3133f; MAT_ELEM(angles,1,157)=83.2562f;
735  MAT_ELEM(angles,0,158)=319.3133f; MAT_ELEM(angles,1,158)=83.2562f;
736  MAT_ELEM(angles,0,159)=232.3101f; MAT_ELEM(angles,1,159)=51.609f;
737  MAT_ELEM(angles,0,160)=160.3101f; MAT_ELEM(angles,1,160)=51.609f;
738  MAT_ELEM(angles,0,161)=8.2677f; MAT_ELEM(angles,1,161)=34.9243f;
739  MAT_ELEM(angles,0,162)=80.2677f; MAT_ELEM(angles,1,162)=34.9243f;
740  MAT_ELEM(angles,0,163)=296.2677f; MAT_ELEM(angles,1,163)=34.9243f;
741  MAT_ELEM(angles,0,164)=0.0f; MAT_ELEM(angles,1,164)=73.733f;
742  MAT_ELEM(angles,0,165)=72.0f; MAT_ELEM(angles,1,165)=73.733f;
743  MAT_ELEM(angles,0,166)=288.0f; MAT_ELEM(angles,1,166)=73.733f;
744  MAT_ELEM(angles,0,167)=211.2307f; MAT_ELEM(angles,1,167)=81.9488f;
745  MAT_ELEM(angles,0,168)=139.2307f; MAT_ELEM(angles,1,168)=81.9488f;
746  MAT_ELEM(angles,0,169)=216.0f; MAT_ELEM(angles,1,169)=42.832f;
747  MAT_ELEM(angles,0,170)=144.0f; MAT_ELEM(angles,1,170)=42.832f;
748  MAT_ELEM(angles,0,171)=0.0f; MAT_ELEM(angles,1,171)=12.9432f;
749  MAT_ELEM(angles,0,172)=72.0f; MAT_ELEM(angles,1,172)=12.9432f;
750  MAT_ELEM(angles,0,173)=288.0f; MAT_ELEM(angles,1,173)=12.9432f;
751  MAT_ELEM(angles,0,174)=337.2786f; MAT_ELEM(angles,1,174)=68.041f;
752  MAT_ELEM(angles,0,175)=49.2786f; MAT_ELEM(angles,1,175)=68.041f;
753  MAT_ELEM(angles,0,176)=193.2786f; MAT_ELEM(angles,1,176)=68.041f;
754  MAT_ELEM(angles,0,177)=265.2786f; MAT_ELEM(angles,1,177)=68.041f;
755  MAT_ELEM(angles,0,178)=121.2786f; MAT_ELEM(angles,1,178)=68.041f;
756  MAT_ELEM(angles,0,179)=189.4537f; MAT_ELEM(angles,1,179)=53.278f;
757  MAT_ELEM(angles,0,180)=261.4537f; MAT_ELEM(angles,1,180)=53.278f;
758  MAT_ELEM(angles,0,181)=117.4537f; MAT_ELEM(angles,1,181)=53.278f;
759  MAT_ELEM(angles,0,182)=333.4537f; MAT_ELEM(angles,1,182)=53.2783f;
760  MAT_ELEM(angles,0,183)=45.4537f; MAT_ELEM(angles,1,183)=53.2783f;
761  MAT_ELEM(angles,0,184)=180.0f; MAT_ELEM(angles,1,184)=76.378f;
762  MAT_ELEM(angles,0,185)=252.0f; MAT_ELEM(angles,1,185)=76.378f;
763  MAT_ELEM(angles,0,186)=108.0f; MAT_ELEM(angles,1,186)=76.378f;
764  MAT_ELEM(angles,0,187)=238.7214f; MAT_ELEM(angles,1,187)=68.041f;
765  MAT_ELEM(angles,0,188)=94.7214f; MAT_ELEM(angles,1,188)=68.041f;
766  MAT_ELEM(angles,0,189)=166.7214f; MAT_ELEM(angles,1,189)=68.041f;
767  MAT_ELEM(angles,0,190)=216.0f; MAT_ELEM(angles,1,190)=12.943f;
768  MAT_ELEM(angles,0,191)=144.0f; MAT_ELEM(angles,1,191)=12.943f;
769  MAT_ELEM(angles,0,192)=26.5463f; MAT_ELEM(angles,1,192)=53.2783f;
770  MAT_ELEM(angles,0,193)=314.5463f; MAT_ELEM(angles,1,193)=53.2783f;
771  MAT_ELEM(angles,0,194)=22.7214f; MAT_ELEM(angles,1,194)=68.041f;
772  MAT_ELEM(angles,0,195)=310.7214f; MAT_ELEM(angles,1,195)=68.041f;
773  MAT_ELEM(angles,0,196)=36.0f; MAT_ELEM(angles,1,196)=76.3782f;
774  MAT_ELEM(angles,0,197)=324.0f; MAT_ELEM(angles,1,197)=76.3782f;
775  MAT_ELEM(angles,0,198)=242.5463f; MAT_ELEM(angles,1,198)=53.278f;
776  MAT_ELEM(angles,0,199)=98.5463f; MAT_ELEM(angles,1,199)=53.278f;
777  MAT_ELEM(angles,0,200)=170.5463f; MAT_ELEM(angles,1,200)=53.278f;
778  MAT_ELEM(angles,0,201)=336.7264f; MAT_ELEM(angles,1,201)=37.1611f;
779  MAT_ELEM(angles,0,202)=48.7264f; MAT_ELEM(angles,1,202)=37.1611f;
780  MAT_ELEM(angles,0,203)=351.0f; MAT_ELEM(angles,1,203)=90.0f;
781  MAT_ELEM(angles,0,204)=63.0f; MAT_ELEM(angles,1,204)=90.0f;
782  MAT_ELEM(angles,0,205)=279.0f; MAT_ELEM(angles,1,205)=90.0f;
783  MAT_ELEM(angles,0,206)=221.1634f; MAT_ELEM(angles,1,206)=66.042f;
784  MAT_ELEM(angles,0,207)=149.1634f; MAT_ELEM(angles,1,207)=66.042f;
785  MAT_ELEM(angles,0,208)=196.498f; MAT_ELEM(angles,1,208)=27.943f;
786  MAT_ELEM(angles,0,209)=268.498f; MAT_ELEM(angles,1,209)=27.943f;
787  MAT_ELEM(angles,0,210)=124.498f; MAT_ELEM(angles,1,210)=27.943f;
788  MAT_ELEM(angles,0,211)=340.498f; MAT_ELEM(angles,1,211)=27.9429f;
789  MAT_ELEM(angles,0,212)=52.498f; MAT_ELEM(angles,1,212)=27.9429f;
790  MAT_ELEM(angles,0,213)=346.0516f; MAT_ELEM(angles,1,213)=81.9568f;
791  MAT_ELEM(angles,0,214)=58.0516f; MAT_ELEM(angles,1,214)=81.9568f;
792  MAT_ELEM(angles,0,215)=274.0516f; MAT_ELEM(angles,1,215)=81.9568f;
793  MAT_ELEM(angles,0,216)=210.8366f; MAT_ELEM(angles,1,216)=66.042f;
794  MAT_ELEM(angles,0,217)=138.8366f; MAT_ELEM(angles,1,217)=66.042f;
795  MAT_ELEM(angles,0,218)=192.7264f; MAT_ELEM(angles,1,218)=37.161f;
796  MAT_ELEM(angles,0,219)=264.7264f; MAT_ELEM(angles,1,219)=37.161f;
797  MAT_ELEM(angles,0,220)=120.7264f; MAT_ELEM(angles,1,220)=37.161f;
798  MAT_ELEM(angles,0,221)=6.0948f; MAT_ELEM(angles,1,221)=50.7685f;
799  MAT_ELEM(angles,0,222)=78.0948f; MAT_ELEM(angles,1,222)=50.7685f;
800  MAT_ELEM(angles,0,223)=294.0948f; MAT_ELEM(angles,1,223)=50.7685f;
801  MAT_ELEM(angles,0,224)=13.9484f; MAT_ELEM(angles,1,224)=81.9568f;
802  MAT_ELEM(angles,0,225)=85.9484f; MAT_ELEM(angles,1,225)=81.9568f;
803  MAT_ELEM(angles,0,226)=301.9484f; MAT_ELEM(angles,1,226)=81.9568f;
804  MAT_ELEM(angles,0,227)=45.0f; MAT_ELEM(angles,1,227)=90.0f;
805  MAT_ELEM(angles,0,228)=333.0f; MAT_ELEM(angles,1,228)=90.0f;
806  MAT_ELEM(angles,0,229)=239.2736f; MAT_ELEM(angles,1,229)=37.161f;
807  MAT_ELEM(angles,0,230)=95.2736f; MAT_ELEM(angles,1,230)=37.161f;
808  MAT_ELEM(angles,0,231)=167.2736f; MAT_ELEM(angles,1,231)=37.161f;
809  MAT_ELEM(angles,0,232)=324.0f; MAT_ELEM(angles,1,232)=7.9294f;
810  MAT_ELEM(angles,0,233)=36.0f; MAT_ELEM(angles,1,233)=7.9294f;
811  MAT_ELEM(angles,0,234)=332.6069f; MAT_ELEM(angles,1,234)=61.2449f;
812  MAT_ELEM(angles,0,235)=315.3931f; MAT_ELEM(angles,1,235)=61.2449f;
813  MAT_ELEM(angles,0,236)=328.9523f; MAT_ELEM(angles,1,236)=69.9333f;
814  MAT_ELEM(angles,0,237)=319.0477f; MAT_ELEM(angles,1,237)=69.9333f;
815  MAT_ELEM(angles,0,238)=44.6069f; MAT_ELEM(angles,1,238)=61.2449f;
816  MAT_ELEM(angles,0,239)=31.0477f; MAT_ELEM(angles,1,239)=69.9333f;
817  MAT_ELEM(angles,0,240)=40.9523f; MAT_ELEM(angles,1,240)=69.9333f;
818  MAT_ELEM(angles,0,241)=188.6069f; MAT_ELEM(angles,1,241)=61.245f;
819  MAT_ELEM(angles,0,242)=171.3931f; MAT_ELEM(angles,1,242)=61.245f;
820  MAT_ELEM(angles,0,243)=180.0f; MAT_ELEM(angles,1,243)=55.506f;
821  MAT_ELEM(angles,0,244)=260.6069f; MAT_ELEM(angles,1,244)=61.245f;
822  MAT_ELEM(angles,0,245)=252.0f; MAT_ELEM(angles,1,245)=55.506f;
823  MAT_ELEM(angles,0,246)=108.0f; MAT_ELEM(angles,1,246)=55.506f;
824  MAT_ELEM(angles,0,247)=324.0f; MAT_ELEM(angles,1,247)=39.6468f;
825  MAT_ELEM(angles,0,248)=36.0f; MAT_ELEM(angles,1,248)=39.6468f;
826  MAT_ELEM(angles,0,249)=9.299f; MAT_ELEM(angles,1,249)=58.6205f;
827  MAT_ELEM(angles,0,250)=278.701f; MAT_ELEM(angles,1,250)=58.6205f;
828  MAT_ELEM(angles,0,251)=166.1881f; MAT_ELEM(angles,1,251)=83.2609f;
829  MAT_ELEM(angles,0,252)=121.8119f; MAT_ELEM(angles,1,252)=83.2609f;
830  MAT_ELEM(angles,0,253)=81.299f; MAT_ELEM(angles,1,253)=58.6205f;
831  MAT_ELEM(angles,0,254)=193.8119f; MAT_ELEM(angles,1,254)=83.2609f;
832  MAT_ELEM(angles,0,255)=238.1881f; MAT_ELEM(angles,1,255)=83.2609f;
833  MAT_ELEM(angles,0,256)=265.8119f; MAT_ELEM(angles,1,256)=83.2609f;
834  MAT_ELEM(angles,0,257)=94.1881f; MAT_ELEM(angles,1,257)=83.2609f;
835  MAT_ELEM(angles,0,258)=225.299f; MAT_ELEM(angles,1,258)=58.621f;
836  MAT_ELEM(angles,0,259)=134.701f; MAT_ELEM(angles,1,259)=58.621f;
837  MAT_ELEM(angles,0,260)=180.0f; MAT_ELEM(angles,1,260)=23.788f;
838  MAT_ELEM(angles,0,261)=252.0f; MAT_ELEM(angles,1,261)=23.788f;
839  MAT_ELEM(angles,0,262)=108.0f; MAT_ELEM(angles,1,262)=23.788f;
840  MAT_ELEM(angles,0,263)=353.9052f; MAT_ELEM(angles,1,263)=50.7685f;
841  MAT_ELEM(angles,0,264)=65.9052f; MAT_ELEM(angles,1,264)=50.7685f;
842  MAT_ELEM(angles,0,265)=281.9052f; MAT_ELEM(angles,1,265)=50.7685f;
843  MAT_ELEM(angles,0,266)=9.0f; MAT_ELEM(angles,1,266)=90.0f;
844  MAT_ELEM(angles,0,267)=81.0f; MAT_ELEM(angles,1,267)=90.0f;
845  MAT_ELEM(angles,0,268)=297.0f; MAT_ELEM(angles,1,268)=90.0f;
846  MAT_ELEM(angles,0,269)=229.9484f; MAT_ELEM(angles,1,269)=81.9568f;
847  MAT_ELEM(angles,0,270)=157.9484f; MAT_ELEM(angles,1,270)=81.9568f;
848  MAT_ELEM(angles,0,271)=235.502f; MAT_ELEM(angles,1,271)=27.943f;
849  MAT_ELEM(angles,0,272)=91.502f; MAT_ELEM(angles,1,272)=27.943f;
850  MAT_ELEM(angles,0,273)=163.502f; MAT_ELEM(angles,1,273)=27.943f;
851  MAT_ELEM(angles,0,274)=19.502f; MAT_ELEM(angles,1,274)=27.9429f;
852  MAT_ELEM(angles,0,275)=307.502f; MAT_ELEM(angles,1,275)=27.9429f;
853  MAT_ELEM(angles,0,276)=354.8366f; MAT_ELEM(angles,1,276)=66.0423f;
854  MAT_ELEM(angles,0,277)=66.8366f; MAT_ELEM(angles,1,277)=66.0423f;
855  MAT_ELEM(angles,0,278)=282.8366f; MAT_ELEM(angles,1,278)=66.0423f;
856  MAT_ELEM(angles,0,279)=202.0516f; MAT_ELEM(angles,1,279)=81.9568f;
857  MAT_ELEM(angles,0,280)=130.0516f; MAT_ELEM(angles,1,280)=81.9568f;
858  MAT_ELEM(angles,0,281)=209.9052f; MAT_ELEM(angles,1,281)=50.768f;
859  MAT_ELEM(angles,0,282)=137.9052f; MAT_ELEM(angles,1,282)=50.768f;
860  MAT_ELEM(angles,0,283)=23.2736f; MAT_ELEM(angles,1,283)=37.1611f;
861  MAT_ELEM(angles,0,284)=311.2736f; MAT_ELEM(angles,1,284)=37.1611f;
862  MAT_ELEM(angles,0,285)=5.1634f; MAT_ELEM(angles,1,285)=66.0423f;
863  MAT_ELEM(angles,0,286)=77.1634f; MAT_ELEM(angles,1,286)=66.0423f;
864  MAT_ELEM(angles,0,287)=293.1634f; MAT_ELEM(angles,1,287)=66.0423f;
865  MAT_ELEM(angles,0,288)=27.0f; MAT_ELEM(angles,1,288)=90.0f;
866  MAT_ELEM(angles,0,289)=315.0f; MAT_ELEM(angles,1,289)=90.0f;
867  MAT_ELEM(angles,0,290)=222.0948f; MAT_ELEM(angles,1,290)=50.768f;
868  MAT_ELEM(angles,0,291)=150.0948f; MAT_ELEM(angles,1,291)=50.768f;
869  MAT_ELEM(angles,0,292)=324.0f; MAT_ELEM(angles,1,292)=23.7881f;
870  MAT_ELEM(angles,0,293)=36.0f; MAT_ELEM(angles,1,293)=23.7881f;
871  MAT_ELEM(angles,0,294)=350.701f; MAT_ELEM(angles,1,294)=58.6205f;
872  MAT_ELEM(angles,0,295)=297.299f; MAT_ELEM(angles,1,295)=58.6205f;
873  MAT_ELEM(angles,0,296)=337.8119f; MAT_ELEM(angles,1,296)=83.2609f;
874  MAT_ELEM(angles,0,297)=310.1881f; MAT_ELEM(angles,1,297)=83.2609f;
875  MAT_ELEM(angles,0,298)=62.701f; MAT_ELEM(angles,1,298)=58.6205f;
876  MAT_ELEM(angles,0,299)=22.1881f; MAT_ELEM(angles,1,299)=83.2609f;
877  MAT_ELEM(angles,0,300)=49.8119f; MAT_ELEM(angles,1,300)=83.2609f;
878  MAT_ELEM(angles,0,301)=206.701f; MAT_ELEM(angles,1,301)=58.621f;
879  MAT_ELEM(angles,0,302)=153.299f; MAT_ELEM(angles,1,302)=58.621f;
880  MAT_ELEM(angles,0,303)=180.0f; MAT_ELEM(angles,1,303)=39.647f;
881  MAT_ELEM(angles,0,304)=252.0f; MAT_ELEM(angles,1,304)=39.647f;
882  MAT_ELEM(angles,0,305)=108.0f; MAT_ELEM(angles,1,305)=39.647f;
883  MAT_ELEM(angles,0,306)=324.0f; MAT_ELEM(angles,1,306)=55.5056f;
884  MAT_ELEM(angles,0,307)=36.0f; MAT_ELEM(angles,1,307)=55.5056f;
885  MAT_ELEM(angles,0,308)=27.3931f; MAT_ELEM(angles,1,308)=61.2449f;
886  MAT_ELEM(angles,0,309)=175.0477f; MAT_ELEM(angles,1,309)=69.933f;
887  MAT_ELEM(angles,0,310)=112.9523f; MAT_ELEM(angles,1,310)=69.933f;
888  MAT_ELEM(angles,0,311)=184.9523f; MAT_ELEM(angles,1,311)=69.933f;
889  MAT_ELEM(angles,0,312)=247.0477f; MAT_ELEM(angles,1,312)=69.933f;
890  MAT_ELEM(angles,0,313)=256.9523f; MAT_ELEM(angles,1,313)=69.933f;
891  MAT_ELEM(angles,0,314)=103.0477f; MAT_ELEM(angles,1,314)=69.933f;
892  MAT_ELEM(angles,0,315)=243.3931f; MAT_ELEM(angles,1,315)=61.245f;
893  MAT_ELEM(angles,0,316)=116.6069f; MAT_ELEM(angles,1,316)=61.245f;
894  MAT_ELEM(angles,0,317)=180.0f; MAT_ELEM(angles,1,317)=7.929f;
895  MAT_ELEM(angles,0,318)=252.0f; MAT_ELEM(angles,1,318)=7.929f;
896  MAT_ELEM(angles,0,319)=99.3931f; MAT_ELEM(angles,1,319)=61.245f;
897  MAT_ELEM(angles,0,320)=108.0f; MAT_ELEM(angles,1,320)=7.929f;
898  }
899  else
900  {
901  angles.initZeros(2,81);
902  MAT_ELEM(angles, 0, 0) = 0.000000f; MAT_ELEM(angles, 1, 0) = 0.000000f;
903  MAT_ELEM(angles, 0, 1) = 36.000000f; MAT_ELEM(angles, 1, 1) = 15.858741f;
904  MAT_ELEM(angles, 0, 2) = 36.000000f; MAT_ELEM(angles, 1, 2) = 31.717482f;
905  MAT_ELEM(angles, 0, 3) = 36.000000f; MAT_ELEM(angles, 1, 3) = 47.576224f;
906  MAT_ELEM(angles, 0, 4) = 36.000000f; MAT_ELEM(angles, 1, 4) = 63.434965f;
907  MAT_ELEM(angles, 0, 5) = 62.494295f; MAT_ELEM(angles, 1, 5) = -76.558393f;
908  MAT_ELEM(angles, 0, 6) = 54.000000f; MAT_ELEM(angles, 1, 6) = 90.000000f;
909  MAT_ELEM(angles, 0, 7) = 45.505705f; MAT_ELEM(angles, 1, 7) = 76.558393f;
910  MAT_ELEM(angles, 0, 8) = 108.000000f; MAT_ELEM(angles, 1, 8) = 15.858741f;
911  MAT_ELEM(angles, 0, 9) = 108.000000f; MAT_ELEM(angles, 1, 9) = 31.717482f;
912  MAT_ELEM(angles, 0, 10) = 108.000000f; MAT_ELEM(angles, 1, 10) = 47.576224f;
913  MAT_ELEM(angles, 0, 11) = 108.000000f; MAT_ELEM(angles, 1, 11) = 63.434965f;
914  MAT_ELEM(angles, 0, 12) = 134.494295f; MAT_ELEM(angles, 1, 12) = -76.558393f;
915  MAT_ELEM(angles, 0, 13) = 126.000000f; MAT_ELEM(angles, 1, 13) = 90.000000f;
916  MAT_ELEM(angles, 0, 14) = 117.505705f; MAT_ELEM(angles, 1, 14) = 76.558393f;
917  MAT_ELEM(angles, 0, 15) = 144.000000f; MAT_ELEM(angles, 1, 15) = -15.858741f;
918  MAT_ELEM(angles, 0, 16) = 144.000000f; MAT_ELEM(angles, 1, 16) = -31.717482f;
919  MAT_ELEM(angles, 0, 17) = 144.000000f; MAT_ELEM(angles, 1, 17) = -47.576224f;
920  MAT_ELEM(angles, 0, 18) = 144.000000f; MAT_ELEM(angles, 1, 18) = -63.434965f;
921  MAT_ELEM(angles, 0, 19) = 170.494295f; MAT_ELEM(angles, 1, 19) = 76.558393f;
922  MAT_ELEM(angles, 0, 20) = 162.000000f; MAT_ELEM(angles, 1, 20) = 90.000000f;
923  MAT_ELEM(angles, 0, 21) = 153.505705f; MAT_ELEM(angles, 1, 21) = -76.558393f;
924  MAT_ELEM(angles, 0, 22) = 72.000000f; MAT_ELEM(angles, 1, 22) = -15.858741f;
925  MAT_ELEM(angles, 0, 23) = 72.000000f; MAT_ELEM(angles, 1, 23) = -31.717482f;
926  MAT_ELEM(angles, 0, 24) = 72.000000f; MAT_ELEM(angles, 1, 24) = -47.576224f;
927  MAT_ELEM(angles, 0, 25) = 72.000000f; MAT_ELEM(angles, 1, 25) = -63.434965f;
928  MAT_ELEM(angles, 0, 26) = 98.494295f; MAT_ELEM(angles, 1, 26) = 76.558393f;
929  MAT_ELEM(angles, 0, 27) = 90.000000f; MAT_ELEM(angles, 1, 27) = 90.000000f;
930  MAT_ELEM(angles, 0, 28) = 81.505705f; MAT_ELEM(angles, 1, 28) = -76.558393f;
931  MAT_ELEM(angles, 0, 29) = 0.000000f; MAT_ELEM(angles, 1, 29) = -15.858741f;
932  MAT_ELEM(angles, 0, 30) = 0.000000f; MAT_ELEM(angles, 1, 30) = -31.717482f;
933  MAT_ELEM(angles, 0, 31) = 0.000000f; MAT_ELEM(angles, 1, 31) = -47.576224f;
934  MAT_ELEM(angles, 0, 32) = 0.000000f; MAT_ELEM(angles, 1, 32) = -63.434965f;
935  MAT_ELEM(angles, 0, 33) = 26.494295f; MAT_ELEM(angles, 1, 33) = 76.558393f;
936  MAT_ELEM(angles, 0, 34) = 18.000000f; MAT_ELEM(angles, 1, 34) = 90.000000f;
937  MAT_ELEM(angles, 0, 35) = 9.505705f; MAT_ELEM(angles, 1, 35) = -76.558393f;
938  MAT_ELEM(angles, 0, 36) = 12.811021f; MAT_ELEM(angles, 1, 36) = 42.234673f;
939  MAT_ELEM(angles, 0, 37) = 18.466996f; MAT_ELEM(angles, 1, 37) = 59.620797f;
940  MAT_ELEM(angles, 0, 38) = 0.000000f; MAT_ELEM(angles, 1, 38) = 90.000000f;
941  MAT_ELEM(angles, 0, 39) = 8.867209f; MAT_ELEM(angles, 1, 39) = 75.219088f;
942  MAT_ELEM(angles, 0, 40) = 72.000000f; MAT_ELEM(angles, 1, 40) = 26.565058f;
943  MAT_ELEM(angles, 0, 41) = 59.188979f; MAT_ELEM(angles, 1, 41) = 42.234673f;
944  MAT_ELEM(angles, 0, 42) = 84.811021f; MAT_ELEM(angles, 1, 42) = 42.234673f;
945  MAT_ELEM(angles, 0, 43) = 53.533003f; MAT_ELEM(angles, 1, 43) = 59.620797f;
946  MAT_ELEM(angles, 0, 44) = 72.000000f; MAT_ELEM(angles, 1, 44) = 58.282544f;
947  MAT_ELEM(angles, 0, 45) = 90.466996f; MAT_ELEM(angles, 1, 45) = 59.620797f;
948  MAT_ELEM(angles, 0, 46) = 72.000000f; MAT_ELEM(angles, 1, 46) = 90.000000f;
949  MAT_ELEM(angles, 0, 47) = 63.132791f; MAT_ELEM(angles, 1, 47) = 75.219088f;
950  MAT_ELEM(angles, 0, 48) = 80.867209f; MAT_ELEM(angles, 1, 48) = 75.219088f;
951  MAT_ELEM(angles, 0, 49) = 144.000000f; MAT_ELEM(angles, 1, 49) = 26.565058f;
952  MAT_ELEM(angles, 0, 50) = 131.188979f; MAT_ELEM(angles, 1, 50) = 42.234673f;
953  MAT_ELEM(angles, 0, 51) = 156.811021f; MAT_ELEM(angles, 1, 51) = 42.234673f;
954  MAT_ELEM(angles, 0, 52) = 125.533003f; MAT_ELEM(angles, 1, 52) = 59.620797f;
955  MAT_ELEM(angles, 0, 53) = 144.000000f; MAT_ELEM(angles, 1, 53) = 58.282544f;
956  MAT_ELEM(angles, 0, 54) = 162.466996f; MAT_ELEM(angles, 1, 54) = 59.620797f;
957  MAT_ELEM(angles, 0, 55) = 144.000000f; MAT_ELEM(angles, 1, 55) = 90.000000f;
958  MAT_ELEM(angles, 0, 56) = 135.132791f; MAT_ELEM(angles, 1, 56) = 75.219088f;
959  MAT_ELEM(angles, 0, 57) = 152.867209f; MAT_ELEM(angles, 1, 57) = 75.219088f;
960  MAT_ELEM(angles, 0, 58) = 180.000000f; MAT_ELEM(angles, 1, 58) = -26.565058f;
961  MAT_ELEM(angles, 0, 59) = 167.188979f; MAT_ELEM(angles, 1, 59) = -42.234673f;
962  MAT_ELEM(angles, 0, 60) = 180.000000f; MAT_ELEM(angles, 1, 60) = -58.282544f;
963  MAT_ELEM(angles, 0, 61) = 161.533003f; MAT_ELEM(angles, 1, 61) = -59.620797f;
964  MAT_ELEM(angles, 0, 62) = 171.132791f; MAT_ELEM(angles, 1, 62) = -75.219088f;
965  MAT_ELEM(angles, 0, 63) = 108.000000f; MAT_ELEM(angles, 1, 63) = -26.565058f;
966  MAT_ELEM(angles, 0, 64) = 120.811021f; MAT_ELEM(angles, 1, 64) = -42.234673f;
967  MAT_ELEM(angles, 0, 65) = 95.188979f; MAT_ELEM(angles, 1, 65) = -42.234673f;
968  MAT_ELEM(angles, 0, 66) = 126.466996f; MAT_ELEM(angles, 1, 66) = -59.620797f;
969  MAT_ELEM(angles, 0, 67) = 108.000000f; MAT_ELEM(angles, 1, 67) = -58.282544f;
970  MAT_ELEM(angles, 0, 68) = 89.533003f; MAT_ELEM(angles, 1, 68) = -59.620797f;
971  MAT_ELEM(angles, 0, 69) = 108.000000f; MAT_ELEM(angles, 1, 69) = 90.000000f;
972  MAT_ELEM(angles, 0, 70) = 116.867209f; MAT_ELEM(angles, 1, 70) = -75.219088f;
973  MAT_ELEM(angles, 0, 71) = 99.132791f; MAT_ELEM(angles, 1, 71) = -75.219088f;
974  MAT_ELEM(angles, 0, 72) = 36.000000f; MAT_ELEM(angles, 1, 72) = -26.565058f;
975  MAT_ELEM(angles, 0, 73) = 48.811021f; MAT_ELEM(angles, 1, 73) = -42.234673f;
976  MAT_ELEM(angles, 0, 74) = 23.188979f; MAT_ELEM(angles, 1, 74) = -42.234673f;
977  MAT_ELEM(angles, 0, 75) = 54.466996f; MAT_ELEM(angles, 1, 75) = -59.620797f;
978  MAT_ELEM(angles, 0, 76) = 36.000000f; MAT_ELEM(angles, 1, 76) = -58.282544f;
979  MAT_ELEM(angles, 0, 77) = 17.533003f; MAT_ELEM(angles, 1, 77) = -59.620797f;
980  MAT_ELEM(angles, 0, 78) = 36.000000f; MAT_ELEM(angles, 1, 78) = 90.000000f;
981  MAT_ELEM(angles, 0, 79) = 44.867209f; MAT_ELEM(angles, 1, 79) = -75.219088f;
982  MAT_ELEM(angles, 0, 80) = 27.132791f; MAT_ELEM(angles, 1, 80) = -75.219088f;
983  }
984 
985  angles *= PI/180.0;
986 }
987 
988 
989 void ProgFSO::anistropyParameter(const MultidimArray<float> &FSC,
990  MultidimArray<float> & aniParam, double thrs)
991 {
992  for (size_t k = 0; k<aniParam.nzyxdim; k++)
993  {
994  if (DIRECT_MULTIDIM_ELEM(FSC, k) >= thrs)
995  {
996  DIRECT_MULTIDIM_ELEM(aniParam, k) += 1.0;
997  }
998  }
999 }
1000 
1001 
1002 void ProgFSO::prepareData(MultidimArray<double> &half1, MultidimArray<double> &half2)
1003 {
1004  std::cout << "Reading data..." << std::endl;
1005  Image<double> imgHalf1, imgHalf2;
1006  imgHalf1.read(fnhalf1);
1007  imgHalf2.read(fnhalf2);
1008 
1009  half1 = imgHalf1();
1010  half2 = imgHalf2();
1011 
1012  // Applying the mask
1013  if (fnmask!="")
1014  {
1015  Image<double> mask;
1016  MultidimArray<double> &pmask = mask();
1017  mask.read(fnmask);
1018  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(pmask)
1019  {
1020  auto valmask = (double) DIRECT_MULTIDIM_ELEM(pmask, n);
1021  DIRECT_MULTIDIM_ELEM(half1, n) = DIRECT_MULTIDIM_ELEM(half1, n) * valmask;
1022  DIRECT_MULTIDIM_ELEM(half2, n) = DIRECT_MULTIDIM_ELEM(half2, n) * valmask;
1023  }
1024  }
1025 
1026  half1.setXmippOrigin();
1027  half2.setXmippOrigin();
1028 
1029  // fftshift must by applied before computing the fft. This will be computed later
1030  CenterFFT(half1, true);
1031  CenterFFT(half2, true);
1032 }
1033 
1034 
1035 void ProgFSO::saveAnisotropyToMetadata(MetaDataVec &mdAnisotropy,
1036  const MultidimArray<double> &freq,
1037  const MultidimArray<float> &anisotropy, const MultidimArray<double> &isotropyMatrix)
1038 {
1039  MDRowVec row;
1040  FOR_ALL_ELEMENTS_IN_ARRAY1D(anisotropy)
1041  {
1042  if (i>0)
1043  {
1044  row.setValue(MDL_RESOLUTION_FREQ, dAi(freq, i));
1045  row.setValue(MDL_RESOLUTION_FSO, static_cast<double>(dAi(anisotropy, i)));
1046  row.setValue(MDL_RESOLUTION_FREQREAL, 1.0/dAi(freq, i));
1047 
1048  row.setValue(MDL_RESOLUTION_ANISOTROPY, incompleteGammaFunction(dAi(isotropyMatrix, i)));
1049  mdAnisotropy.addRow(row);
1050  }
1051  }
1052  mdAnisotropy.write(fnOut+"/fso.xmd");
1053 }
1054 
1055 
1056 void ProgFSO::directionalFilter(MultidimArray<std::complex<double>> &FThalf1,
1057  MultidimArray<double> &threeDfsc,
1058  MultidimArray<double> &filteredMap, int m1sizeX, int m1sizeY, int m1sizeZ)
1059  {
1060  Image<double> imgHalf1, imgHalf2;
1061  imgHalf1.read(fnhalf1);
1062  imgHalf2.read(fnhalf2);
1063 
1064  auto &half1 = imgHalf1();
1065  auto &half2 = imgHalf2();
1066 
1067  FourierTransformer transformer1(FFTW_BACKWARD);
1068  transformer1.FourierTransform(half1, FThalf1);//, false);
1069  FourierTransformer transformer2(FFTW_BACKWARD);
1070  MultidimArray<std::complex<double>> FThalf2;
1071  FThalf2.resizeNoCopy(FThalf1);
1072  transformer1.FourierTransform(half2, FThalf2, false);
1073 
1074  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(threeDfsc)
1075  {
1076  DIRECT_MULTIDIM_ELEM(FThalf1, n) += DIRECT_MULTIDIM_ELEM(FThalf2, n);
1077  DIRECT_MULTIDIM_ELEM(FThalf1, n) *= DIRECT_MULTIDIM_ELEM(threeDfsc, n);
1078  }
1079  filteredMap.resizeNoCopy(m1sizeX, m1sizeY, m1sizeZ);
1080  transformer1.inverseFourierTransform(FThalf1, filteredMap);
1081  }
1082 
1083 
1084 void ProgFSO::directionalFilterHalves(MultidimArray<std::complex<double>> &FThalf1,
1085  MultidimArray<double> &threeDfsc)
1086 {
1087  Image<double> imgHalf1, imgHalf2;
1088  imgHalf1.read(fnhalf1);
1089  imgHalf2.read(fnhalf2);
1090 
1091  auto &half1 = imgHalf1();
1092  auto &half2 = imgHalf2();
1093 
1094  if (fnmask!="")
1095  {
1096  Image<double> msk;
1097  msk.read(fnmask);
1098  auto &pmsk = msk();
1099  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(threeDfsc)
1100  {
1101  double mskValue = DIRECT_MULTIDIM_ELEM(pmsk, n);
1102  DIRECT_MULTIDIM_ELEM(half1, n) *= mskValue;
1103  DIRECT_MULTIDIM_ELEM(half2, n) *= mskValue;
1104  }
1105  }
1106 
1107 
1108  FourierTransformer transformer1(FFTW_BACKWARD);
1109  transformer1.FourierTransform(half1, FThalf1);//, false);
1110  FourierTransformer transformer2(FFTW_BACKWARD);
1111  MultidimArray<std::complex<double>> FThalf2;
1112  FThalf2.resizeNoCopy(FThalf1);
1113  transformer2.FourierTransform(half2, FThalf2);
1114 
1115  // std::random_device r;
1116  // std::mt19937 gen(r());
1117  // std::uniform_int_distribution<> dis(1, 9);
1118 
1119  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(threeDfsc)
1120  {
1121  double fscValue;
1122  fscValue = DIRECT_MULTIDIM_ELEM(threeDfsc, n);
1123  // if (fscValue < 0.1)
1124  // {
1125  // DIRECT_MULTIDIM_ELEM(FThalf1, n) *= 0.05;//(rand() % 19 - 9)*1e-38;
1126  // DIRECT_MULTIDIM_ELEM(FThalf2, n) *= 0.05;//(rand() % 19 - 9)*1e-38;
1127  // }
1128  // else
1129  // {
1130  DIRECT_MULTIDIM_ELEM(FThalf1, n) *= 1;//fscValue;
1131  DIRECT_MULTIDIM_ELEM(FThalf2, n) *= 1;//fscValue;
1132  }
1133 
1134  MultidimArray<double> filteredMap1, filteredMap2;
1135 
1136  filteredMap1.resizeNoCopy(xvoldim, yvoldim, zvoldim);
1137  filteredMap2.resizeNoCopy(xvoldim, yvoldim, zvoldim);
1138 
1139  transformer1.inverseFourierTransform(FThalf1, filteredMap1);
1140  transformer2.inverseFourierTransform(FThalf2, filteredMap2);
1141 
1142  // int N_smoothing = 5;
1143 
1144  // int siz_z = zvoldim*0.5;
1145  // int siz_y = yvoldim*0.5;
1146  // int siz_x = xvoldim*0.5;
1147 
1148 
1149  // int limit_distance_x = (siz_x-N_smoothing);
1150  // int limit_distance_y = (siz_y-N_smoothing);
1151  // int limit_distance_z = (siz_z-N_smoothing);
1152 
1153  // long n=0;
1154  // for(int k=0; k<zvoldim; ++k)
1155  // {
1156  // double uz = (k - siz_z);
1157  // for(int i=0; i<yvoldim; ++i)
1158  // {
1159  // double uy = (i - siz_y);
1160  // for(int j=0; j<xvoldim; ++j)
1161  // {
1162  // double ux = (j - siz_x);
1163 
1164  // if (abs(ux)>=limit_distance_x)
1165  // {
1166  // double val = 0.5*(1+cos(PI*(limit_distance_x - abs(ux))/N_smoothing));
1167  // DIRECT_MULTIDIM_ELEM(filteredMap1, n) *= val;
1168  // DIRECT_MULTIDIM_ELEM(filteredMap2, n) *= val;
1169  // }
1170  // if (abs(uy)>=limit_distance_y)
1171  // {
1172  // double val = 0.5*(1+cos(PI*(limit_distance_y - abs(uy))/N_smoothing));
1173  // DIRECT_MULTIDIM_ELEM(filteredMap1, n) *= val;
1174  // DIRECT_MULTIDIM_ELEM(filteredMap2, n) *= val;
1175  // }
1176  // if (abs(uz)>=limit_distance_z)
1177  // {
1178  // double val = 0.5*(1+cos(PI*(limit_distance_z - abs(uz))/N_smoothing));
1179  // DIRECT_MULTIDIM_ELEM(filteredMap1, n) *= val;
1180  // DIRECT_MULTIDIM_ELEM(filteredMap2, n) *= val;
1181  // }
1182  // ++n;
1183  // }
1184  // }
1185  // }
1186 
1187  Image<double> saveImg(filteredMap1);
1188  saveImg.write(fnOut+"/filteredHalfMap1.mrc");
1189  saveImg() = filteredMap2;
1190  saveImg.write(fnOut+"/filteredHalfMap2.mrc");
1191 
1192 }
1193 
1194 void ProgFSO::resolutionDistribution(MultidimArray<double> &resDirFSC, FileName &fn)
1195  {
1196  Matrix2D<int> anglesResolution;
1197  const size_t Nrot = 360;
1198  const size_t Ntilt = 91;
1199  size_t objIdOut;
1200 
1201  MetaDataVec mdOut;
1202  Matrix2D<double> w, wt;
1203  w.initZeros(Nrot, Ntilt);
1204  wt = w;
1205  const float cosAngle = cosf(ang_con);
1206  const float aux = 4.0/((cosAngle -1)*(cosAngle -1));
1207  // Directional resolution is store in a metadata
1208 
1209  for (int i=0; i<Nrot; i++)
1210  {
1211  float rotmatrix = i*PI/180.0;
1212  float cr = cosf(rotmatrix);
1213  float sr = sinf(rotmatrix);
1214 
1215  for (int j=0; j<Ntilt; j++)
1216  {
1217  float tiltmatrix = j*PI/180.0;
1218  // position on the spehere
1219  float st = sinf(tiltmatrix);
1220  float xx = st*cr;
1221  float yy = st*sr;
1222  float zz = cosf(tiltmatrix);
1223 
1224  // initializing the weights
1225  double w = 0;
1226  double wt = 0;
1227 
1228  for (size_t k = 0; k<angles.mdimx; k++)
1229  {
1230 
1231  float rot = MAT_ELEM(angles, 0, k);
1232  float tilt = MAT_ELEM(angles, 1, k);
1233 
1234  // position of the direction on the sphere
1235  float x_dir = sinf(tilt)*cosf(rot);
1236  float y_dir = sinf(tilt)*sinf(rot);
1237  float z_dir = cosf(tilt);
1238 
1239 
1240  float cosine = fabs(x_dir*xx + y_dir*yy + z_dir*zz);
1241  if (cosine>=cosAngle)
1242  {
1243  cosine = expf( -((cosine -1)*(cosine -1))*aux );
1244  w += cosine*( dAi(resDirFSC, k) );
1245  wt += cosine;
1246  }
1247  }
1248 
1249  double wRes = w/wt;
1250  {
1251  MDRowVec row;
1252  row.setValue(MDL_ANGLE_ROT, (double) i);
1253  row.setValue(MDL_ANGLE_TILT, (double) j);
1254  row.setValue(MDL_RESOLUTION_FRC, wRes);
1255  mdOut.addRow(row);
1256  }
1257  }
1258  }
1259 
1260  mdOut.write(fn);
1261  }
1262 
1263 
1264 void ProgFSO::getCompleteFourier(MultidimArray<double> &V, MultidimArray<double> &newV,
1265  int m1sizeX, int m1sizeY, int m1sizeZ)
1266  {
1267  newV.resizeNoCopy(m1sizeX, m1sizeY, m1sizeZ);
1268  int ndim=3;
1269  if (m1sizeX==1)
1270  {
1271  ndim=2;
1272  if (m1sizeY==1)
1273  ndim=1;
1274  }
1275  double *ptrSource=nullptr;
1276  double *ptrDest=nullptr;
1277  FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(newV)
1278  {
1279  ptrDest=(double*)&DIRECT_A3D_ELEM(newV,k,i,j);
1280  if (j<XSIZE(V))
1281  {
1282  ptrSource=(double*)&DIRECT_A3D_ELEM(V,k,i,j);
1283  *ptrDest=*ptrSource;
1284  }
1285  else
1286  {
1287  ptrSource=(double*)&DIRECT_A3D_ELEM(V,
1288  (m1sizeZ-k)%m1sizeZ,
1289  (m1sizeY-i)%m1sizeY,
1290  m1sizeX-j);
1291  *ptrDest=*ptrSource;
1292  }
1293  }
1294  }
1295 
1296 
1297 void ProgFSO::createFullFourier(MultidimArray<double> &fourierHalf, FileName &fnMap,
1298  int m1sizeX, int m1sizeY, int m1sizeZ)
1299  {
1300  MultidimArray<double> fullMap;
1301  getCompleteFourier(fourierHalf, fullMap, m1sizeX, m1sizeY, m1sizeZ);
1302  CenterFFT(fullMap, true);
1303  Image<double> saveImg;
1304  saveImg() = fullMap;
1305  saveImg.write(fnMap);
1306  }
1307 
1308 
1309 void ProgFSO::run()
1310  {
1311  std::cout << "Starting ... " << std::endl << std::endl;
1312 
1313  MultidimArray<double> half1, half2;
1314  MultidimArray<double> &phalf1 = half1, &phalf2 = half2;
1315 
1316  //This read the data and applies a fftshift to in the next step compute the fft
1317  prepareData(half1, half2);
1318 
1319  //Computing the FFT
1320  FourierTransformer transformer2(FFTW_BACKWARD), transformer1(FFTW_BACKWARD);
1321  transformer1.setThreadsNumber(Nthreads);
1322  transformer2.setThreadsNumber(Nthreads);
1323 
1324  transformer1.FourierTransform(half1, FT1, false);
1325  transformer2.FourierTransform(half2, FT2, false);
1326 
1327  // Defining frequencies freq_fourier_x,y,z and freqMap
1328  // The number of frequencies in each shell freqElem is determined
1329  defineFrequencies(FT1, phalf1);
1330 
1331  half1.clear();
1332  half2.clear();
1333 
1334  // Storing the shell of both maps as vectors global
1335  // The global FSC is also computed
1336  MultidimArray<double> freq;
1337  arrangeFSC_and_fscGlobal(sampling, thrs, freq);
1338 
1339  FT2.clear();
1340 
1341  // Generating the set of directions to be analyzed
1342  // And converting the cone angle to radians
1343  generateDirections(angles, true);
1344  ang_con = ang_con*PI/180.0;
1345 
1346  // Preparing the metadata for storing the FSO
1347  MultidimArray<double> resDirFSC(angles.mdimx);//, aniParam;
1348  //aniParam.initZeros(xvoldim/2+1);
1349 
1350  // Computing directional FSC and 3DFSC
1351  // Create one copy for each thread
1352  std::vector<MultidimArray<float>> threeD_FSCs(Nthreads);
1353  std::vector<MultidimArray<float>> normalizationMaps(Nthreads);
1354  std::vector<MultidimArray<float>> aniParams(Nthreads);
1355  for (auto & ma : threeD_FSCs) {
1356  ma.initZeros(real_z1z2);
1357  }
1358  for (auto & ma : normalizationMaps) {
1359  ma.initZeros(real_z1z2);
1360  }
1361  for (auto & ma : aniParams) {
1362  ma.initZeros(xvoldim/2+1);
1363  }
1364 
1365  // Binham Test
1366  // We will have as many vectors as threads
1367  // Each thread considers a vector with the matrices of the Bingham test
1368  // instantiate vector object of type vector<Matrix2D<float>>
1369  std::vector<std::vector<Matrix2D<double>>> isotropyMatrices;
1370 
1371  // resize the vector to `freq.nzyxdim` elements of type vector<float>, each having size `C`
1372  isotropyMatrices.resize(Nthreads, std::vector<Matrix2D<double>>(freq.nzyxdim));
1373 
1374  // Each component of the vector is a 0-matrix (null matrix)
1375  for (auto & im : isotropyMatrices)
1376  {
1377  for (auto & m : im)
1378  m.initZeros(3,3);
1379  }
1380 
1381  ctpl::thread_pool threadPool;
1382  threadPool.resize(Nthreads);
1383  auto futures = std::vector<std::future<void>>();
1384  futures.reserve(angles.mdimx);
1385 
1386  for (size_t k = 0; k<angles.mdimx; k++)
1387  {
1388  futures.emplace_back(threadPool.push(
1389  [k, this, &resDirFSC, &threeD_FSCs, &normalizationMaps, &aniParams, &isotropyMatrices](int thrId){
1390  float rot = MAT_ELEM(angles, 0, k);
1391  float tilt = MAT_ELEM(angles, 1, k);
1392 
1393  double resInterp = -1;
1394  MultidimArray<float> fsc;
1395 
1396  auto &threeD_FSC = threeD_FSCs.at(thrId);
1397  auto &normalizationMap = normalizationMaps.at(thrId);
1398  auto &freqMat = isotropyMatrices.at(thrId);
1399  // Estimating the direction FSC along the direction given by rot and tilt
1400  fscDir_fast(fsc, rot, tilt, threeD_FSC, normalizationMap, thrs, resInterp, freqMat, k);
1401 
1402  printf ("Direction %zu/%zu -> %.2f A \n", k, angles.mdimx, resInterp);
1403 
1404  dAi(resDirFSC, k) = resInterp;
1405 
1406  // Updating the FSO curve
1407  anistropyParameter(fsc, aniParams.at(thrId), thrs);
1408  }));
1409  }
1410 
1411  // wait till done
1412  for (auto &f : futures) {
1413  f.get();
1414  }
1415 
1416  std::cout << "----- Directional resolution estimated -----" << std::endl << std::endl;
1417  std::cout << "Preparing results ..." << std::endl;
1418 
1419  // merge results from multiple threads
1420  for (size_t i = 1; i < aniParams.size(); ++i) {
1421  aniParams.at(0) += aniParams.at(i);
1422  }
1423 
1424  // For the Bingham Test
1425  // merge results from multiple threads
1426  for (size_t i = 0; i < freq.nzyxdim; ++i)
1427  {
1428  for (size_t j = 1; j<Nthreads; ++j)
1429  {
1430  isotropyMatrices.at(0)[i] += isotropyMatrices.at(j)[i];
1431  }
1432  isotropyMatrices.at(0)[i] /= aniParams.at(0)[i];
1433  }
1434 
1435  //auto &freqMat = isotropyMatrices.at(0);
1436  auto &aniParam = aniParams.at(0);
1437  MultidimArray<double> isotropyMatrix;
1438  isotropyMatrix.initZeros(freq.nzyxdim);
1439 
1440  for (size_t i = 0; i< freq.nzyxdim; ++i)
1441  {
1442  double trT2 = 0;
1443  if (aniParams.at(0)[i] >0)
1444  {
1445  // Let us define T = 1/n*Sum(wi*xi*xi) => Tr(T^2) = x*x + y*y + z*z;
1446  //This is the Bingham Test (1/2)(p-1)*(p-2)*n*Sum(Tr(T^2) - 1/p)
1447  //std::cout << isotropyMatrices.at(0)[i] << aniParams.at(0)[i] << std::endl;
1448  Matrix2D<double> T2;
1449  Matrix2D<double> T;
1450  T = isotropyMatrices.at(0)[i];//freqMat[i]/DIRECT_MULTIDIM_ELEM(aniParam, i);
1451  T2 = T*T;
1452  trT2 = T2.trace();
1453  double pdim = 3;
1454  dAi(isotropyMatrix, i) = 0.5*pdim*(pdim+2)*(2*aniParams.at(0)[i])*(trT2 - 1./pdim);
1455  // The factor 2 is because we need to compute the point in the whole sphere, but currently
1456  // we are measuing half of the sphere due to the symmetry of the problem
1457  // S = 0.5 * p * (p+2) * n * (trace(T2) - 1/p); Where p is the
1458  // dimension of the sphere (p=3) and n the number of points to
1459  // determine if they are uniformly or not. T = x*x';
1460 
1461  // The hypohtesis test considers p = 0.05.
1462  // p=0.05; nu = 5; x = chi2inv(p,nu);
1463  }
1464  }
1465 
1466 
1467  // ANISOTROPY CURVE
1468  aniParams.at(0) /= (double) angles.mdimx;
1469  for (size_t k = 0; k<5; k++)
1470  DIRECT_MULTIDIM_ELEM(aniParams.at(0), k) = 1.0;
1471  MetaDataVec mdani;
1472  saveAnisotropyToMetadata(mdani, freq, aniParams.at(0), isotropyMatrix);
1473 
1474 
1475  // DIRECTIONAL RESOLUTION DISTRIBUTION ON THE PROJECTION SPHERE
1476  FileName fn;
1477  fn = fnOut+"/Resolution_Distribution.xmd";
1478 
1479  resolutionDistribution(resDirFSC, fn);
1480 
1481 
1482  // 3DFSC and ANISOTROPIC FILTER
1483  if (do_3dfsc_filter)
1484  {
1485  // HALF 3DFSC MAP
1486  MultidimArray<double> d3_FSCMap;
1487  MultidimArray<double> d3_aniFilter;
1488  d3_FSCMap.resizeNoCopy(FT1);
1489  d3_FSCMap.initConstant(0);
1490  d3_aniFilter.resizeNoCopy(FT1);
1491  d3_aniFilter.initConstant(0);
1492 
1493  // merge results from multiple threads
1494  for (size_t i = 1; i < Nthreads; ++i) {
1495  threeD_FSCs.at(0) += threeD_FSCs.at(i);
1496  normalizationMaps.at(0) += normalizationMaps.at(i);
1497  }
1498  auto &threeD_FSC = threeD_FSCs.at(0);
1499  auto &normalizationMap = normalizationMaps.at(0);
1500 
1501  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(threeD_FSC)
1502  {
1503  double value = DIRECT_MULTIDIM_ELEM(threeD_FSC, n) /= DIRECT_MULTIDIM_ELEM(normalizationMap, n);
1504  if (std::isnan(value))
1505  value = 1.0;
1506 
1507  size_t idx = DIRECT_MULTIDIM_ELEM(arr2indx, n);
1508  DIRECT_MULTIDIM_ELEM(d3_FSCMap, idx) = value;
1509 
1510  if (DIRECT_MULTIDIM_ELEM(threeD_FSC, n)> thrs)//&& (DIRECT_MULTIDIM_ELEM(aniFilter, n) <1))
1511  {
1512  DIRECT_MULTIDIM_ELEM(d3_aniFilter, idx) = 1; //;DIRECT_MULTIDIM_ELEM(aniFilter, n) = 1;
1513  }
1514  else
1515  {
1516  DIRECT_MULTIDIM_ELEM(d3_aniFilter, idx) = value;//DIRECT_MULTIDIM_ELEM(aniFilter, n);
1517  }
1518 
1519  }
1520 
1521  // This code fix the empty line line in Fourier space
1522  size_t auxVal;
1523  auxVal = YSIZE(d3_FSCMap)/2;
1524 
1525  size_t j = 0;
1526  for(size_t i=(auxVal+1); i<YSIZE(d3_FSCMap); ++i)
1527  {
1528  for(size_t k=0; k<ZSIZE(d3_FSCMap); ++k)
1529  {
1530  DIRECT_A3D_ELEM(d3_FSCMap,k,i,0) = DIRECT_A3D_ELEM(d3_FSCMap,k,i,1);
1531  DIRECT_A3D_ELEM(d3_aniFilter,k,i,0) = DIRECT_A3D_ELEM(d3_aniFilter,k,i,1);
1532  }
1533  }
1534  size_t i=0;
1535  for(size_t k=0; k<ZSIZE(d3_FSCMap); ++k)
1536  {
1537  DIRECT_A3D_ELEM(d3_FSCMap,k,0,0) = DIRECT_A3D_ELEM(d3_FSCMap,k,0, 1);
1538  DIRECT_A3D_ELEM(d3_aniFilter,k,0,0) = DIRECT_A3D_ELEM(d3_aniFilter,k,0,1);
1539  }
1540 
1541 
1542  Image<double> img;
1543  img() = d3_FSCMap;
1544  img.write("threeD_FSC.mrc");
1545 
1546  double sigma = 3;
1547 
1548  realGaussianFilter(d3_aniFilter, sigma);
1549 
1550  // DIRECTIONAL FILTERED MAP
1551  MultidimArray<double> filteredMap;
1552  //directionalFilter(FT1, d3_aniFilter, filteredMap, xvoldim, yvoldim, zvoldim);
1553  directionalFilterHalves(FT1, d3_aniFilter);
1554 
1555  //FULL 3DFSC MAP
1556  fn = fnOut+"/3dFSC.mrc";
1557  createFullFourier(d3_FSCMap, fn, xvoldim, yvoldim, zvoldim);
1558  }
1559 
1560 
1561 
1562  std::cout << "-------------Finished-------------" << std::endl;
1563 }
1564 
1565 
dAi
#define dAi(v, i)
Definition: multidim_array_base.h:530
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
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
MDL_RESOLUTION_FSO
Fourier shell occupancy (double)
Definition: metadata_label.h:391
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
FourierTransformer::inverseFourierTransform
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
FourierTransformer::setThreadsNumber
void setThreadsNumber(int tNumber)
Definition: xmipp_fftw.h:119
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
ctpl::thread_pool::push
auto push(F &&f, Rest &&... rest) -> std::future< decltype(f(0, rest...))>
Definition: ctpl.h:152
resolution_fso.h
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
MDL_RESOLUTION_ANISOTROPY
Resolution anisotropy used to store the significance of the Bingham Test (double) ...
Definition: metadata_label.h:383
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
MDRowVec::setValue
void setValue(const MDObject &object) override
Definition: metadata_row_vec.cpp:199
Matrix2D::trace
T trace() const
Definition: matrix2d.cpp:1304
ImageBase::write
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
Definition: xmipp_image_base.cpp:292
w
doublereal * w
Definition: numerical_recipes.cpp:2477
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
FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:322
ctpl::thread_pool::resize
void resize(int nThreads)
Definition: ctpl.h:70
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
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
ctpl_stl.h
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
MetaDataVec::setValue
bool setValue(const MDObject &mdValueIn, size_t id)
Definition: metadata_vec.cpp:282
MetaDataVec::addObject
size_t addObject() override
Definition: metadata_vec.cpp:435
MDL_RESOLUTION_FREQREAL
Frequency in A (double)
Definition: metadata_label.h:390
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
for
for(j=1;j<=i__1;++j)
Definition: numerical_recipes.cpp:2310
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
Ntilt
#define Ntilt
Definition: project.cpp:557
DIRECT_A3D_ELEM
#define DIRECT_A3D_ELEM(v, k, i, j)
Definition: multidim_array_base.h:239
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
MDL_RESOLUTION_FREQ
Frequency in 1/A (double)
Definition: metadata_label.h:388
m
int m
Definition: image_find_center.cpp:57
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
realGaussianFilter
void realGaussianFilter(MultidimArray< double > &img, double sigma)
Definition: fourier_filter.cpp:849
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
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
ctpl::thread_pool
Definition: ctpl.h:48
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
ProgFSO::readParams
void readParams()
Definition: resolution_fso.cpp:95
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
Nrot
#define Nrot
Definition: project.cpp:556
ProgFSO::defineParams
void defineParams()
Definition: resolution_fso.cpp:36
MDL_RESOLUTION_FRC
Fourier shell correlation (double)
Definition: metadata_label.h:386
Generated by   doxygen 1.8.13