Xmipp  v3.23.11-Nereus
resolution_directional.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Jose Luis Vilas, jlvilas@cnb.csic.es
4  * Carlos Oscar S. Sorzano coss@cnb.csic.es (2018)
5  *
6  * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
21  * 02111-1307 USA
22  *
23  * All comments concerning this program package may be sent to the
24  * e-mail address 'xmipp@cnb.csic.es'
25  ***************************************************************************/
26 
27 #include "resolution_directional.h"
28 //#define DEBUG
29 //#define DEBUG_MASK
30 //#define DEBUG_DIR
31 //#define DEBUG_FILTER
32 //#define MONO_AMPLITUDE
33 //define DEBUG_SYMMETRY
34 
35 void ProgResDir::readParams()
36 {
37  fnVol = getParam("--vol");
38  fnOut = getParam("-o");
39  fnMask = getParam("--mask");
40  sampling = getDoubleParam("--sampling_rate");
41  R = getDoubleParam("--volumeRadius");
42  significance = getDoubleParam("--significance");
43  res_step = getDoubleParam("--resStep");
44  fnradial = getParam("--radialRes");
45  fnazimuthal = getParam("--azimuthalRes");
46  fnRadialAvG = getParam("--radialAvG");
47  fnHighestResolution = getParam("--highestResolutionVol");
48  fnLowestResolution = getParam("--lowestResolutionVol");
49  fnDoa1 = getParam("--doa1");
50  fnDoa2 = getParam("--doa2");
51  fnMDThr = getParam("--radialAzimuthalThresholds");
52  fnMonoRes = getParam("--monores");
53  fnprefMin = getParam("--prefMin");
54  fnZscore = getParam("--zScoremap");
55  Nthr = getIntParam("--threads");
56  fastCompute = checkParam("--fast");
57 }
58 
59 
60 void ProgResDir::defineParams()
61 {
62  addUsageLine("This function determines the local resolution of a map");
63  addParamsLine(" --vol <vol_file=\"\"> : Input volume");
64  addParamsLine(" --mask <vol_file=\"\"> : Mask defining the macromolecule");
65  addParamsLine(" -o <output=\"MGresolution.vol\"> : Local resolution volume (in Angstroms)");
66  addParamsLine(" [--sampling_rate <s=1>] : Sampling rate (A/px)");
67  addParamsLine(" [--resStep <s=0.5>] : Resolution step (precision) in A");
68  addParamsLine(" [--volumeRadius <s=100>] : This parameter determines the radius of a sphere where the volume is");
69  addParamsLine(" [--significance <s=0.95>] : The level of confidence for the hypothesis test.");
70  addParamsLine(" --radialRes <vol_file=\"\"> : Output radial resolution map");
71  addParamsLine(" --azimuthalRes <vol_file=\"\"> : Output azimuthal resolution map");
72  addParamsLine(" --highestResolutionVol <vol_file=\"\"> : Output highest resolution map");
73  addParamsLine(" --lowestResolutionVol <vol_file=\"\"> : Output lowest resolution map");
74  addParamsLine(" --doa1 <vol_file=\"\"> : Output doa1 map");
75  addParamsLine(" --doa2 <vol_file=\"\"> : Output doa2 map");
76  addParamsLine(" --radialAzimuthalThresholds <vol_file=\"\"> : Radial and azimuthal threshold for representation resolution maps");
77  addParamsLine(" --radialAvG <md_file=\"\"> : Radial Average of higesht, lowest, azimuthal, and radial, and local resolution map");
78  addParamsLine(" --monores <vol_file=\"\"> : Local resolution map");
79  addParamsLine(" --prefMin <vol_file=\"\"> : Metadata of highest resolution per direction");
80  addParamsLine(" [--threads <s=4>] : Number of threads");
81  addParamsLine(" --zScoremap <vol_file=\"\"> : Local zScore map, voxel with zscore higher than 3 are weird");
82  addParamsLine(" [--fast] : Fast computation");
83 }
84 
85 void ProgResDir::produceSideInfo()
86 {
87 
88  std::cout << "Starting..." << std::endl;
89 
90  Image<double> V;
91  V.read(fnVol);
92 
93  V().setXmippOrigin();
94 
95  //Sweeping the projection sphere
96  std::cout << "Obtaining angular projections..." << std::endl;
97  generateGridProjectionMatching(angles);
98 
99  FourierTransformer transformer;
100 
101  MultidimArray<double> &inputVol = V();
102  VRiesz.resizeNoCopy(inputVol);
103  N_freq = ZSIZE(inputVol);
104  maxRes = ZSIZE(inputVol);
105  minRes = 2*sampling;
106 
107  transformer_inv.setThreadsNumber(Nthr);
108 
109  transformer.FourierTransform(inputVol, fftV);
110  iu.initZeros(fftV);
111 
112  // Frequency volume
113  double uz, uy, ux, uz2, u2, uz2y2;
114  long n=0;
115 
116  for(size_t k=0; k<ZSIZE(fftV); ++k)
117  {
118  FFT_IDX2DIGFREQ(k,ZSIZE(inputVol),uz);
119  uz2=uz*uz;
120 
121  for(size_t i=0; i<YSIZE(fftV); ++i)
122  {
123  FFT_IDX2DIGFREQ(i,YSIZE(inputVol),uy);
124  uz2y2=uz2+uy*uy;
125 
126  for(size_t j=0; j<XSIZE(fftV); ++j)
127  {
128  FFT_IDX2DIGFREQ(j,XSIZE(inputVol), ux);
129  u2=uz2y2+ux*ux;
130  if ((k != 0) || (i != 0) || (j != 0))
131  DIRECT_MULTIDIM_ELEM(iu,n) = 1.0/sqrt(u2);
132  else
133  DIRECT_MULTIDIM_ELEM(iu,n) = 1e38;
134  ++n;
135  }
136  }
137  }
138 
139  // Prepare mask
140  MultidimArray<int> &pMask=mask();
141 
142  if (fnMask != "")
143  {
144  mask.read(fnMask);
145  mask().setXmippOrigin();
146  }
147  else
148  {
149  std::cout << "Error: a mask ought to be provided" << std::endl;
150  exit(0);
151  }
152 
153  N_smoothing = 7;
154  NVoxelsOriginalMask = 0;
155  double radius = 0;
156  FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
157  {
158  if (A3D_ELEM(pMask, k, i, j) == 1)
159  {
160  if ((k*k + i*i + j*j)>radius)
161  radius = k*k + i*i + j*j;
162  }
163 // std::cout << "i j k " << i << " " << j << " " << k << std::endl;
164 
165  if (A3D_ELEM(pMask, k, i, j) == 1)
166  ++NVoxelsOriginalMask;
167  if (i*i+j*j+k*k > (R-N_smoothing)*(R-N_smoothing))
168  A3D_ELEM(pMask, k, i, j) = -1;
169  }
170  Rparticle = round(sqrt(radius));
171  std::cout << "particle radius = " << Rparticle << std::endl;
172  size_t xrows = angles.mdimx;
173 
174  resolutionMatrix.initConstant(xrows, NVoxelsOriginalMask, maxRes);
175 
176 
177  #ifdef DEBUG_MASK
178  std::cout << "-------------DEBUG-----------" <<std::endl;
179  std::cout << "Next number ought to be the same than number of directions"
180  << std::endl;
181  std::cout << "number of angles" << xrows << std::endl;
182  std::cout << "---------END-DEBUG--------" <<std::endl;
183  mask.write("mask.vol");
184  #endif
185 
186  freq_fourier.initZeros(ZSIZE(inputVol));
187  int size = ZSIZE(inputVol);
188  maxRes = size;
189  minRes = 1;
190  V.clear();
191 
192 
193  double u;
194  int size_fourier(ZSIZE(fftV));
195 
196  VEC_ELEM(freq_fourier,0) = 1e-38;
197  for(size_t k=1; k<size_fourier; ++k)
198  {
199  FFT_IDX2DIGFREQ(k,size, u);
200  VEC_ELEM(freq_fourier,k) = u;
201  }
202 }
203 
204 
205 void ProgResDir::generateGridProjectionMatching(Matrix2D<double> &angles)
206 {
207  if (fastCompute == false)
208  {
209  angles.initZeros(2,81);
210  MAT_ELEM(angles, 0, 0) = 0.000000; MAT_ELEM(angles, 1, 0) = 0.000000;
211  MAT_ELEM(angles, 0, 1) = 36.000000; MAT_ELEM(angles, 1, 1) = 15.858741;
212  MAT_ELEM(angles, 0, 2) = 36.000000; MAT_ELEM(angles, 1, 2) = 31.717482;
213  MAT_ELEM(angles, 0, 3) = 36.000000; MAT_ELEM(angles, 1, 3) = 47.576224;
214  MAT_ELEM(angles, 0, 4) = 36.000000; MAT_ELEM(angles, 1, 4) = 63.434965;
215  MAT_ELEM(angles, 0, 5) = 62.494295; MAT_ELEM(angles, 1, 5) = -76.558393;
216  MAT_ELEM(angles, 0, 6) = 54.000000; MAT_ELEM(angles, 1, 6) = 90.000000;
217  MAT_ELEM(angles, 0, 7) = 45.505705; MAT_ELEM(angles, 1, 7) = 76.558393;
218  MAT_ELEM(angles, 0, 8) = 108.000000; MAT_ELEM(angles, 1, 8) = 15.858741;
219  MAT_ELEM(angles, 0, 9) = 108.000000; MAT_ELEM(angles, 1, 9) = 31.717482;
220  MAT_ELEM(angles, 0, 10) = 108.000000; MAT_ELEM(angles, 1, 10) = 47.576224;
221  MAT_ELEM(angles, 0, 11) = 108.000000; MAT_ELEM(angles, 1, 11) = 63.434965;
222  MAT_ELEM(angles, 0, 12) = 134.494295; MAT_ELEM(angles, 1, 12) = -76.558393;
223  MAT_ELEM(angles, 0, 13) = 126.000000; MAT_ELEM(angles, 1, 13) = 90.000000;
224  MAT_ELEM(angles, 0, 14) = 117.505705; MAT_ELEM(angles, 1, 14) = 76.558393;
225  MAT_ELEM(angles, 0, 15) = 144.000000; MAT_ELEM(angles, 1, 15) = -15.858741;
226  MAT_ELEM(angles, 0, 16) = 144.000000; MAT_ELEM(angles, 1, 16) = -31.717482;
227  MAT_ELEM(angles, 0, 17) = 144.000000; MAT_ELEM(angles, 1, 17) = -47.576224;
228  MAT_ELEM(angles, 0, 18) = 144.000000; MAT_ELEM(angles, 1, 18) = -63.434965;
229  MAT_ELEM(angles, 0, 19) = 170.494295; MAT_ELEM(angles, 1, 19) = 76.558393;
230  MAT_ELEM(angles, 0, 20) = 162.000000; MAT_ELEM(angles, 1, 20) = 90.000000;
231  MAT_ELEM(angles, 0, 21) = 153.505705; MAT_ELEM(angles, 1, 21) = -76.558393;
232  MAT_ELEM(angles, 0, 22) = 72.000000; MAT_ELEM(angles, 1, 22) = -15.858741;
233  MAT_ELEM(angles, 0, 23) = 72.000000; MAT_ELEM(angles, 1, 23) = -31.717482;
234  MAT_ELEM(angles, 0, 24) = 72.000000; MAT_ELEM(angles, 1, 24) = -47.576224;
235  MAT_ELEM(angles, 0, 25) = 72.000000; MAT_ELEM(angles, 1, 25) = -63.434965;
236  MAT_ELEM(angles, 0, 26) = 98.494295; MAT_ELEM(angles, 1, 26) = 76.558393;
237  MAT_ELEM(angles, 0, 27) = 90.000000; MAT_ELEM(angles, 1, 27) = 90.000000;
238  MAT_ELEM(angles, 0, 28) = 81.505705; MAT_ELEM(angles, 1, 28) = -76.558393;
239  MAT_ELEM(angles, 0, 29) = 0.000000; MAT_ELEM(angles, 1, 29) = -15.858741;
240  MAT_ELEM(angles, 0, 30) = 0.000000; MAT_ELEM(angles, 1, 30) = -31.717482;
241  MAT_ELEM(angles, 0, 31) = 0.000000; MAT_ELEM(angles, 1, 31) = -47.576224;
242  MAT_ELEM(angles, 0, 32) = 0.000000; MAT_ELEM(angles, 1, 32) = -63.434965;
243  MAT_ELEM(angles, 0, 33) = 26.494295; MAT_ELEM(angles, 1, 33) = 76.558393;
244  MAT_ELEM(angles, 0, 34) = 18.000000; MAT_ELEM(angles, 1, 34) = 90.000000;
245  MAT_ELEM(angles, 0, 35) = 9.505705; MAT_ELEM(angles, 1, 35) = -76.558393;
246  MAT_ELEM(angles, 0, 36) = 12.811021; MAT_ELEM(angles, 1, 36) = 42.234673;
247  MAT_ELEM(angles, 0, 37) = 18.466996; MAT_ELEM(angles, 1, 37) = 59.620797;
248  MAT_ELEM(angles, 0, 38) = 0.000000; MAT_ELEM(angles, 1, 38) = 90.000000;
249  MAT_ELEM(angles, 0, 39) = 8.867209; MAT_ELEM(angles, 1, 39) = 75.219088;
250  MAT_ELEM(angles, 0, 40) = 72.000000; MAT_ELEM(angles, 1, 40) = 26.565058;
251  MAT_ELEM(angles, 0, 41) = 59.188979; MAT_ELEM(angles, 1, 41) = 42.234673;
252  MAT_ELEM(angles, 0, 42) = 84.811021; MAT_ELEM(angles, 1, 42) = 42.234673;
253  MAT_ELEM(angles, 0, 43) = 53.533003; MAT_ELEM(angles, 1, 43) = 59.620797;
254  MAT_ELEM(angles, 0, 44) = 72.000000; MAT_ELEM(angles, 1, 44) = 58.282544;
255  MAT_ELEM(angles, 0, 45) = 90.466996; MAT_ELEM(angles, 1, 45) = 59.620797;
256  MAT_ELEM(angles, 0, 46) = 72.000000; MAT_ELEM(angles, 1, 46) = 90.000000;
257  MAT_ELEM(angles, 0, 47) = 63.132791; MAT_ELEM(angles, 1, 47) = 75.219088;
258  MAT_ELEM(angles, 0, 48) = 80.867209; MAT_ELEM(angles, 1, 48) = 75.219088;
259  MAT_ELEM(angles, 0, 49) = 144.000000; MAT_ELEM(angles, 1, 49) = 26.565058;
260  MAT_ELEM(angles, 0, 50) = 131.188979; MAT_ELEM(angles, 1, 50) = 42.234673;
261  MAT_ELEM(angles, 0, 51) = 156.811021; MAT_ELEM(angles, 1, 51) = 42.234673;
262  MAT_ELEM(angles, 0, 52) = 125.533003; MAT_ELEM(angles, 1, 52) = 59.620797;
263  MAT_ELEM(angles, 0, 53) = 144.000000; MAT_ELEM(angles, 1, 53) = 58.282544;
264  MAT_ELEM(angles, 0, 54) = 162.466996; MAT_ELEM(angles, 1, 54) = 59.620797;
265  MAT_ELEM(angles, 0, 55) = 144.000000; MAT_ELEM(angles, 1, 55) = 90.000000;
266  MAT_ELEM(angles, 0, 56) = 135.132791; MAT_ELEM(angles, 1, 56) = 75.219088;
267  MAT_ELEM(angles, 0, 57) = 152.867209; MAT_ELEM(angles, 1, 57) = 75.219088;
268  MAT_ELEM(angles, 0, 58) = 180.000000; MAT_ELEM(angles, 1, 58) = -26.565058;
269  MAT_ELEM(angles, 0, 59) = 167.188979; MAT_ELEM(angles, 1, 59) = -42.234673;
270  MAT_ELEM(angles, 0, 60) = 180.000000; MAT_ELEM(angles, 1, 60) = -58.282544;
271  MAT_ELEM(angles, 0, 61) = 161.533003; MAT_ELEM(angles, 1, 61) = -59.620797;
272  MAT_ELEM(angles, 0, 62) = 171.132791; MAT_ELEM(angles, 1, 62) = -75.219088;
273  MAT_ELEM(angles, 0, 63) = 108.000000; MAT_ELEM(angles, 1, 63) = -26.565058;
274  MAT_ELEM(angles, 0, 64) = 120.811021; MAT_ELEM(angles, 1, 64) = -42.234673;
275  MAT_ELEM(angles, 0, 65) = 95.188979; MAT_ELEM(angles, 1, 65) = -42.234673;
276  MAT_ELEM(angles, 0, 66) = 126.466996; MAT_ELEM(angles, 1, 66) = -59.620797;
277  MAT_ELEM(angles, 0, 67) = 108.000000; MAT_ELEM(angles, 1, 67) = -58.282544;
278  MAT_ELEM(angles, 0, 68) = 89.533003; MAT_ELEM(angles, 1, 68) = -59.620797;
279  MAT_ELEM(angles, 0, 69) = 108.000000; MAT_ELEM(angles, 1, 69) = 90.000000;
280  MAT_ELEM(angles, 0, 70) = 116.867209; MAT_ELEM(angles, 1, 70) = -75.219088;
281  MAT_ELEM(angles, 0, 71) = 99.132791; MAT_ELEM(angles, 1, 71) = -75.219088;
282  MAT_ELEM(angles, 0, 72) = 36.000000; MAT_ELEM(angles, 1, 72) = -26.565058;
283  MAT_ELEM(angles, 0, 73) = 48.811021; MAT_ELEM(angles, 1, 73) = -42.234673;
284  MAT_ELEM(angles, 0, 74) = 23.188979; MAT_ELEM(angles, 1, 74) = -42.234673;
285  MAT_ELEM(angles, 0, 75) = 54.466996; MAT_ELEM(angles, 1, 75) = -59.620797;
286  MAT_ELEM(angles, 0, 76) = 36.000000; MAT_ELEM(angles, 1, 76) = -58.282544;
287  MAT_ELEM(angles, 0, 77) = 17.533003; MAT_ELEM(angles, 1, 77) = -59.620797;
288  MAT_ELEM(angles, 0, 78) = 36.000000; MAT_ELEM(angles, 1, 78) = 90.000000;
289  MAT_ELEM(angles, 0, 79) = 44.867209; MAT_ELEM(angles, 1, 79) = -75.219088;
290  MAT_ELEM(angles, 0, 80) = 27.132791; MAT_ELEM(angles, 1, 80) = -75.219088;
291  }
292  else
293  {
294  angles.initZeros(2,47);
295  MAT_ELEM(angles, 0,1) =0; MAT_ELEM(angles, 1,1) =0;
296  MAT_ELEM(angles, 0,2) =36; MAT_ELEM(angles, 1,2) =21.145;
297  MAT_ELEM(angles, 0,3) =36; MAT_ELEM(angles, 1,3) =42.29;
298  MAT_ELEM(angles, 0,4) =36; MAT_ELEM(angles, 1,4) =63.435;
299  MAT_ELEM(angles, 0,5) =59.6043; MAT_ELEM(angles, 1,5) =-81.0207;
300  MAT_ELEM(angles, 0,6) =48.3957; MAT_ELEM(angles, 1,6) =81.0207;
301  MAT_ELEM(angles, 0,7) =108; MAT_ELEM(angles, 1,7) =21.145;
302  MAT_ELEM(angles, 0,8) =108; MAT_ELEM(angles, 1,8) =42.29;
303  MAT_ELEM(angles, 0,9) =108; MAT_ELEM(angles, 1,9) =63.435;
304  MAT_ELEM(angles, 0,10) =131.6043; MAT_ELEM(angles, 1,10) =-81.0207;
305  MAT_ELEM(angles, 0,11) =120.3957; MAT_ELEM(angles, 1,11) =81.0207;
306  MAT_ELEM(angles, 0,12) =144; MAT_ELEM(angles, 1,12) =-21.145;
307  MAT_ELEM(angles, 0,13) =144; MAT_ELEM(angles, 1,13) =-42.29;
308  MAT_ELEM(angles, 0,14) =144; MAT_ELEM(angles, 1,14) =-63.435;
309  MAT_ELEM(angles, 0,15) =167.6043; MAT_ELEM(angles, 1,15) =81.0207;
310  MAT_ELEM(angles, 0,16) =156.3957; MAT_ELEM(angles, 1,16) =-81.0207;
311  MAT_ELEM(angles, 0,17) =72; MAT_ELEM(angles, 1,17) =-21.145;
312  MAT_ELEM(angles, 0,18) =72; MAT_ELEM(angles, 1,18) =-42.29;
313  MAT_ELEM(angles, 0,19) =72; MAT_ELEM(angles, 1,19) =-63.435;
314  MAT_ELEM(angles, 0,20) =95.6043; MAT_ELEM(angles, 1,20) =81.0207;
315  MAT_ELEM(angles, 0,21) =84.3957; MAT_ELEM(angles, 1,21) =-81.0207;
316  MAT_ELEM(angles, 0,22) =0; MAT_ELEM(angles, 1,22) =-21.145;
317  MAT_ELEM(angles, 0,23) =0; MAT_ELEM(angles, 1,23) =-42.29;
318  MAT_ELEM(angles, 0,24) =0; MAT_ELEM(angles, 1,24) =-63.435;
319  MAT_ELEM(angles, 0,25) =23.6043; MAT_ELEM(angles, 1,25) =81.0207;
320  MAT_ELEM(angles, 0,26) =12.3957; MAT_ELEM(angles, 1,26) =-81.0207;
321  MAT_ELEM(angles, 0,27) =12.3756; MAT_ELEM(angles, 1,27) =58.8818;
322  MAT_ELEM(angles, 0,28) =72; MAT_ELEM(angles, 1,28) =36.349;
323  MAT_ELEM(angles, 0,29) =59.6244; MAT_ELEM(angles, 1,29) =58.8818;
324  MAT_ELEM(angles, 0,30) =84.3756; MAT_ELEM(angles, 1,30) =58.8818;
325  MAT_ELEM(angles, 0,31) =72; MAT_ELEM(angles, 1,31) =80.2161;
326  MAT_ELEM(angles, 0,32) =144; MAT_ELEM(angles, 1,32) =36.349;
327  MAT_ELEM(angles, 0,33) =131.6244; MAT_ELEM(angles, 1,33) =58.8818;
328  MAT_ELEM(angles, 0,34) =156.3756; MAT_ELEM(angles, 1,34) =58.8818;
329  MAT_ELEM(angles, 0,35) =144; MAT_ELEM(angles, 1,35) =80.2161;
330  MAT_ELEM(angles, 0,36) =180; MAT_ELEM(angles, 1,36) =-36.349;
331  MAT_ELEM(angles, 0,37) =167.6244; MAT_ELEM(angles, 1,37) =-58.8818;
332  MAT_ELEM(angles, 0,38) =180; MAT_ELEM(angles, 1,38) =-80.2161;
333  MAT_ELEM(angles, 0,39) =108; MAT_ELEM(angles, 1,39) =-36.349;
334  MAT_ELEM(angles, 0,40) =120.3756; MAT_ELEM(angles, 1,40) =-58.8818;
335  MAT_ELEM(angles, 0,41) =95.6244; MAT_ELEM(angles, 1,41) =-58.8818;
336  MAT_ELEM(angles, 0,42) =108; MAT_ELEM(angles, 1,42) =-80.2161;
337  MAT_ELEM(angles, 0,43) =36; MAT_ELEM(angles, 1,43) =-36.349;
338  MAT_ELEM(angles, 0,44) =48.3756; MAT_ELEM(angles, 1,44) =-58.8818;
339  MAT_ELEM(angles, 0,45) =23.6244; MAT_ELEM(angles, 1,45) =-58.8818;
340  MAT_ELEM(angles, 0,46) =36; MAT_ELEM(angles, 1,46) =-80.2161;
341  }
342 
343 }
344 
345 
346 
347 void ProgResDir::amplitudeMonogenicSignal3D_fast(const MultidimArray< std::complex<double> > &myfftV,
348  double freq, double freqH, double freqL, MultidimArray<double> &amplitude, int count, int dir, FileName fnDebug)
349 {
350  fftVRiesz.initZeros(myfftV);
351 // MultidimArray<double> coneVol;
352 // coneVol.initZeros(myfftV);
353  fftVRiesz_aux.initZeros(myfftV);
354  std::complex<double> J(0,1);
355 
356  // Filter the input volume and add it to amplitude
357  long n=0;
358  double ideltal=PI/(freq-freqH);
359 
360  for(size_t k=0; k<ZSIZE(myfftV); ++k)
361  {
362  for(size_t i=0; i<YSIZE(myfftV); ++i)
363  {
364  for(size_t j=0; j<XSIZE(myfftV); ++j)
365  {
366  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
367 // double iun = *ptriun;
368  double un=1.0/iun;
369  if (freqH<=un && un<=freq)
370  {
371  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
372 // DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= DIRECT_MULTIDIM_ELEM(conefilter, n);
373  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-freq)*ideltal));//H;
374  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = -J;
375  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= DIRECT_MULTIDIM_ELEM(fftVRiesz, n);
376  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= iun;
377 // DIRECT_MULTIDIM_ELEM(coneVol, n) = DIRECT_MULTIDIM_ELEM(conefilter, n);
378  } else if (un>freq)
379  {
380  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
381 // DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= DIRECT_MULTIDIM_ELEM(conefilter, n);
382  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = -J;
383  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= DIRECT_MULTIDIM_ELEM(fftVRiesz, n);
384  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= iun;
385 // DIRECT_MULTIDIM_ELEM(coneVol, n) = real(DIRECT_MULTIDIM_ELEM(myfftV, n)*conj(DIRECT_MULTIDIM_ELEM(myfftV, n)));
386  }
387  ++n;
388  }
389  }
390  }
391 
392 
393 // #ifdef DEBUG_DIR
396 // Image<double> direction;
397 // direction = coneVol;
398 // direction.write(formatString("cone_%i_%i.vol", dir+1, count));
400 // #endif
401 
402  transformer_inv.inverseFourierTransform(fftVRiesz, amplitude);
403 
404 // #ifdef DEBUG_DIR
405 // Image<double> filteredvolume;
406 // filteredvolume = VRiesz;
407 // filteredvolume.write(formatString("Volumen_filtrado_%i_%i.vol", dir+1,count));
408 // #endif
409 
410 
411 // amplitude.resizeNoCopy(VRiesz);
412  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
413  DIRECT_MULTIDIM_ELEM(amplitude,n) *= DIRECT_MULTIDIM_ELEM(amplitude,n);
414 
415 
416  // Calculate first component of Riesz vector
417  double ux;
418  n=0;
419  for(size_t k=0; k<ZSIZE(myfftV); ++k)
420  {
421  for(size_t i=0; i<YSIZE(myfftV); ++i)
422  {
423  for(size_t j=0; j<XSIZE(myfftV); ++j)
424  {
425  ux = VEC_ELEM(freq_fourier,j);
426  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = ux*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
427  ++n;
428  }
429  }
430  }
431 
432  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
433 
434  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
435  {
436  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
437  }
438 
439  // Calculate second and third component of Riesz vector
440  n=0;
441  double uy, uz;
442  for(size_t k=0; k<ZSIZE(myfftV); ++k)
443  {
444  uz = VEC_ELEM(freq_fourier,k);
445  for(size_t i=0; i<YSIZE(myfftV); ++i)
446  {
447  uy = VEC_ELEM(freq_fourier,i);
448  for(size_t j=0; j<XSIZE(myfftV); ++j)
449  {
450  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = uz*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
451  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = uy*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
452  ++n;
453  }
454  }
455  }
456 
457  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
458 
459  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
460  DIRECT_MULTIDIM_ELEM(amplitude,n) += DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
461 
462  transformer_inv.inverseFourierTransform(fftVRiesz_aux, VRiesz);
463 
464 
465 // amplitude.setXmippOrigin();
466  int z_size = ZSIZE(amplitude);
467  int siz = z_size*0.5;
468 
469  double limit_radius = (siz-N_smoothing);
470  n=0;
471  for(int k=0; k<z_size; ++k)
472  {
473  uz = (k - siz);
474  uz *= uz;
475  for(int i=0; i<z_size; ++i)
476  {
477  uy = (i - siz);
478  uy *= uy;
479  for(int j=0; j<z_size; ++j)
480  {
481  ux = (j - siz);
482  ux *= ux;
483  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
484  DIRECT_MULTIDIM_ELEM(amplitude,n)=sqrt(DIRECT_MULTIDIM_ELEM(amplitude,n));
485  double radius = sqrt(ux + uy + uz);
486  if ((radius>=limit_radius) && (radius<=siz))
487  DIRECT_MULTIDIM_ELEM(amplitude, n) *= 0.5*(1+cos(PI*(limit_radius-radius)/N_smoothing));
488  else if (radius>siz)
489  DIRECT_MULTIDIM_ELEM(amplitude, n) = 0;
490  ++n;
491  }
492  }
493  }
494 
495  //TODO: change (k - z_size*0.5)
496 
497 // #ifdef MONO_AMPLITUDE
498 // Image<double> saveImg2;
499 // saveImg2 = amplitude;
500 // if (fnDebug.c_str() != "")
501 // {
502 // FileName iternumber = formatString("smoothed_volume_%i_%i.vol", dir+1, count);
503 // saveImg2.write(fnDebug+iternumber);
504 // }
505 // saveImg2.clear();
506 // #endif
507 
508 
509  transformer_inv.FourierTransform(amplitude, fftVRiesz, false);
510 
511  double raised_w = PI/(freqL-freq);
512 
513  n=0;
514  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(fftVRiesz)
515  {
516  double un=1.0/DIRECT_MULTIDIM_ELEM(iu,n);
517  if (freqL>=un && un>=freq)
518  {
519  DIRECT_MULTIDIM_ELEM(fftVRiesz,n) *= 0.5*(1 + cos(raised_w*(un-freq)));
520  }
521  else
522  {
523  if (un>freqL)
524  {
525  DIRECT_MULTIDIM_ELEM(fftVRiesz,n) = 0;
526  }
527  }
528  }
529 
530  transformer_inv.inverseFourierTransform();
531 
532 // #ifdef MONO_AMPLITUDE
533 
534 // if (fnDebug.c_str() != "")
535 // {
536 // Image<double> saveImg2;
537 // saveImg2 = amplitude;
538 // FileName iternumber = formatString("_Filtered_Amplitude_%i_%i.vol", dir+1, count);
539 // saveImg2.write(fnDebug+iternumber);
540 // }
541 // #endif // DEBUG
542 }
543 
544 
545 void ProgResDir::defineCone(MultidimArray< std::complex<double> > &myfftV,
546  MultidimArray< std::complex<double> > &conefilter, double rot, double tilt)
547 {
548 // conefilter.initZeros(myfftV);
549  conefilter = myfftV;
550  // Filter the input volume and add it to amplitude
551 
552 // MultidimArray<double> conetest;
553 // conetest.initZeros(myfftV);
554 // #ifdef DEBUG_DIR
555 // MultidimArray<double> coneVol;
556 // coneVol.initZeros(iu);
557 // #endif
558 
559  double x_dir, y_dir, z_dir;
560 
561  x_dir = sin(tilt*PI/180)*cos(rot*PI/180);
562  y_dir = sin(tilt*PI/180)*sin(rot*PI/180);
563  z_dir = cos(tilt*PI/180);
564 
565 // double ang_con = 10*PI/180;
566  double ang_con = 15*PI/180;
567 
568  double uz, uy, ux;
569  long n = 0;
570  for(size_t k=0; k<ZSIZE(myfftV); ++k)
571  {
572  uz = VEC_ELEM(freq_fourier,k);
573  uz *= z_dir;
574  for(size_t i=0; i<YSIZE(myfftV); ++i)
575  {
576  uy = VEC_ELEM(freq_fourier,i);
577  uy *= y_dir;
578  for(size_t j=0; j<XSIZE(myfftV); ++j)
579  {
580  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
581  ux = VEC_ELEM(freq_fourier,j);
582  ux *= x_dir;
583 
584  //BE CAREFULL with the order
585  //double dotproduct = (uy*x_dir + ux*y_dir + uz*z_dir)*iun;
586  iun *= (ux + uy + uz);
587  double acosine = acos(fabs(iun));
588 // DIRECT_MULTIDIM_ELEM(conetest, n) = log(real(conj(DIRECT_MULTIDIM_ELEM(myfftV, n))*DIRECT_MULTIDIM_ELEM(myfftV, n)));
589  if (acosine>ang_con)
590  {
591  DIRECT_MULTIDIM_ELEM(conefilter, n) = 0;
592 // DIRECT_MULTIDIM_ELEM(conetest, n) = 0;
593  }
594 /*
595  //4822.53 mean a smoothed cone angle of 20 degrees
596  double arg_exp = acosine*acosine*acosine*acosine*4822.53;
597  DIRECT_MULTIDIM_ELEM(conefilter, n) *= exp(-arg_exp);
598 // DIRECT_MULTIDIM_ELEM(conetest, n) = exp(-arg_exp);
599  */
600  ++n;
601  }
602  }
603  }
604 //
605 // Image<double> saveImg2;
606 // saveImg2 = conetest;
607 // saveImg2.write("cono.vol");
608 
609 }
610 
611 void ProgResDir::diagSymMatrix3x3(Matrix2D<double> A,
612  Matrix1D<double> &eigenvalues, Matrix2D<double> &eigenvectors)
613 {
614  Matrix2D<double> B;
615  B.initZeros(3,3);
616 
617  MAT_ELEM(B,0,0) = 1;
618  MAT_ELEM(B,1,1) = 1;
619  MAT_ELEM(B,2,2) = 1;
620 
621  generalizedEigs(A, B, eigenvalues, eigenvectors);
622 }
623 
624 void ProgResDir::resolution2eval_(int &fourier_idx, double min_step,
625  double &resolution, double &last_resolution,
626  int &last_fourier_idx,
627  double &freq, double &freqL, double &freqH,
628  bool &continueIter, bool &breakIter, bool &doNextIteration)
629 {
630  int volsize = ZSIZE(VRiesz);
631 
632  FFT_IDX2DIGFREQ(fourier_idx, volsize, freq);
633 
634  resolution = sampling/freq;
635 // std::cout << "res = " << resolution << std::endl;
636 // std::cout << "min_step = " << min_step << std::endl;
637  if (resolution>8)
638  min_step =1;
639 
640 
641 
642  if ( fabs(resolution - last_resolution)<min_step )
643  {
644  freq = sampling/(last_resolution-min_step);
645  DIGFREQ2FFT_IDX(freq, volsize, fourier_idx);
646  FFT_IDX2DIGFREQ(fourier_idx, volsize, freq);
647 
648  if (fourier_idx == last_fourier_idx)
649  {
650  continueIter = true;
651  ++fourier_idx;
652  return;
653  }
654  }
655 
656  resolution = sampling/freq;
657  last_resolution = resolution;
658 
659  double step = 0.05*resolution;
660 
661  double resolution_L, resolution_H;
662 
663  if ( step < min_step)
664  {
665  resolution_L = resolution - min_step;
666  resolution_H = resolution + min_step;
667  }
668  else
669  {
670  resolution_L = 0.95*resolution;
671  resolution_H = 1.05*resolution;
672  }
673 
674  freqH = sampling/resolution_H;
675  freqL = sampling/resolution_L;
676 
677  if (freqH>0.5 || freqH<0)
678  freqH = 0.5;
679 
680  if (freqL>0.5 || freqL<0)
681  freqL = 0.5;
682  int fourier_idx_H, fourier_idx_L;
683 
684  DIGFREQ2FFT_IDX(freqH, volsize, fourier_idx_H);
685  DIGFREQ2FFT_IDX(freqL, volsize, fourier_idx_L);
686 
687  if (fourier_idx_H == fourier_idx)
688  fourier_idx_H = fourier_idx - 1;
689 
690  if (fourier_idx_L == fourier_idx)
691  fourier_idx_L = fourier_idx + 1;
692 
693  FFT_IDX2DIGFREQ(fourier_idx_H, volsize, freqH);
694  FFT_IDX2DIGFREQ(fourier_idx_L, volsize, freqL);
695 
696 // std::cout << "freq_H = " << freqH << std::endl;
697 // std::cout << "freq_L = " << freqL << std::endl;
698 
699  if (freq>0.49 || freq<0)
700  {
701  std::cout << "Nyquist limit reached" << std::endl;
702  breakIter = true;
703  doNextIteration = false;
704  return;
705  }
706  else
707  {
708  breakIter = false;
709  doNextIteration = true;
710  }
711 // std::cout << "resolution = " << resolution << " resolutionL = " <<
712 // sampling/(freqL) << " resolutionH = " << sampling/freqH
713 // << " las_res = " << last_resolution << std::endl;
714  last_fourier_idx = fourier_idx;
715  ++fourier_idx;
716 }
717 
718 
719 void ProgResDir::removeOutliers(Matrix2D<double> &resolutionMat)
720 {
721  std::cout << "Removing outliers..." << std::endl;
722 
723  double x1, y1, z1, x2, y2, z2, distance, resolution, sigma,
724  rot, tilt, threshold, sigma2, lastMinDistance;
725  double meandistance = 0, distance_2 = 0;
726  int numberdirections = angles.mdimx, N=0, count = 0;
727 
728  double criticalZ = icdf_gauss(significance);
729  double threshold_gauss;
730  double counter = 0;
731 
732  double ang = 20.0;
733 
734  Matrix2D<double> neigbour_dir;
735 
736  for (int k = 0; k<NVoxelsOriginalMask; ++k)
737  {
738  meandistance = 0;
739  sigma = 0;
740  distance_2 = 0;
741  counter = 0;
742 
743 // std::vector<double> neighbours;
744  neigbour_dir.initZeros(numberdirections, 2);
745  //Computing closest neighbours and its mean distance
746  for (int i = 0; i<numberdirections; ++i)
747  {
748 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
749 // {
750 // std::cout << k << " " << MAT_ELEM(resolutionMat, i, k) << " " << MAT_ELEM(trigProducts, 0, i) << " " <<
751 // MAT_ELEM(trigProducts, 1, i) << " " << MAT_ELEM(trigProducts, 2, i) << ";" << std::endl;
752 // }
753 
754  double resi = MAT_ELEM(resolutionMat, i, k);
755  if (resi>0)
756  {
757  x1 = MAT_ELEM(trigProducts, 0, i);
758  y1 = MAT_ELEM(trigProducts, 1, i);
759  z1 = MAT_ELEM(trigProducts, 2, i);
760 
761  for (int j = 0; j<numberdirections; ++j)
762  {
763  if (i != j)
764  {
765  x2 = MAT_ELEM(trigProducts, 0, j);
766  y2 = MAT_ELEM(trigProducts, 1, j);
767  z2 = MAT_ELEM(trigProducts, 2, j);
768  double resj = MAT_ELEM(resolutionMat, j, k);
769  if (resj>0)
770  {
771  distance = (180/PI)*acos(x1*x2 + y1*y2 + z1*z2);
772 
773  if (distance < ang)
774  {
775  x2 *= resj;
776  y2 *= resj;
777  z2 *= resj;
778  distance = sqrt( (resi*x1-x2)*(resi*x1-x2) + (resi*y1-y2)*(resi*y1-y2) + (resi*z1-z2)*(resi*z1-z2) );
779  MAT_ELEM(neigbour_dir, i, 0) += distance;
780  MAT_ELEM(neigbour_dir, i, 1) += 1;
781 // neighbours.push_back(distance);
782  meandistance += distance;
783  ++counter;
784  sigma += distance*distance;
785  }
786  }
787  }
788  }
789  }
790  }
791 
792  double thresholdDirection;
793 
794  meandistance= meandistance/counter;
795  sigma = sigma/counter - meandistance*meandistance;
796 
797 // std::sort(neighbours.begin(), neighbours.end());
798 // thresholdDirection = neighbours[size_t(neighbours.size()*significance)];
799  threshold_gauss = meandistance + criticalZ*sqrt(sigma);
800 
801 // neighbours.clear();
802 
803  //A direction is an outlier if is significative higher than overal distibution
804 
805  for (int i = 0; i<numberdirections; ++i)
806  {
807  double meandistance;
808  if (MAT_ELEM(neigbour_dir, i, 0)>0)
809  {
810  meandistance = MAT_ELEM(neigbour_dir, i, 0)/MAT_ELEM(neigbour_dir, i, 1);
811 
812  if ((meandistance>threshold_gauss) || MAT_ELEM(neigbour_dir, i, 0) <= 1)
813  {
814  MAT_ELEM(resolutionMat, i, k)=-1;
815  }
816  }
817 
818  }
819 
820 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
821 // std::cout << "threshold_gauss--------------=" << threshold_gauss << std::endl;
822 
823  }
824 }
825 
826 
827 
828 void ProgResDir::ellipsoidFitting(Matrix2D<double> &resolutionMat,
829  Matrix2D<double> &axis)
830 {
831 
832  std::cout << "FITTIG" << std::endl;
833  double x, y, z, a, b, c, resolution, rot, tilt;
834  int numberdirections = angles.mdimx;
835  std::vector<double> list_distances;
836 
837  //std::cout << "xrows = " << xrows << std::endl;
838 
839  //MAT_ELEM(resolutionMat, direccion, resolucion)
840 
841  Matrix2D<double> ellipMat;
842  int dimMatrix = 0;
843  size_t mycounter;
844  Matrix2D<double> pseudoinv, quadricMatrix;
845  Matrix1D<double> onesVector, leastSquares, residuals, residualssorted;
846  Matrix2D<double> eigenvectors;
847  Matrix1D<double> eigenvalues;
848  //rows 1 2 3 (length axis a b c- where a is the smallest one)
849  //rows 4 5 6 x y z coordinates of the first eigenvector
850  //rows 7 8 9 x y z coordinates of the second eigenvector
851  //rows 10 11 12 x y z coordinates of the third eigenvector
852  axis.initZeros(12, NVoxelsOriginalMask);
853 
854  quadricMatrix.initZeros(3,3);
855  for (int k = 0; k<NVoxelsOriginalMask; ++k)
856  {
857  dimMatrix = 0;
858  for (int i = 0; i<numberdirections; ++i)
859  {
860  if (MAT_ELEM(resolutionMat, i, k) > 0)
861  ++dimMatrix;
862  }
863 
864  ellipMat.initZeros(dimMatrix, 6);
865  mycounter = 0; //It is required to store the matrix ellipMat
866  for (int i = 0; i<numberdirections; ++i)
867  {
868  resolution = MAT_ELEM(resolutionMat, i, k);
869 
870  if (resolution>0)
871  {
872  x = resolution*MAT_ELEM(trigProducts, 0, i);
873  y = resolution*MAT_ELEM(trigProducts, 1, i);
874  z = resolution*MAT_ELEM(trigProducts, 2, i);
875 
876  MAT_ELEM(ellipMat, mycounter, 0) = x*x;
877  MAT_ELEM(ellipMat, mycounter, 1) = y*y;
878  MAT_ELEM(ellipMat, mycounter, 2) = z*z;
879  MAT_ELEM(ellipMat, mycounter, 3) = 2*x*y;
880  MAT_ELEM(ellipMat, mycounter, 4) = 2*x*z;
881  MAT_ELEM(ellipMat, mycounter, 5) = 2*y*z;
882  ++mycounter;
883  }
884  }
885 
886 
887  ellipMat.inv(pseudoinv);
888 
889  onesVector.initConstant(mycounter, 1.0);
890  leastSquares = pseudoinv*onesVector;
891 
892  //Removing outliers
893  residuals = ellipMat*leastSquares - onesVector;
894  residualssorted = residuals.sort();
895 
896  double threshold_plus = VEC_ELEM(residualssorted, size_t(residualssorted.size()*significance));
897  double threshold_minus = VEC_ELEM(residualssorted, size_t(residualssorted.size()*(1.0-significance)));
898 
899  mycounter = 0;
900  size_t ellipsoidcounter = 0;
901  for (int i = 0; i<numberdirections; ++i)
902  {
903  resolution = MAT_ELEM(resolutionMat, i, k);
904 
905  if (resolution>0)
906  {
907  if ( (VEC_ELEM(residuals, mycounter) > threshold_plus) ||
908  (VEC_ELEM(residuals, mycounter) < threshold_minus) )
909  {
910  MAT_ELEM(resolutionMat, i, k) = -1;
911  }
912  else
913  {
914  ++ellipsoidcounter;
915  }
916  ++mycounter;
917  }
918  }
919 
920  ellipMat.initZeros(ellipsoidcounter, 6);
921  mycounter = 0; //It is required to store the matrix ellipMat
922  for (int i = 0; i<numberdirections; ++i)
923  {
924  resolution = MAT_ELEM(resolutionMat, i, k);
925 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
926 // {
927 // std::cout << k << " " << resolution << " " << MAT_ELEM(trigProducts, 0, i)
928 // << " " << MAT_ELEM(trigProducts, 1, i)
929 // << " " << MAT_ELEM(trigProducts, 2, i) << ";" << std::endl;
930 // }
931 
932  if (resolution>0)
933  {
934  x = resolution*MAT_ELEM(trigProducts, 0, i);
935  y = resolution*MAT_ELEM(trigProducts, 1, i);
936  z = resolution*MAT_ELEM(trigProducts, 2, i);
937 
938  MAT_ELEM(ellipMat, mycounter, 0) = x*x;
939  MAT_ELEM(ellipMat, mycounter, 1) = y*y;
940  MAT_ELEM(ellipMat, mycounter, 2) = z*z;
941  MAT_ELEM(ellipMat, mycounter, 3) = 2*x*y;
942  MAT_ELEM(ellipMat, mycounter, 4) = 2*x*z;
943  MAT_ELEM(ellipMat, mycounter, 5) = 2*y*z;
944  ++mycounter;
945  }
946  }
947 
948  // defining ellipsoid
949 
950  ellipMat.inv(pseudoinv);
951 
952  onesVector.initConstant(mycounter, 1.0);
953  leastSquares = pseudoinv*onesVector;
954 
955 
956  MAT_ELEM(quadricMatrix, 0, 0) = VEC_ELEM(leastSquares, 0);
957  MAT_ELEM(quadricMatrix, 0, 1) = VEC_ELEM(leastSquares, 3);
958  MAT_ELEM(quadricMatrix, 0, 2) = VEC_ELEM(leastSquares, 4);
959  MAT_ELEM(quadricMatrix, 1, 0) = VEC_ELEM(leastSquares, 3);
960  MAT_ELEM(quadricMatrix, 1, 1) = VEC_ELEM(leastSquares, 1);
961  MAT_ELEM(quadricMatrix, 1, 2) = VEC_ELEM(leastSquares, 5);
962  MAT_ELEM(quadricMatrix, 2, 0) = VEC_ELEM(leastSquares, 4);
963  MAT_ELEM(quadricMatrix, 2, 1) = VEC_ELEM(leastSquares, 5);
964  MAT_ELEM(quadricMatrix, 2, 2) = VEC_ELEM(leastSquares, 2);
965 
966  diagSymMatrix3x3(quadricMatrix, eigenvalues, eigenvectors);
967 
968  if (VEC_ELEM(eigenvalues, 0)<0)
969  {//This is de the vectorial product
970  VEC_ELEM(eigenvalues, 0) = VEC_ELEM(eigenvalues, 1);
971  MAT_ELEM(axis,3, k) = MAT_ELEM(eigenvectors,0,1)*MAT_ELEM(eigenvectors,2,2)-
972  MAT_ELEM(eigenvectors,2,1)*MAT_ELEM(eigenvectors,1,2);
973  MAT_ELEM(axis,4, k) = MAT_ELEM(eigenvectors,2,1)*MAT_ELEM(eigenvectors,0,2)-
974  MAT_ELEM(eigenvectors,0,1)*MAT_ELEM(eigenvectors,2,2);
975  MAT_ELEM(axis,5, k) = MAT_ELEM(eigenvectors,0,1)*MAT_ELEM(eigenvectors,1,2)-
976  MAT_ELEM(eigenvectors,1,1)*MAT_ELEM(eigenvectors,0,2);
977  }
978 
979  a = 1/sqrt(VEC_ELEM(eigenvalues, 0));
980  b = 1/sqrt(VEC_ELEM(eigenvalues, 1));
981  c = 1/sqrt(VEC_ELEM(eigenvalues, 2));
982 
983 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
984 // {
985 // std::cout << "a = " << a << std::endl;
986 // std::cout << "b = " << b << std::endl;
987 // std::cout << "c = " << c << std::endl;
988 // std::cout << "=--------------=" << std::endl;
989 // }
990 
991  MAT_ELEM(axis,0, k) = a;
992  MAT_ELEM(axis,1, k) = b;
993  MAT_ELEM(axis,2, k) = c;
994 
995  MAT_ELEM(axis,3, k) = MAT_ELEM(eigenvectors,0,0);
996  MAT_ELEM(axis,4, k) = MAT_ELEM(eigenvectors,1,0);
997  MAT_ELEM(axis,5, k) = MAT_ELEM(eigenvectors,2,0);
998 
999  MAT_ELEM(axis,6, k) = MAT_ELEM(eigenvectors,0,1);
1000  MAT_ELEM(axis,7, k) = MAT_ELEM(eigenvectors,1,1);
1001  MAT_ELEM(axis,8, k) = MAT_ELEM(eigenvectors,2,1);
1002 
1003  MAT_ELEM(axis,9, k) = MAT_ELEM(eigenvectors,0,2);
1004  MAT_ELEM(axis,10, k) = MAT_ELEM(eigenvectors,1,2);
1005  MAT_ELEM(axis,11, k) = MAT_ELEM(eigenvectors,2,2);
1006  }
1007 }
1008 
1009 
1010 void ProgResDir::radialAverageInMask(MultidimArray<int> &mask,
1011  MultidimArray<double> &inputVol_1, MultidimArray<double> &inputVol_2,
1012  MultidimArray<double> &inputVol_3, MultidimArray<double> &inputVol_4,
1013  MultidimArray<double> &inputVol_5, MetaDataVec &md)
1014 {
1015  double u_inf, u_sup, u;
1016 
1017  int step = 1;
1018 
1019  double N, NN;
1020  MultidimArray<double> radialAvg(XSIZE(inputVol_1)*0.5);
1021  int uk, uj, ui;
1022 
1023  int siz = XSIZE(inputVol_1);
1024  size_t objId;
1025 
1026  inputVol_1.setXmippOrigin();
1027  inputVol_2.setXmippOrigin();
1028  inputVol_3.setXmippOrigin();
1029  inputVol_4.setXmippOrigin();
1030  inputVol_5.setXmippOrigin();
1031  MultidimArray<double> zVolume;
1032  zVolume.initZeros(inputVol_1);
1033  zVolume.setXmippOrigin();
1034 
1035  MultidimArray<int> &pMask = mask;
1036  pMask.setXmippOrigin();
1037 
1038  Matrix2D<double> std_mean_Radial_1, std_mean_Radial_2, std_mean_Radial_3,
1039  std_mean_Radial_4, std_mean_Radial_5;
1040 
1041  std_mean_Radial_1.initZeros(2, (size_t) siz*0.5 + 1);
1042  std_mean_Radial_2.initZeros(2, (size_t) siz*0.5 + 1);
1043  std_mean_Radial_3.initZeros(2, (size_t) siz*0.5 + 1);
1044  std_mean_Radial_4.initZeros(2, (size_t) siz*0.5 + 1);
1045  std_mean_Radial_5.initZeros(2, (size_t) siz*0.5 + 1);
1046 
1047  for(size_t kk=1; kk<siz*0.5; ++kk)
1048  {
1049  double cum_mean_1 = 0;
1050  double cum_mean_2 = 0;
1051  double cum_mean_3 = 0;
1052  double cum_mean_4 = 0;
1053  double cum_mean_5 = 0;
1054 
1055  double cum2_mean_1 = 0;
1056  double cum2_mean_2 = 0;
1057  double cum2_mean_3 = 0;
1058  double cum2_mean_4 = 0;
1059  double cum2_mean_5 = 0;
1060 
1061  N = 0;
1062  u_sup = kk + step;
1063  u_inf = kk - step;
1064 
1065  FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
1066  {
1067  if (A3D_ELEM(pMask, k, i, j)>0)
1068  {
1069  u = sqrt(k*k + i*i + j*j);
1070  if ((u<u_sup) && (u>=u_inf))
1071  {
1072  double aux;
1073  aux = A3D_ELEM(inputVol_1, k, i, j);
1074  cum_mean_1 += aux;
1075  cum2_mean_1 += aux*aux;
1076  aux = A3D_ELEM(inputVol_2, k, i, j);
1077  cum_mean_2 += aux;
1078  cum2_mean_2 += aux*aux;
1079  aux = A3D_ELEM(inputVol_3, k, i, j);
1080  cum_mean_3 += aux;
1081  cum2_mean_3 += aux*aux;
1082  aux = A3D_ELEM(inputVol_4, k, i, j);
1083  cum_mean_4 += aux;
1084  cum2_mean_4 += aux*aux;
1085  aux = A3D_ELEM(inputVol_5, k, i, j);
1086  cum_mean_5 += aux;
1087  cum2_mean_5 += aux*aux;
1088 
1089  N = N + 1;
1090  }
1091  }
1092  }
1093 
1094  double sigma1, sigma2, sigma3, sigma4, sigma5;
1095 
1096 
1097 
1098 // if ((cum_mean_1==0) || (cum_mean_2==0) || (cum_mean_3==0) || (cum_mean_4==0) || (cum_mean_5==0))
1099 // {
1100 // md.setValue(MDL_IDX, kk, objId);
1101 // md.setValue(MDL_VOLUME_SCORE1, cum_mean_1, objId);
1102 // md.setValue(MDL_VOLUME_SCORE2, cum_mean_2, objId);
1103 // md.setValue(MDL_VOLUME_SCORE3, cum_mean_3, objId);
1104 // md.setValue(MDL_VOLUME_SCORE4, cum_mean_4, objId);
1105 // md.setValue(MDL_AVG, cum_mean_5, objId);
1106 // }
1107 // else
1108 // {
1109  if ((cum_mean_1>0) || (cum_mean_2>0) || (cum_mean_3>0) || (cum_mean_4>0) || (cum_mean_5>0))
1110  {
1111  objId = md.addObject();
1112  cum_mean_1 = (cum_mean_1/N);
1113  cum_mean_2 = (cum_mean_2/N);
1114  cum_mean_3 = (cum_mean_3/N);
1115  cum_mean_4 = (cum_mean_4/N);
1116  cum_mean_5 = (cum_mean_5/N);
1117 
1118  MAT_ELEM(std_mean_Radial_1,1, kk) = cum_mean_1;
1119  MAT_ELEM(std_mean_Radial_2,1, kk) = cum_mean_2;
1120  MAT_ELEM(std_mean_Radial_3,1, kk) = cum_mean_3;
1121  MAT_ELEM(std_mean_Radial_4,1, kk) = cum_mean_4;
1122  MAT_ELEM(std_mean_Radial_5,1, kk) = cum_mean_5;
1123 
1124  md.setValue(MDL_IDX, kk, objId);
1125  md.setValue(MDL_VOLUME_SCORE1, cum_mean_1, objId);
1126  md.setValue(MDL_VOLUME_SCORE2, cum_mean_2, objId);
1127  md.setValue(MDL_VOLUME_SCORE3, cum_mean_3, objId);
1128  md.setValue(MDL_VOLUME_SCORE4, cum_mean_4, objId);
1129  md.setValue(MDL_AVG, cum_mean_5, objId);
1130 
1131  MAT_ELEM(std_mean_Radial_1,0, kk) = sqrt(cum2_mean_1/N - cum_mean_1*cum_mean_1);
1132  MAT_ELEM(std_mean_Radial_2,0, kk) = sqrt(cum2_mean_2/N - cum_mean_2*cum_mean_2);
1133  MAT_ELEM(std_mean_Radial_3,0, kk) = sqrt(cum2_mean_3/N - cum_mean_3*cum_mean_3);
1134  MAT_ELEM(std_mean_Radial_4,0, kk) = sqrt(cum2_mean_4/N - cum_mean_4*cum_mean_4+0.001);
1135  MAT_ELEM(std_mean_Radial_5,0, kk) = sqrt(cum2_mean_5/N - cum_mean_5*cum_mean_5);
1136 
1137 // std::cout << "cum2_mean_4/(N) = " << cum2_mean_4/(N) << " " << "cum_mean_4*cum_mean_4) = " << cum_mean_4*cum_mean_4 << std::endl;
1138 // std::cout << "MAT_ELEM(std_mean_Radial_1,0, kk) = " << MAT_ELEM(std_mean_Radial_1,0, kk) << " " << MAT_ELEM(std_mean_Radial_2,0, kk) << " " << MAT_ELEM(std_mean_Radial_3,0, kk) << " " << MAT_ELEM(std_mean_Radial_4,0, kk) << " " << MAT_ELEM(std_mean_Radial_5,0, kk)<< std::endl;
1139 // std::cout << "MAT_ELEM(std_mean_Radial_1,1, kk) = " << MAT_ELEM(std_mean_Radial_1,1, kk) << " " << MAT_ELEM(std_mean_Radial_2,1, kk) << " " << MAT_ELEM(std_mean_Radial_3,1, kk) << " " << MAT_ELEM(std_mean_Radial_4,1, kk) << " " << MAT_ELEM(std_mean_Radial_5,1, kk)<< std::endl;
1140  }
1141 
1142 
1143  double lastz;
1144  FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
1145  {
1146  if (A3D_ELEM(pMask, k, i, j)>0)
1147  {
1148  lastz = -1;
1149  u = sqrt(k*k + i*i + j*j);
1150  if ((u<u_sup) && (u>=u_inf) && (MAT_ELEM(std_mean_Radial_1,1, kk)>0))
1151  {
1152  double z, aux;
1153  aux = abs((A3D_ELEM(inputVol_1, k, i, j) - MAT_ELEM(std_mean_Radial_1,1, kk))/MAT_ELEM(std_mean_Radial_1,0, kk)) + + 0.002;
1154  if (aux > lastz)
1155  lastz = aux;
1156  aux = abs((A3D_ELEM(inputVol_2, k, i, j) - MAT_ELEM(std_mean_Radial_2,1, kk))/MAT_ELEM(std_mean_Radial_2,0, kk)) + 0.002;
1157  if (aux > lastz)
1158  lastz = aux;
1159  aux = abs((A3D_ELEM(inputVol_3, k, i, j) - MAT_ELEM(std_mean_Radial_3,1, kk))/MAT_ELEM(std_mean_Radial_3,0, kk)) + 0.002;
1160  if (aux > lastz)
1161  lastz = aux;
1162 // aux = abs((A3D_ELEM(inputVol_4, k, i, j) - MAT_ELEM(std_mean_Radial_4,1, kk))/MAT_ELEM(std_mean_Radial_4,0, kk));
1163 // if (aux > lastz)
1164 // lastz = aux;
1165  aux = abs((A3D_ELEM(inputVol_5, k, i, j) - MAT_ELEM(std_mean_Radial_5,1, kk))/MAT_ELEM(std_mean_Radial_5,0, kk)) + 0.002;
1166  if (aux > lastz)
1167  lastz = aux;
1168 
1169  if (lastz>5.0) //This line considers zscores higher than 5 as 5(saturated at 5sigma)
1170  lastz = 5;
1171 
1172  A3D_ELEM(zVolume, k, i, j) = lastz;
1173  }
1174  }
1175  }
1176  }
1177  Image<double> zVolumesave;
1178  zVolumesave = zVolume;
1179  zVolumesave.write(fnZscore);
1180 }
1181 
1182 void ProgResDir::radialAzimuthalResolution(Matrix2D<double> &resolutionMat,
1183  MultidimArray<int> &pmask,
1184  MultidimArray<double> &radial,
1185  MultidimArray<double> &azimuthal,
1186 // MultidimArray<double> &meanResolution,
1187  MultidimArray<double> &lowestResolution,
1188  MultidimArray<double> &highestResolution,
1189  MultidimArray<double> &doaResolution_1,
1190  MultidimArray<double> &doaResolution_2,
1191  double &radial_Thr, double &azimuthal_Thr,
1192  MetaDataVec &mdprefDirs)
1193 {
1194 
1195  radial.initZeros(pmask);
1196  azimuthal.initZeros(pmask);
1197  lowestResolution.initZeros(pmask);
1198  highestResolution.initZeros(pmask);
1199  doaResolution_1.initZeros(pmask);
1200  doaResolution_2.initZeros(pmask);
1201 
1202  double radial_angle = 45*PI/180;
1203  double azimuthal_resolution = 0;
1204  double radial_resolution = 0;
1205  double azimuthal_angle = 70*PI/180;
1206  double resolution, dotproduct, x, y, z, iu, arcos;
1207  int xrows = angles.mdimx;
1208  int idx = 0;
1209 
1210  double count_radial = 0.0, count_azimuthal = 0.0;
1211 
1212  Matrix1D<int> PrefferredDirHist, resolutionMeanVector;
1213  PrefferredDirHist.initZeros(xrows);
1214  resolutionMeanVector.initZeros(xrows);
1215 
1216 // meanResolution.initZeros(pmask);
1217  size_t objId;
1218  FOR_ALL_ELEMENTS_IN_ARRAY3D(pmask)
1219  {
1220  //i defines the direction and k the voxel
1221  if (A3D_ELEM(pmask,k,i,j) > 0 )
1222  {
1223  iu = 1/sqrt(i*i + j*j + k*k);
1224  count_radial = 0;
1225  count_azimuthal = 0;
1226  std::vector<double> ResList;
1227 
1228  double lastRes = 100; //A non-sense value
1229 
1230  for (int ii = 0; ii<xrows; ++ii)
1231  {
1232  resolution = MAT_ELEM(resolutionMat, ii, idx);
1233 
1234  if (resolution>0)
1235  {
1236  ResList.push_back(resolution);
1237  x = MAT_ELEM(trigProducts, 0, ii);
1238  y = MAT_ELEM(trigProducts, 1, ii);
1239  z = MAT_ELEM(trigProducts, 2, ii);
1240 
1241  dotproduct = (x*i + y*j + z*k)*iu;
1242  arcos = acos(fabs(dotproduct));
1243  if (arcos>=azimuthal_angle)
1244  {
1245  count_azimuthal = count_azimuthal + 1;
1246  azimuthal_resolution += resolution;
1247  }
1248  if (arcos<=radial_angle)
1249  {
1250  count_radial = count_radial + 1;
1251  radial_resolution += resolution;
1252  }
1253 
1254  }
1255 
1256  }
1257 
1258 // std::cout << "count_radial = " << count_radial << std::endl;
1259 // std::cout << "count_azimuthal = " << count_azimuthal << std::endl;
1260 // std::cout << " " << std::endl;
1261 
1262 // A3D_ELEM(meanResolution,k,i,j) = meanRes[(size_t) floor(0.5*meanRes.size())];
1263  std::sort(ResList.begin(),ResList.end());
1264 
1265  double Mres, mres, medianResol, res75, res25;
1266 
1267  Mres = ResList[ (size_t) floor(0.95*ResList.size()) ];
1268  A3D_ELEM(lowestResolution,k,i,j) = Mres;
1269 
1270  mres = ResList[ (size_t) floor(0.05*ResList.size()) ];
1271 
1272  res75 = ResList[ (size_t) floor(0.83*ResList.size()) ];
1273  res25 = ResList[ (size_t) floor(0.17*ResList.size()) ];
1274 
1275  A3D_ELEM(highestResolution,k,i,j) = mres;
1276 
1277  A3D_ELEM(doaResolution_1,k,i,j) = 0.5*(res75 - res25);//( (Mres - mres)/(Mres + mres) );
1278 
1279  A3D_ELEM(doaResolution_2,k,i,j) = 0.5*( Mres + mres );
1280 
1281  ResList.clear();
1282 
1283  //Prefferred directions
1284 
1285  for (int ii = 0; ii<xrows; ++ii)
1286  {
1287  resolution = MAT_ELEM(resolutionMat, ii, idx);
1288 
1289  if (resolution>0)
1290  {
1291  if ((mres>(resolution-0.1)) && (mres<(resolution+0.1)))
1292  {
1293  VEC_ELEM(PrefferredDirHist,ii) += 1;
1294  VEC_ELEM(resolutionMeanVector,ii) += resolution;
1295  }
1296  }
1297 
1298  }
1299  ++idx;
1300  }
1301 
1302  if (count_radial<1)
1303  A3D_ELEM(radial,k,i,j) = A3D_ELEM(doaResolution_2,k,i,j);
1304  else
1305  A3D_ELEM(radial,k,i,j) = radial_resolution/count_radial;
1306 
1307  if (count_azimuthal<1)
1308  A3D_ELEM(azimuthal,k,i,j) = A3D_ELEM(doaResolution_2,k,i,j);
1309  else
1310  A3D_ELEM(azimuthal,k,i,j) = azimuthal_resolution/count_azimuthal;
1311 
1312 
1313  azimuthal_resolution = 0;
1314  radial_resolution = 0;
1315  }
1316 
1317 
1318  for (size_t ii = 0; ii<xrows; ++ii)
1319  {
1320  objId = mdprefDirs.addObject();
1321  size_t con;
1322  con = (size_t) VEC_ELEM(PrefferredDirHist,ii);
1323 // std::cout << ii << " " << con << std::endl;
1324  double rot = MAT_ELEM(angles, 0, ii);
1325  double tilt = MAT_ELEM(angles, 1, ii);
1326 
1327  if (tilt<0)
1328  {
1329  tilt = abs(tilt);
1330  rot = rot + 180;
1331  }
1332  double meanRes;
1333 
1334  meanRes = VEC_ELEM(resolutionMeanVector,ii)/((double) VEC_ELEM(PrefferredDirHist,ii));
1335 
1336  mdprefDirs.setValue(MDL_ANGLE_ROT, rot, objId);
1337  mdprefDirs.setValue(MDL_ANGLE_TILT, tilt, objId);
1338  mdprefDirs.setValue(MDL_WEIGHT, (double) con, objId);
1339  mdprefDirs.setValue(MDL_RESOLUTION_FREQ, meanRes, objId);
1340  mdprefDirs.setValue(MDL_X, (double) ii, objId);
1341  mdprefDirs.setValue(MDL_COUNT, con, objId);
1342  }
1343  mdprefDirs.write(fnprefMin);
1344 
1345  std::vector<double> radialList, azimuthalList;
1346 
1347  FOR_ALL_ELEMENTS_IN_ARRAY3D(radial)
1348  {
1349  //i defines the direction and k the voxel
1350  if (A3D_ELEM(pmask,k,i,j) > 0 )
1351  {
1352  radialList.push_back(A3D_ELEM(radial,k,i,j));
1353  azimuthalList.push_back(A3D_ELEM(azimuthal,k,i,j));
1354  }
1355  }
1356 
1357  std::sort(radialList.begin(),radialList.end());
1358  std::sort(azimuthalList.begin(),azimuthalList.end());
1359 
1360  radial_Thr = radialList[(size_t) floor(radialList.size()*0.9)];
1361  azimuthal_Thr = azimuthalList[(size_t) floor(azimuthalList.size()*0.9)];
1362 }
1363 
1364 //TODO: change this function to be more efficient
1365 double ProgResDir::firstMonoResEstimation(MultidimArray< std::complex<double> > &myfftV,
1366  double freq, double freqH, MultidimArray<double> &amplitude)
1367 {
1368  fftVRiesz.initZeros(myfftV);
1369  amplitude.resizeNoCopy(VRiesz);
1370  std::complex<double> J(0,1);
1371 
1372  // Filter the input volume and add it to amplitude
1373  long n=0;
1374  double iw=1.0/freq;
1375  double iwl=1.0/freqH;
1376  double ideltal=PI/(freq-freqH);
1377 
1378  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(myfftV)
1379  {
1380  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1381  double un=1.0/iun;
1382  if (freqH<=un && un<=freq)
1383  {
1384  //double H=0.5*(1+cos((un-w1)*ideltal));
1385  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1386  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-freq)*ideltal));//H;
1387  } else if (un>freq)
1388  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1389  }
1390 
1391  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1392 
1393  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1394  DIRECT_MULTIDIM_ELEM(amplitude,n)=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1395 
1396  // Calculate first component of Riesz vector
1397  fftVRiesz.initZeros(myfftV);
1398  double uz, uy, ux;
1399  n=0;
1400 
1401  for(size_t k=0; k<ZSIZE(myfftV); ++k)
1402  {
1403  for(size_t i=0; i<YSIZE(myfftV); ++i)
1404  {
1405  for(size_t j=0; j<XSIZE(myfftV); ++j)
1406  {
1407  ux = VEC_ELEM(freq_fourier,j);
1408  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1409  double un=1.0/iun;
1410  if (freqH<=un && un<=freq)
1411  {
1412  //double H=0.5*(1+cos((un-w1)*ideltal));
1413  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-J*ux*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1414  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *=0.5*(1+cos((un-w1)*ideltal));//H;
1415  //Next lines are an optimization of the commented ones
1416  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1417  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1418  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= -ux*iun*0.5*(1+cos((un-freq)*ideltal));//H;
1419  } else if (un>freq)
1420  {
1421  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-ux*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1422  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1423  }
1424  ++n;
1425  }
1426  }
1427  }
1428  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1429  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1430  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1431 
1432  // Calculate second component of Riesz vector
1433  fftVRiesz.initZeros(myfftV);
1434  n=0;
1435 
1436  for(size_t k=0; k<ZSIZE(myfftV); ++k)
1437  {
1438  for(size_t i=0; i<YSIZE(myfftV); ++i)
1439  {
1440  uy = VEC_ELEM(freq_fourier,i);
1441  for(size_t j=0; j<XSIZE(myfftV); ++j)
1442  {
1443  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1444  double un=1.0/iun;
1445  if (freqH<=un && un<=freq)
1446  {
1447  //double H=0.5*(1+cos((un-w1)*ideltal));
1448  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-J*uy*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1449  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *=0.5*(1+cos((un-w1)*ideltal));//H;
1450  //Next lines are an optimization of the commented ones
1451  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1452  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1453  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= -uy*iun*0.5*(1+cos((un-freq)*ideltal));//H;
1454  } else if (un>freq)
1455  {
1456  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-uy*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1457  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1458  }
1459  ++n;
1460  }
1461  }
1462  }
1463  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1464  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1465  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1466 
1467  // Calculate third component of Riesz vector
1468  fftVRiesz.initZeros(myfftV);
1469  n=0;
1470  for(size_t k=0; k<ZSIZE(myfftV); ++k)
1471  {
1472  uz = VEC_ELEM(freq_fourier,k);
1473  for(size_t i=0; i<YSIZE(myfftV); ++i)
1474  {
1475  for(size_t j=0; j<XSIZE(myfftV); ++j)
1476  {
1477  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1478  double un=1.0/iun;
1479  if (freqH<=un && un<=freq)
1480  {
1481  //double H=0.5*(1+cos((un-w1)*ideltal));
1482  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-J*uz*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1483  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-w1)*ideltal));//H;
1484  //Next lines are an optimization of the commented ones
1485  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1486  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1487  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= -uz*iun*0.5*(1+cos((un-freq)*ideltal));//H;
1488  } else if (un>freq)
1489  {
1490  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-uz*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1491  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1492  }
1493  ++n;
1494  }
1495  }
1496  }
1497  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1498  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1499  {
1500  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1501  DIRECT_MULTIDIM_ELEM(amplitude,n)=sqrt(DIRECT_MULTIDIM_ELEM(amplitude,n));
1502  }
1503 
1504  // Low pass filter the monogenic amplitude
1505  // Prepare low pass filter
1506  FourierFilter lowPassFilter, FilterBand;
1507  lowPassFilter.FilterShape = RAISED_COSINE;
1508  lowPassFilter.raised_w = 0.01;
1509  lowPassFilter.do_generate_3dmask = false;
1510  lowPassFilter.FilterBand = LOWPASS;
1511  lowPassFilter.w1 = freq;
1512  amplitude.setXmippOrigin();
1513  lowPassFilter.applyMaskSpace(amplitude);
1514 
1515  double sumS=0, sumS2=0, sumN=0, sumN2=0, NN = 0, NS = 0;
1516  MultidimArray<int> &pMask = mask();
1517  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1518  {
1519  double amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitude, n);
1520  if (DIRECT_MULTIDIM_ELEM(pMask, n)==0)
1521  {
1522  sumN += amplitudeValue;
1523  sumN2 += amplitudeValue*amplitudeValue;
1524  ++NN;
1525  }
1526  }
1527 
1528  double mean_noise = sumN/NN;
1529  return mean_noise;
1530 
1531 }
1532 
1533 void ProgResDir::run()
1534 {
1535  produceSideInfo();
1536  bool continueIter = false, breakIter = false;
1537  double criticalZ=icdf_gauss(significance);
1538 
1539  double step;
1540  step = res_step;
1541 
1542  std::cout << "Analyzing directions " << std::endl;
1543 
1544  double w = 0.0;
1545  double wH = 0.0;
1546  int volsize = ZSIZE(VRiesz);
1547 
1548  //Checking with MonoRes at 50A;
1549  int aux_idx;
1550 
1551  if (maxRes>18)
1552  {
1553  DIGFREQ2FFT_IDX(sampling/18, volsize, aux_idx);
1554 
1555  FFT_IDX2DIGFREQ(aux_idx, volsize, w);
1556  FFT_IDX2DIGFREQ(aux_idx+1, volsize, wH); //Frequency chosen for a first estimation
1557  }
1558  else
1559  {
1560  FFT_IDX2DIGFREQ(3, volsize, w);
1561  FFT_IDX2DIGFREQ(4, volsize, w);
1562  aux_idx = 3;
1563  }
1564  //std::cout << "fourier idx = " << aux_idx << std::endl;
1565  //std::cout << "Calling MonoRes core as a first estimation at " << sampling/w << "A." << std::endl;
1566 
1567  MultidimArray<double> amplitudeMS;
1568  double AvgNoise;
1569  AvgNoise = firstMonoResEstimation(fftV, w, wH, amplitudeMS)/9.0;
1570 
1571  N_directions=angles.mdimx;
1572 
1573  std::cout << "N_directions = " << N_directions << std::endl;
1574 
1575  double cone_angle = 45.0; //(degrees)
1576  cone_angle = PI*cone_angle/180;
1577 
1578  trigProducts.initZeros(3, N_directions);
1579 
1580  Image<double> outputResolution;
1581 
1582  for (size_t dir=0; dir<N_directions; dir++)
1583  {
1584  outputResolution().initZeros(VRiesz);
1585 // MultidimArray<double> &pOutputResolution = outputResolution();
1586  double freq, freqL, freqH, counter, resolution_2;
1587  MultidimArray<int> mask_aux = mask();
1588  MultidimArray<int> &pMask = mask_aux;
1589  std::vector<double> list;
1590  double resolution; //A huge value for achieving last_resolution < resolution
1591 
1592  double max_meanS = -1e38;
1593  double cut_value = 0.025;
1594 
1595  bool doNextIteration=true;
1596 
1597  int fourier_idx, last_fourier_idx = -1, iter = 0, fourier_idx_2;
1598  fourier_idx = aux_idx;
1599  int count_res = 0;
1600  double rot = MAT_ELEM(angles, 0, dir);
1601  double tilt = MAT_ELEM(angles, 1, dir);
1602  MAT_ELEM(trigProducts, 0, dir) = sin(tilt*PI/180)*cos(rot*PI/180);
1603  MAT_ELEM(trigProducts, 1, dir) = sin(tilt*PI/180)*sin(rot*PI/180);
1604  MAT_ELEM(trigProducts, 2, dir) = cos(tilt*PI/180);
1605  std::cout << "--------------NEW DIRECTION--------------" << std::endl;
1606  std::cout << "direction = " << dir+1 << " rot = " << rot << " tilt = " << tilt << std::endl;
1607 
1608 
1609  std::vector<float> noiseValues;
1610  FileName fnDebug;
1611  double last_resolution = 0;
1612 
1613  defineCone(fftV, conefilter, rot, tilt);
1614  maskMatrix.initConstant(1, NVoxelsOriginalMask, 1);
1615  do
1616  {
1617  continueIter = false;
1618  breakIter = false;
1619  //std::cout << "--------------Frequency--------------" << std::endl;
1620 
1621  resolution2eval_(fourier_idx, step,
1622  resolution, last_resolution, last_fourier_idx,
1623  freq, freqL, freqH,
1624  continueIter, breakIter, doNextIteration);
1625 
1626  if (breakIter)
1627  break;
1628 
1629  if (continueIter)
1630  continue;
1631 
1632  list.push_back(resolution);
1633 
1634  if (iter<2)
1635  resolution_2 = list[0];
1636  else
1637  resolution_2 = list[iter - 2];
1638 
1639  fnDebug = "Signal";
1640 
1641  amplitudeMonogenicSignal3D_fast(conefilter, freq, freqH, freqL, amplitudeMS, iter, dir, fnDebug);
1642 
1643  double sumS=0, sumS2=0, sumN=0, sumN2=0, NN = 0, NS = 0;
1644  noiseValues.clear();
1645 
1646 
1647  double x_dir = sin(tilt*PI/180)*cos(rot*PI/180);
1648  double y_dir = sin(tilt*PI/180)*sin(rot*PI/180);
1649  double z_dir = cos(tilt*PI/180);
1650 
1651  double uz, uy, ux;
1652 
1653  int n=0;
1654  int z_size = ZSIZE(amplitudeMS);
1655  int x_size = XSIZE(amplitudeMS);
1656  int y_size = YSIZE(amplitudeMS);
1657 
1658  size_t idx_mask;
1659  idx_mask = 0;
1660 
1661  double amplitudeValue;
1662 
1663  for(int k=0; k<z_size; ++k)
1664  {
1665  for(int i=0; i<y_size; ++i)
1666  {
1667  for(int j=0; j<x_size; ++j)
1668  {
1669  if (DIRECT_MULTIDIM_ELEM(pMask, n)>=1)
1670  {
1671  if (MAT_ELEM(maskMatrix, 0, idx_mask) >0)
1672  {
1673  amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitudeMS, n);
1674  sumS += amplitudeValue;
1675  ++NS;
1676  }
1677  ++idx_mask;
1678 
1679  }
1680  else if (DIRECT_MULTIDIM_ELEM(pMask, n)==0)
1681  {
1682  uz = (k - z_size*0.5);
1683  ux = (j - x_size*0.5);
1684  uy = (i - y_size*0.5);
1685 
1686  double rad = sqrt(ux*ux + uy*uy + uz*uz);
1687  double iun = 1/rad;
1688 
1689  //BE CAREFULL with the order
1690  double dotproduct = (uy*y_dir + ux*x_dir + uz*z_dir)*iun;
1691 
1692  double acosine = acos(dotproduct);
1693 
1694  //TODO: change efficiency the if condition by means of fabs(cos(angle))
1695  if (((acosine<cone_angle) || (acosine>(PI-cone_angle)) )
1696  && (rad>Rparticle))
1697  {
1698 // DIRECT_MULTIDIM_ELEM(coneVol, n) = 1;
1699  amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitudeMS, n);
1700  noiseValues.push_back((float) amplitudeValue);
1701  sumN += amplitudeValue;
1702  sumN2 += amplitudeValue*amplitudeValue;
1703  ++NN;
1704  }
1705  }
1706  ++n;
1707  }
1708  }
1709  }
1710 
1711  #ifdef DEBUG_DIR
1712  if (iter == 0)
1713  {
1714  Image<double> img;
1715 
1716  FileName iternumber;
1717  iternumber = formatString("cone_noise_%i_%i.vol", dir, iter);
1718  img = coneVol;
1719  img.write(iternumber);
1720  }
1721  #endif
1722 
1723  if ( (NS/(double) NVoxelsOriginalMask)<cut_value ) //when the 2.5% is reached then the iterative process stops
1724  {
1725  std::cout << "Search of resolutions stopped due to mask has been completed" << std::endl;
1726  doNextIteration =false;
1727  Nvoxels = 0;
1728 
1729  #ifdef DEBUG_MASK
1730  mask.write("partial_mask.vol");
1731  #endif
1732  }
1733  else
1734  {
1735  if (NS == 0)
1736  {
1737  std::cout << "There are no points to compute inside the mask" << std::endl;
1738  std::cout << "If the number of computed frequencies is low, perhaps the provided"
1739  "mask is not enough tight to the volume, in that case please try another mask" << std::endl;
1740  break;
1741  }
1742 
1743  double meanS=sumS/NS;
1744  // double sigma2S=sumS2/NS-meanS*meanS;
1745  double meanN=sumN/NN;
1746  double sigma2N=sumN2/NN-meanN*meanN;
1747 
1748  if (meanS>max_meanS)
1749  max_meanS = meanS;
1750 
1751  if (meanS<0.001*AvgNoise)//0001*max_meanS)
1752  {
1753  //std::cout << " meanS= " << meanS << " sigma2S= " << sigma2S << " NS = " << NS << std::endl;
1754  //std::cout << " meanN= " << meanN << " sigma2N= " << sigma2N << " NN= " << NN << std::endl;
1755  std::cout << "Search of resolutions stopped due to too low signal" << std::endl;
1756  std::cout << "\n"<< std::endl;
1757  doNextIteration = false;
1758  }
1759  else
1760  {
1761  // Check local resolution
1762  double thresholdNoise;
1763  //thresholdNoise = meanN+criticalZ*sqrt(sigma2N);
1764 
1765  std::sort(noiseValues.begin(),noiseValues.end());
1766  thresholdNoise = (double) noiseValues[size_t(noiseValues.size()*significance)];
1767 
1768  //std::cout << "thr="<< thresholdNoise << " " << meanN+criticalZ*sqrt(sigma2N) << " " << NN << std::endl;
1769  noiseValues.clear();
1770 
1771  std::cout << "Iteration = " << iter << ", Resolution= " << resolution << ", Signal = " << meanS << ", Noise = " << meanN << ", Threshold = " << thresholdNoise <<std::endl;
1772 
1773 
1774  size_t maskPos = 0;
1775  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitudeMS)
1776  {
1777  if (DIRECT_MULTIDIM_ELEM(pMask, n)>=1)
1778  {
1779  if (MAT_ELEM(maskMatrix, 0, maskPos) >=1)
1780  {
1781  if (DIRECT_MULTIDIM_ELEM(amplitudeMS, n)>thresholdNoise)
1782  {
1783  MAT_ELEM(resolutionMatrix, dir, maskPos) = resolution;
1784  MAT_ELEM(maskMatrix, 0, maskPos) = 1;
1785  }
1786  else
1787  {
1788  MAT_ELEM(maskMatrix, 0, maskPos) += 1;
1789  if (MAT_ELEM(maskMatrix, 0, maskPos) >2)
1790  {
1791  MAT_ELEM(maskMatrix, 0, maskPos) = 0;
1792  MAT_ELEM(resolutionMatrix, dir, maskPos) = resolution_2;
1793  }
1794  }
1795  }
1796  ++maskPos;
1797  }
1798  }
1799 
1800  if (doNextIteration)
1801  if (resolution <= (minRes-0.001))
1802  doNextIteration = false;
1803  }
1804  }
1805  ++iter;
1806  last_resolution = resolution;
1807  }while(doNextIteration);
1808 
1809 // amplitudeMS.clear();
1810 // fftVRiesz.clear();
1811 
1812 // size_t maskPos=0;
1813 // Image<double> ResolutionVol;
1814 // MultidimArray<double> &pResolutionVol = ResolutionVol();
1815 //
1816 // pResolutionVol.initZeros(amplitudeMS);
1817 // FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(pResolutionVol)
1818 // {
1819 // if (DIRECT_MULTIDIM_ELEM(mask(), n) == 1)
1820 // {
1821 // double myres = MAT_ELEM(resolutionMatrix, dir, maskPos);
1822 // DIRECT_MULTIDIM_ELEM(pResolutionVol, n) = myres;
1825 // ++maskPos;
1826 // }
1827 // }
1828 // #ifdef DEBUG_DIR
1829 // Image<double> saveImg;
1830 // saveImg = pResolutionVol;
1831 // FileName fnres = formatString("resolution_dir_%i.vol", dir+1);
1832 // saveImg.write(fnres);
1833 // saveImg.clear();
1834 // #endif
1835 // pResolutionVol.clear();
1836  list.clear();
1837 
1838  std::cout << "----------------direction-finished----------------" << std::endl;
1839  }
1840 
1841 
1842 
1843 
1845 
1846  int maskPos = 0;
1847 
1848  //Remove outliers
1849  removeOutliers(resolutionMatrix);
1850 
1851  MultidimArray<double> radial, azimuthal, lowestResolution, highestResolution, doavol1, doavol2;
1852  MetaDataVec prefDir;
1853 
1854  double radialThr, azimuthalThr;
1855  radialAzimuthalResolution(resolutionMatrix, mask(), radial, azimuthal,
1856  lowestResolution, highestResolution, doavol1, doavol2, radialThr, azimuthalThr, prefDir);
1857 
1858  Image<double> imgdoa;
1859  imgdoa = radial;
1860  imgdoa.write(fnradial);
1861  imgdoa = azimuthal;
1862  imgdoa.write(fnazimuthal);
1863  imgdoa = lowestResolution;
1864  imgdoa.write(fnLowestResolution);
1865  imgdoa = highestResolution;
1866  imgdoa.write(fnHighestResolution);
1867  imgdoa = doavol1;
1868  imgdoa.write(fnDoa1);
1869  imgdoa = doavol2;
1870  imgdoa.write(fnDoa2);
1871 
1872  MetaDataVec mdRadialAzimuthalThr;
1873  size_t objIdx;
1874  objIdx = mdRadialAzimuthalThr.addObject();
1875  mdRadialAzimuthalThr.setValue(MDL_RESOLUTION_FREQ, radialThr, objIdx);
1876  mdRadialAzimuthalThr.setValue(MDL_RESOLUTION_FREQ2, azimuthalThr, objIdx);
1877 
1878  mdRadialAzimuthalThr.write(fnMDThr);
1879 
1880  //std::cout << "radial = " << radialThr << " azimuthal = " << azimuthalThr << std::endl;
1881  std::cout << "Calculating the radial and azimuthal resolution " << std::endl;
1882 
1883 
1884  MetaDataVec mdRadial, mdAvg, mdHighest, mdLowest;
1885 
1886  Image<double> monores;
1887  monores.read(fnMonoRes);
1888  MultidimArray<double> monoresVol;
1889  monoresVol = monores();
1890  radialAverageInMask(mask(), radial, azimuthal, highestResolution, lowestResolution, monoresVol, mdAvg);
1891  mdAvg.write(fnRadialAvG);
1892 }
ProgResDir::transformer_inv
FourierTransformer transformer_inv
Definition: resolution_directional.h:129
radialAvg
void radialAvg(Image< double > &op)
Definition: image_operate.cpp:269
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
MDL_VOLUME_SCORE3
< Score 2 for volumes
Definition: metadata_label.h:478
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
FourierTransformer::inverseFourierTransform
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
resolution_directional.h
FourierTransformer::setThreadsNumber
void setThreadsNumber(int tNumber)
Definition: xmipp_fftw.h:119
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
c
doublereal * c
Definition: numerical_recipes.cpp:2230
Matrix2D::initConstant
void initConstant(T val)
Definition: matrix2d.h:602
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ProgResDir::defineCone
void defineCone(MultidimArray< std::complex< double > > &myfftV, MultidimArray< std::complex< double > > &conefilter, double rot, double tilt)
Definition: resolution_directional.cpp:545
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
y
static double * y
Definition: numerical_recipes.cpp:8487
ProgResDir::fftVRiesz
MultidimArray< std::complex< double > > fftVRiesz
Definition: resolution_directional.h:127
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
icdf_gauss
double icdf_gauss(double p)
Definition: xmipp_funcs.cpp:212
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
w
doublereal * w
Definition: numerical_recipes.cpp:2477
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
ProgResDir::radialAzimuthalResolution
void radialAzimuthalResolution(Matrix2D< double > &resolutionMat, MultidimArray< int > &pmask, MultidimArray< double > &radial, MultidimArray< double > &azimuthal, MultidimArray< double > &lowestResolution, MultidimArray< double > &highestResolution, MultidimArray< double > &doaResolution_1, MultidimArray< double > &doaResolution_2, double &radial_Thr, double &azimuthal_Thr, MetaDataVec &mdprefDirs)
Definition: resolution_directional.cpp:1182
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
DIGFREQ2FFT_IDX
#define DIGFREQ2FFT_IDX(freq, size, idx)
Definition: xmipp_fft.h:61
ProgResDir::angles
Matrix2D< double > angles
Definition: resolution_directional.h:131
ProgResDir::firstMonoResEstimation
double firstMonoResEstimation(MultidimArray< std::complex< double > > &myfftV, double w1, double w1l, MultidimArray< double > &amplitude)
Definition: resolution_directional.cpp:1365
ProgResDir::resolutionMatrix
Matrix2D< double > resolutionMatrix
Definition: resolution_directional.h:131
ProgResDir::ellipsoidFitting
void ellipsoidFitting(Matrix2D< double > &resolutionMat, Matrix2D< double > &axis)
Definition: resolution_directional.cpp:828
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
ProgResDir::radialAverageInMask
void radialAverageInMask(MultidimArray< int > &mask, MultidimArray< double > &inputVol_1, MultidimArray< double > &inputVol_2, MultidimArray< double > &inputVol_3, MultidimArray< double > &inputVol_4, MultidimArray< double > &inputVol_5, MetaDataVec &md)
Definition: resolution_directional.cpp:1010
ProgResDir::fnLowestResolution
FileName fnLowestResolution
Definition: resolution_directional.h:55
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
threshold
void threshold(double *phi, unsigned long nvox, double limit)
Definition: lib_vwk.cpp:524
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
ProgResDir::resolution2eval_
void resolution2eval_(int &fourier_idx, double min_step, double &resolution, double &last_resolution, int &last_fourier_idx, double &freq, double &freqL, double &freqH, bool &continueIter, bool &breakIter, bool &doNextIteration)
Definition: resolution_directional.cpp:624
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
b
doublereal * b
Definition: numerical_recipes.cpp:2230
MDL_COUNT
Number of elements of a type (int) [this is a genereic type do not use to transfer information to ano...
Definition: metadata_label.h:92
axis
char axis
Definition: image_operate.cpp:257
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
MDL_RESOLUTION_FREQ2
Frequency in 1/A squared (double)
Definition: metadata_label.h:389
ProgResDir::amplitudeMonogenicSignal3D_fast
void amplitudeMonogenicSignal3D_fast(const MultidimArray< std::complex< double > > &myfftV, double w1, double w1l, double wH, MultidimArray< double > &amplitude, int count, int dir, FileName fnDebug)
Definition: resolution_directional.cpp:347
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_VOLUME_SCORE4
< Score 3 for volumes
Definition: metadata_label.h:479
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
RAISED_COSINE
#define RAISED_COSINE
Definition: fourier_filter.h:72
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
generalizedEigs
void generalizedEigs(const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix1D< double > &D, Matrix2D< double > &P)
Definition: matrix2d.cpp:267
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
z
double z
Definition: numerical_recipes.cpp:8490
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
ProgResDir::iu
MultidimArray< double > iu
Definition: resolution_directional.h:128
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
ProgResDir::VRiesz
MultidimArray< double > VRiesz
Definition: resolution_directional.h:128
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
ProgResDir::freq_fourier
Matrix1D< double > freq_fourier
Definition: resolution_directional.h:132
j
#define j
Definition: numerical_recipes.cpp:2493
MDL_RESOLUTION_FREQ
Frequency in 1/A (double)
Definition: metadata_label.h:388
ProgResDir::diagSymMatrix3x3
void diagSymMatrix3x3(Matrix2D< double > A, Matrix1D< double > &eigenvalues, Matrix2D< double > &P)
Definition: resolution_directional.cpp:611
distance
TYPE distance(struct Point_T *p, struct Point_T *q)
Definition: point.cpp:28
ProgResDir::fnHighestResolution
FileName fnHighestResolution
Definition: resolution_directional.h:55
ProgResDir::trigProducts
Matrix2D< double > trigProducts
Definition: resolution_directional.h:131
alglib::round
int round(double x)
Definition: ap.cpp:7245
ProgResDir::conefilter
MultidimArray< std::complex< double > > conefilter
Definition: resolution_directional.h:130
MDL_VOLUME_SCORE1
Definition: metadata_label.h:476
formatString
String formatString(const char *format,...)
Definition: xmipp_strings.cpp:602
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
MDL_X
X component (double)
Definition: metadata_label.h:495
MDL_IDX
Index within a list (size_t)
Definition: metadata_label.h:245
ProgResDir::removeOutliers
void removeOutliers(Matrix2D< double > &resolutionMat)
Definition: resolution_directional.cpp:719
u
doublereal * u
Definition: numerical_recipes.cpp:2230
ProgResDir::fftVRiesz_aux
MultidimArray< std::complex< double > > fftVRiesz_aux
Definition: resolution_directional.h:127
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
MDL_VOLUME_SCORE2
< Score 1 for volumes
Definition: metadata_label.h:477
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
MDL_AVG
average value (double)
Definition: metadata_label.h:70
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
ProgResDir::maskMatrix
Matrix2D< double > maskMatrix
Definition: resolution_directional.h:131
Matrix1D::sort
Matrix1D< T > sort() const
Definition: matrix1d.cpp:850
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
ProgResDir::fftV
MultidimArray< std::complex< double > > fftV
Definition: resolution_directional.h:130
PI
#define PI
Definition: tools.h:43
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331
Matrix1D::initConstant
void initConstant(T val)
Definition: matrix1d.cpp:83
n
int * n
Definition: numerical_recipes.cpp:2229
a
doublereal * a
Definition: numerical_recipes.cpp:2230
Image::clear
void clear()
Definition: xmipp_image.h:144
LOWPASS
#define LOWPASS
Definition: fourier_filter.h:77
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305
ProgResDir::generateGridProjectionMatching
void generateGridProjectionMatching(Matrix2D< double > &angles)
Definition: resolution_directional.cpp:205
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710
MDL_WEIGHT
< Score 4 for volumes
Definition: metadata_label.h:480
Generated by   doxygen 1.8.13