Xmipp  v3.23.11-Nereus
sampling.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Roberto Marabini
4  *
5  * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20  * 02111-1307 USA
21  *
22  * All comments concerning this program package may be sent to the
23  * e-mail address 'xmipp@cnb.csic.es'
24  ***************************************************************************/
25 #include <fstream>
26 #include "sampling.h"
27 #include "core/geometry.h"
28 #include "core/xmipp_image_macros.h"
29 #include "core/metadata_vec.h"
30 
31 /* Default Constructor */
32 Sampling::Sampling()
33 {
34  Vertex aux;
35  aux.rot = -PI / 5.;
36  aux.tilt = PI / 2. - cte_w ;
37  aux.psi = 0.;
38  vertices_angles.push_back(aux);
39  aux.rot = PI / 5.;
40  aux.tilt = PI / 2. - cte_w ;
41  aux.psi = 0.;
42  vertices_angles.push_back(aux);
43  aux.rot = 3. * PI / 5.;
44  aux.tilt = PI / 2. - cte_w ;
45  aux.psi = 0.;
46  vertices_angles.push_back(aux);
47  aux.rot = 5. * PI / 5.;
48  aux.tilt = PI / 2. - cte_w ;
49  aux.psi = 0.;
50  vertices_angles.push_back(aux);
51  aux.rot = -3. * PI / 5.;
52  aux.tilt = PI / 2. - cte_w ;
53  aux.psi = 0.;
54  vertices_angles.push_back(aux);
55  aux.rot = 0. / 5.;
56  aux.tilt = -cte_w + PI / 2. ;
57  aux.psi = 0.;
58  vertices_angles.push_back(aux);
59  aux.rot = 2. * PI / 5.;
60  aux.tilt = -cte_w + PI / 2. ;
61  aux.psi = 0.;
62  vertices_angles.push_back(aux);
63  aux.rot = 4. * PI / 5.;
64  aux.tilt = -cte_w + PI / 2. ;
65  aux.psi = 0.;
66  vertices_angles.push_back(aux);
67  aux.rot = -4. * PI / 5.;
68  aux.tilt = -cte_w + PI / 2. ;
69  aux.psi = 0.;
70  vertices_angles.push_back(aux);
71  aux.rot = -2. * PI / 5.;
72  aux.tilt = -cte_w + PI / 2. ;
73  aux.psi = 0.;
74  vertices_angles.push_back(aux);
75 
76  vertices_vectors.push_back(vectorR3(0., 0., 1.));
77  vertices_vectors.push_back(vectorR3(0.723606900230461, -0.525731185781806, 0.447213343087301));
78  vertices_vectors.push_back(vectorR3(0.723606900230461, 0.525731185781806, 0.447213343087301));
79  vertices_vectors.push_back(vectorR3(-0.276393239417711, 0.850650928976665, 0.447213343087301));
80  vertices_vectors.push_back(vectorR3(-0.8944273172062, 0., 0.447213343087301));
81  vertices_vectors.push_back(vectorR3(-0.276393239417711, -0.850650928976665, 0.447213343087301));
82  vertices_vectors.push_back(vectorR3(0.8944273172062, 0., -0.447213343087301));
83  vertices_vectors.push_back(vectorR3(0.276393242471372, 0.850650927984471, -0.447213343087301));
84  vertices_vectors.push_back(vectorR3(-0.723606898343194, 0.525731188379405, -0.447213343087301));
85  vertices_vectors.push_back(vectorR3(-0.723606898343194, -0.525731188379405, -0.447213343087301));
86  vertices_vectors.push_back(vectorR3(0.276393242471372, -0.850650927984471, -0.447213343087301));
87  vertices_vectors.push_back(vectorR3(0., 0., -1.));
88 
89  sampling_noise=0.0;
90  cos_neighborhood_radius=-1.01;
91 
92  numberSamplesAsymmetricUnit=-1;
93  exp_data_fileNames.clear();
94 
95  verbose=1;
96  //#define DEBUG1
97 #ifdef DEBUG1
98 
99  for (int i = 0;
100  i < vertices_vectors.size();
101  i++)
102  std::cout << vertices_vectors[].transpose() << std::endl;
103 #endif
104 #undef DEBUG1
105 }
106 
107 bool Sampling::operator==(const Sampling& op) const
108 {
109 
110  return (XMIPP_EQUAL_REAL(sampling_rate_rad, op.sampling_rate_rad) &&
111  XMIPP_EQUAL_REAL(cos_neighborhood_radius, op.cos_neighborhood_radius) &&
112  no_redundant_sampling_points_angles == op.no_redundant_sampling_points_angles &&
113  no_redundant_sampling_points_vector == op.no_redundant_sampling_points_vector &&
114  no_redundant_sampling_points_index == op.no_redundant_sampling_points_index &&
115  my_neighbors == op.my_neighbors //&&
116  //exp_data_fileNames == op.exp_data_fileNames
117  );
118 
119 }
120 
121 void Sampling::setSampling(double sampling)
122 {
123  sampling_rate_rad = DEG2RAD(sampling);
124  number_of_samples = ROUND(cte_w / sampling_rate_rad)+1;
125  if (number_of_samples < 3)
126  {
127  std::cerr << "maximum value of angular sampling rate is "
128  << cte_w*0.5*180./PI
129  << std::endl;
130  exit(1);
131  }
132 }
133 
134 void Sampling::setNoise(double noise_deviation, int my_seed)
135 {
136  sampling_noise = noise_deviation;
137  init_random_generator(my_seed);
138 }
139 
140 void Sampling::setNeighborhoodRadius(double neighborhood)
141 {
142  if(neighborhood<0)
143  cos_neighborhood_radius=-1.01;
144  else if(neighborhood>180.001)
145  REPORT_ERROR(ERR_ARG_INCORRECT,"Neighborhood can not be greater than 180. \n \
146  Use any negative value to cover the whole space ");
147  else
148  {
149  neighborhood_radius_rad = DEG2RAD(neighborhood);
150  cos_neighborhood_radius = cos(neighborhood_radius_rad);
151  }
152 }
153 
154 /* Compute edge sampling points using Baumgardner 1995 */
155 void Sampling::computeSamplingPoints(bool only_half_sphere,
156  double max_tilt,
157  double min_tilt)
158 {
160  std::vector <Matrix1D<double> > edge_vector_start;
162  std::vector <Matrix1D<double> > edge_vector_end;
163  // I need 10 auxiliary vector for edges
164  Matrix1D<double> starting_point, ending_point;
165  double max_z;
166  double min_z;
167  sampling_points_angles.clear();
168  sampling_points_vector.clear();
169 
170  max_z=cos(PI * max_tilt / 180.);
171  min_z=cos(PI * min_tilt / 180.);
172  if ( (min_tilt > max_tilt) ||
173  (min_tilt <0 ) ||
174  (max_tilt < 0) ||
175  (max_tilt > 180.)
176  )
177  {//OK
178  std::cerr << "tilt angles cannot be negative or min_tilt > max_tilt" << std::endl;
179  std::cerr << "min_tilt=" << min_tilt << "max_tilt=" << max_tilt << std::endl;
180  exit(0);
181  }
182 
183  if (min_z>max_z)
184  {
185  double aux=min_z;
186  min_z=max_z;
187  max_z=aux;
188  }
189 
190  //01a
191  starting_point = vertices_vectors[0];
192  ending_point = vertices_vectors[1];
193  fillEdge(starting_point, ending_point, edge_vector_start, false);
194  starting_point = vertices_vectors[6];
195  ending_point = vertices_vectors[1];
196  fillEdge(starting_point, ending_point, edge_vector_start, true);
197  //01b
198  starting_point = vertices_vectors[0];
199  ending_point = vertices_vectors[2];
200  fillEdge(starting_point, ending_point, edge_vector_end, false);
201  starting_point = vertices_vectors[6];
202  ending_point = vertices_vectors[2];
203  fillEdge(starting_point, ending_point, edge_vector_end, true);
204  //02a
205  starting_point = vertices_vectors[0];
206  ending_point = vertices_vectors[2];
207  fillEdge(starting_point, ending_point, edge_vector_start, false);
208  starting_point = vertices_vectors[7];
209  ending_point = vertices_vectors[2];
210  fillEdge(starting_point, ending_point, edge_vector_start, true);
211  //02b
212  starting_point = vertices_vectors[0];
213  ending_point = vertices_vectors[3];
214  fillEdge(starting_point, ending_point, edge_vector_end, false);
215  starting_point = vertices_vectors[7];
216  ending_point = vertices_vectors[3];
217  fillEdge(starting_point, ending_point, edge_vector_end, true);
218 
219  //03a
220  starting_point = vertices_vectors[0];
221  ending_point = vertices_vectors[3];
222  fillEdge(starting_point, ending_point, edge_vector_start, false);
223  starting_point = vertices_vectors[8];
224  ending_point = vertices_vectors[3];
225  fillEdge(starting_point, ending_point, edge_vector_start, true);
226  //03b
227  starting_point = vertices_vectors[0];
228  ending_point = vertices_vectors[4];
229  fillEdge(starting_point, ending_point, edge_vector_end, false);
230  starting_point = vertices_vectors[8];
231  ending_point = vertices_vectors[4];
232  fillEdge(starting_point, ending_point, edge_vector_end, true);
233 
234  //04a
235  starting_point = vertices_vectors[0];
236  ending_point = vertices_vectors[4];
237  fillEdge(starting_point, ending_point, edge_vector_start, false);
238  starting_point = vertices_vectors[9];
239  ending_point = vertices_vectors[4];
240  fillEdge(starting_point, ending_point, edge_vector_start, true);
241  //04b
242  starting_point = vertices_vectors[0];
243  ending_point = vertices_vectors[5];
244  fillEdge(starting_point, ending_point, edge_vector_end, false);
245  starting_point = vertices_vectors[9];
246  ending_point = vertices_vectors[5];
247  fillEdge(starting_point, ending_point, edge_vector_end, true);
248 
249  //05a
250  starting_point = vertices_vectors[0];
251  ending_point = vertices_vectors[5];
252  fillEdge(starting_point, ending_point, edge_vector_start, false);
253  starting_point = vertices_vectors[10];
254  ending_point = vertices_vectors[5];
255  fillEdge(starting_point, ending_point, edge_vector_start, true);
256  //05b
257  starting_point = vertices_vectors[0];
258  ending_point = vertices_vectors[1];
259  fillEdge(starting_point, ending_point, edge_vector_end, false);
260  starting_point = vertices_vectors[10];
261  ending_point = vertices_vectors[1];
262  fillEdge(starting_point, ending_point, edge_vector_end, true);
263 
264  //06a
265  starting_point = vertices_vectors[11];
266  ending_point = vertices_vectors[10];
267  fillEdge(starting_point, ending_point, edge_vector_start, false);
268  starting_point = vertices_vectors[5];
269  ending_point = vertices_vectors[10];
270  fillEdge(starting_point, ending_point, edge_vector_start, true);
271  //06b
272  starting_point = vertices_vectors[11];
273  ending_point = vertices_vectors[9];
274  fillEdge(starting_point, ending_point, edge_vector_end, false);
275  starting_point = vertices_vectors[5];
276  ending_point = vertices_vectors[9];
277  fillEdge(starting_point, ending_point, edge_vector_end, true);
278 
279  //07a
280  starting_point = vertices_vectors[11];
281  ending_point = vertices_vectors[9];
282  fillEdge(starting_point, ending_point, edge_vector_start, false);
283  starting_point = vertices_vectors[4];
284  ending_point = vertices_vectors[9];
285  fillEdge(starting_point, ending_point, edge_vector_start, true);
286  //07b
287  starting_point = vertices_vectors[11];
288  ending_point = vertices_vectors[8];
289  fillEdge(starting_point, ending_point, edge_vector_end, false);
290  starting_point = vertices_vectors[4];
291  ending_point = vertices_vectors[8];
292  fillEdge(starting_point, ending_point, edge_vector_end, true);
293 
294  //08a
295  starting_point = vertices_vectors[11];
296  ending_point = vertices_vectors[8];
297  fillEdge(starting_point, ending_point, edge_vector_start, false);
298  starting_point = vertices_vectors[3];
299  ending_point = vertices_vectors[8];
300  fillEdge(starting_point, ending_point, edge_vector_start, true);
301  //08b
302  starting_point = vertices_vectors[11];
303  ending_point = vertices_vectors[7];
304  fillEdge(starting_point, ending_point, edge_vector_end, false);
305  starting_point = vertices_vectors[3];
306  ending_point = vertices_vectors[7];
307  fillEdge(starting_point, ending_point, edge_vector_end, true);
308 
309  //09a
310  starting_point = vertices_vectors[11];
311  ending_point = vertices_vectors[7];
312  fillEdge(starting_point, ending_point, edge_vector_start, false);
313  starting_point = vertices_vectors[2];
314  ending_point = vertices_vectors[7];
315  fillEdge(starting_point, ending_point, edge_vector_start, true);
316  //09b
317  starting_point = vertices_vectors[11];
318  ending_point = vertices_vectors[6];
319  fillEdge(starting_point, ending_point, edge_vector_end, false);
320  starting_point = vertices_vectors[2];
321  ending_point = vertices_vectors[6];
322  fillEdge(starting_point, ending_point, edge_vector_end, true);
323 
324  //10a
325  starting_point = vertices_vectors[11];
326  ending_point = vertices_vectors[6];
327  fillEdge(starting_point, ending_point, edge_vector_start, false);
328  starting_point = vertices_vectors[1];
329  ending_point = vertices_vectors[6];
330  fillEdge(starting_point, ending_point, edge_vector_start, true);
331  //10b
332  starting_point = vertices_vectors[11];
333  ending_point = vertices_vectors[10];
334  fillEdge(starting_point, ending_point, edge_vector_end, false);
335  starting_point = vertices_vectors[1];
336  ending_point = vertices_vectors[10];
337  fillEdge(starting_point, ending_point, edge_vector_end, true);
338 
339  //#define DEBUG2
340 #ifdef DEBUG2
341 
342  std::ofstream filestr1;
343  filestr1.open ("computeSamplingPoints_1.bild");
344 
345  for (int i = 0;
346  i < edge_vector_start.size();
347  i++)
348  {
349  filestr1 << ".sphere "
350  << edge_vector_start[i](0) << " "
351  << edge_vector_start[i](1) << " "
352  << edge_vector_start[i](2) << " "
353  << " .025" << std::endl;
354  filestr1 << ".sphere "
355  << edge_vector_end[i](0) << " "
356  << edge_vector_end[i](1) << " "
357  << edge_vector_end[i](2) << " "
358  << " .025" << std::endl;
359  }
360  filestr1.close();
361 
362  //std::cout << ending_point.transpose() << " 1.1 1.5 " << std::endl;
363 #endif
364 #undef DEBUG2
365  // add main corners that are not part of the basic octahedrons
366  {
367  int i=11;
368  if ( (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
369  || ZZ(vertices_vectors[i]) < min_z
370  || ZZ(vertices_vectors[i]) > max_z
371  )
372  //if (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
373  ;
374  else
375  sampling_points_vector.push_back(vertices_vectors[i]);
376  }
377  {
378  int i=0;
379  if ( (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
380  || ZZ(vertices_vectors[i]) < min_z
381  || ZZ(vertices_vectors[i]) > max_z
382  )
383  //if (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
384  ;
385  else
386  sampling_points_vector.push_back(vertices_vectors[i]);
387  }
388 
389  // add edges
390  for (size_t i = 0;
391  i < edge_vector_start.size();
392  i++)
393  {
394  if (i < number_of_samples * 10 - 15)
395  {
396  if ( (only_half_sphere && ZZ(edge_vector_start[i]) < 0.0)
397  || ZZ(edge_vector_start[i]) < min_z
398  || ZZ(edge_vector_start[i]) > max_z
399  )
400  //if (only_half_sphere && ZZ(edge_vector_start[i]) < 0.0)
401  continue;
402  else
403  sampling_points_vector.push_back(edge_vector_start[i]);
404  }
405  else
406  {
407  if ( (only_half_sphere && ZZ(edge_vector_end[i]) < 0.0)
408  || ZZ(edge_vector_end[i]) < min_z
409  || ZZ(edge_vector_end[i]) > max_z
410  )
411  //if (only_half_sphere && ZZ(edge_vector_end[i]) < 0.0)
412  continue;
413  else
414  sampling_points_vector.push_back(edge_vector_end[i]);
415  }
416  }
417  //#define DEBUG3
418 #ifdef DEBUG3
419  std::ofstream filestr;
420  filestr.open ("computeSamplingPoints_2.bild");
421  for (int i = 0;
422  i < sampling_points_vector.size();
423  i++)
424  {
425  filestr << ".color 1 0 0" << std::endl
426  << ".sphere "
427  << sampling_points_vector[i](0) << " "
428  << sampling_points_vector[i](1) << " "
429  << sampling_points_vector[i](2) << " "
430  << " .025" << std::endl;
431  //sampling_points_vector[i];//.add(0.0, sampling_noise, "gaussian");
432  sampling_points_vector[i].selfNormalize();
433  filestr << ".color 0 0 1" << std::endl
434  << ".sphere "
435  << sampling_points_vector[i](0) << " "
436  << sampling_points_vector[i](1) << " "
437  << sampling_points_vector[i](2) << " "
438  <<" .027" << std::endl;
439  }
440  filestr.close();
441 
442 #endif
443 #undef DEBUG3
444 
445  // add in between points
446  int j = 0;
447  bool j_flag = false;
448  for (size_t i = 0;
449  i < edge_vector_start.size();
450  i++)
451  {
452  if ((j % (number_of_samples - 1)) == 0 && j != 0)
453  {
454  j = 0;
455  j_flag = true;
456  }
457  if ((j % (number_of_samples - 2)) == 0 && j != 0 && j_flag == true)
458  {
459  j = 0;
460  j_flag = false;
461  }
462  fillDistance(edge_vector_start[i],
463  edge_vector_end[i],
464  sampling_points_vector,
465  (j + 1) % number_of_samples,
466  only_half_sphere,min_z,max_z);
467  j++;
468  }
469  //#define DEBUG3
470 #ifdef DEBUG3
471  std::ofstream filestr3;
472  filestr3.open ("computeSamplingPoints_3.bild");
473  for (int i = 0;
474  i < sampling_points_vector.size();
475  i++)
476  {
477  filestr3 << ".color 1 0 0" << std::endl
478  << ".sphere "
479  << sampling_points_vector[i](0) << " "
480  << sampling_points_vector[i](1) << " "
481  << sampling_points_vector[i](2) << " "
482  << " .025" << std::endl;
483  }
484  filestr3.close();
485 #endif
486 #undef DEBUG3
487 
488  //noisify angles
489  if(sampling_noise!=0.0)
490  {
491  for (size_t n = 0;
492  n < sampling_points_vector.size();
493  n++)
494  {
495  FOR_ALL_ELEMENTS_IN_MATRIX1D(sampling_points_vector[n])
496  (sampling_points_vector[n])(i) += rnd_gaus(0., sampling_noise);
497  sampling_points_vector[n].selfNormalize();
498  }
499  }
500 
501  //#define DEBUG3
502 #ifdef DEBUG3
503  for (int i = 0;
504  i < sampling_points_vector.size();
505  i++)
506  {
507  std::cout << ".color 0 1 0" << std::endl
508  << ".sphere " << sampling_points_vector[i].transpose() <<
509  " .03" << std::endl;
510  }
511 #endif
512 #undef DEBUG3
513 
514  // store sampling points as angles
515  Matrix1D<double> aux(3), aux1(3);
516  ZZ(aux) = 0.;
517  for (size_t i = 0;
518  i < sampling_points_vector.size();
519  i++)
520  {
521  XX(aux) = atan2(YY(sampling_points_vector[i]),
522  XX(sampling_points_vector[i]));
523  YY(aux) = acos(ZZ(sampling_points_vector[i]));
524  if (YY(aux) < 0.)
525  YY(aux) += PI;
526  sampling_points_angles.push_back(aux*180. / PI);
527  }
528 #ifdef NEVERDEFINED
529  //sort points
530  int k;
531  int bound = sampling_points_angles.size() - 1;
532  int t;
533  int last_swap;
534  Matrix1D<double> aux_angle(3);
535  Matrix1D<double> aux_vector(3);
536 
537  while (bound)
538  {
539  last_swap = 0;
540  for (k = 0; k < bound; k++)
541  {
542  aux = sampling_points_angles[k]; /* aux is a maximum of A[0]..A[k] */
543  aux1 = sampling_points_vector[k]; /* aux is a maximum of A[0]..A[k] */
544  if (sortFunc(aux, sampling_points_angles[k+1]))
545  {
546  sampling_points_angles[k] = sampling_points_angles[k+1];
547  // sampling_points_angles[k] = sampling_points_angles[k+1];
548  sampling_points_angles[k+1] = aux; /*swap*/
549  sampling_points_vector[k] = sampling_points_vector[k+1];
550  sampling_points_vector[k] = sampling_points_vector[k+1];
551  sampling_points_vector[k+1] = aux1; /*swap*/
552  last_swap = k; /* mark the last swap position */
553  }//if
554  }//for
555  bound = last_swap; /*elements after bound already sorted */
556  }//while
557 
558 
559 #endif
560  //#define DEBUG3
561 #ifdef DEBUG3
562  std::ofstream filestr4;
563  filestr4.open ("computeSamplingPoints_4.bild");
564 
565  filestr4 << ".color yellow"
566  << std::endl
567  ;
568  for (int i = 0;
569  i < sampling_points_vector.size();
570  i++)
571  {
572  filestr4 << ".sphere "
573  << sampling_points_vector[i](0) << " "
574  << sampling_points_vector[i](1) << " "
575  << sampling_points_vector[i](2) << " "
576  << " .025" << std::endl;
577  }
578  filestr4.close();
579 #endif
580 #undef DEBUG3
581 
582  numberSamplesAsymmetricUnit=sampling_points_vector.size();
583 }
584 
585 // return 1 if a should go first 0 is equal -1 if before
586 int Sampling::sortFunc(const Matrix1D<double> &t, const Matrix1D<double> &a)
587 {
588  if (YY(t) - 0.000001 > YY(a))
589  {
590  return 1;
591  }
592  else if (YY(t) + 0.000001 < YY(a))
593  {
594  return 0;
595  };
596  if (XX(t) > XX(a))
597  {
598  return 1;
599  }
600  else
601  {
602  return 0;
603  };
604 }
605 
606 void Sampling::fillEdge(const Matrix1D<double> &starting_point,
607  const Matrix1D<double> &ending_point,
608  std::vector <Matrix1D<double> > & edge_vector,
609  bool END_FLAG
610  )
611 {
612  Matrix1D<double> v_aux(3);
613 
614  double alpha;
615  double beta;
616  double gamma;
617  // skip first corener, already computed;
618  double upsilon = acos(dotProduct(starting_point, ending_point));
619  for (size_t i1 = 1; i1 < number_of_samples; i1++)
620  {
621  gamma = (double)i1 / (number_of_samples - 1);
622  alpha = sin((1. - gamma) * upsilon) / (sin(upsilon));
623  beta = sin(gamma * upsilon) / sin(upsilon);
624  v_aux = alpha * starting_point + beta * ending_point;
625  v_aux.selfNormalize();
626  if (beta > 0.9999 && END_FLAG)
627  continue;
628  edge_vector.push_back(v_aux);
629  }
630 }
631 void Sampling::fillDistance(const Matrix1D<double> &starting_point,
632  const Matrix1D<double> &ending_point,
633  std::vector <Matrix1D<double> > &
634  sampling_points_vector,
635  int my_number_of_samples,
636  bool only_half_sphere,
637  double min_z,
638  double max_z
639  )
640 {
641  Matrix1D<double> v_aux(3);
642 
643  double alpha;
644  double beta;
645  double gamma;
646  // skip first corener, already computed;
647  double upsilon = acos(dotProduct(starting_point, ending_point));
648 
649  for (int i1 = 1; i1 < my_number_of_samples; i1++)
650  {
651  gamma = (double)i1 / my_number_of_samples;
652  alpha = sin((1. - gamma) * upsilon) / (sin(upsilon));
653  beta = sin(gamma * upsilon) / sin(upsilon);
654  v_aux = alpha * starting_point + beta * ending_point;
655  v_aux.selfNormalize();
656  if ( (only_half_sphere && ZZ(v_aux) < 0.0)
657  || ZZ(v_aux) < min_z
658  || ZZ(v_aux) > max_z
659  )
660  //if (only_half_sphere && ZZ(v_aux) < 0.0)
661  continue;
662  else
663  sampling_points_vector.push_back(v_aux);
664  }
665  /*
666  //Remove points not in tilt range, check only the edges
667  for (int i = 0;
668  i < auxCounter;
669  i++)
670  {
671  if ( ZZ(sampling_points_vector[i]) < min_z
672  || ZZ(sampling_points_vector[i]) > max_z
673  )
674  sampling_points_vector[i].remove();
675  }
676  */
677 }
678 
679 #define CLEAR_VECTORS() \
680 no_redundant_sampling_points_vector.clear();\
681 no_redundant_sampling_points_angles.clear();\
682 no_redundant_sampling_points_index.clear()
683 
684 #define CREATE_INDEXES() \
685  size_t __size = no_redundant_sampling_points_angles.size();\
686  no_redundant_sampling_points_index.resize(__size, 0);\
687  size_t * __ptrIndex = &(no_redundant_sampling_points_index[0]);\
688  for (size_t i = 1; i < __size; ++i)\
689  __ptrIndex[i] = i
690 
691 void Sampling::removeRedundantPoints(const int symmetry, int sym_order)
692 {
693  Matrix2D<double> L(4, 4), R(4, 4);
694  Matrix2D<double> aux(3, 3);
695  Matrix1D<double> row1(3), row2(3);
696 
697  //int j_end=0;
698  Matrix1D<double> row(3);
699 
700  CLEAR_VECTORS();
701 
702  if (symmetry == pg_CN)
703  {//OK
704  for (size_t i = 0; i < sampling_points_angles.size(); i++)
705  {
706  if (XX(sampling_points_angles[i]) >= (-180. / sym_order) &&
707  XX(sampling_points_angles[i]) <= (180. / sym_order))
708  {
709  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
710  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
711  }
712  }// for i
713  }
714  else if (symmetry == pg_CI ||
715  symmetry == pg_CS )
716  {//OK
717  for (size_t i = 0; i < sampling_points_angles.size(); i++)
718  {
719  if (YY(sampling_points_angles[i]) <= 90)
720  {
721  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
722  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
723  }
724  }// for i
725  }
726  else if (symmetry == pg_CNV )
727  {//OK
728  for (size_t i = 0; i < sampling_points_angles.size(); i++)
729  {
730  if (XX(sampling_points_angles[i]) >= 0. / sym_order &&
731  XX(sampling_points_angles[i]) <= 180. / sym_order)
732  {
733  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
734  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
735  }
736  }// for i
737  }
738  else if (symmetry == pg_CNH )
739  {//OK
740  for (size_t i = 0; i < sampling_points_angles.size(); i++)
741  {
742  if (XX(sampling_points_angles[i]) >= -180. / sym_order &&
743  XX(sampling_points_angles[i]) <= 180. / sym_order &&
744  YY(sampling_points_angles[i]) <= 90.
745  )
746  {
747  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
748  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
749  }
750  }// for i
751  }
752  else if (symmetry == pg_SN )
753  {//OK
754  for (size_t i = 0; i < sampling_points_angles.size(); i++)
755  {
756  if (XX(sampling_points_angles[i]) >= -180.*2. / sym_order &&
757  XX(sampling_points_angles[i]) <= 180.*2. / sym_order &&
758  YY(sampling_points_angles[i]) <= 90.
759  )
760  {
761  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
762  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
763  }
764  }// for i
765  }
766  else if (symmetry == pg_DN )
767  {
768  for (size_t i = 0; i < sampling_points_angles.size(); i++)
769  {
770  if (XX(sampling_points_angles[i]) >= -180. / sym_order + 90. &&
771  XX(sampling_points_angles[i]) <= 180. / sym_order + 90. &&
772  YY(sampling_points_angles[i]) <= 90.
773  )
774  {
775  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
776  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
777  }
778  }// for i
779  }
780  else if (symmetry == pg_DNV )
781  {
782  for (size_t i = 0; i < sampling_points_angles.size(); i++)
783  {
784  if (XX(sampling_points_angles[i]) >= 0. + 90. &&
785  XX(sampling_points_angles[i]) <= 180. / sym_order +90. &&
786  YY(sampling_points_angles[i]) <= 90.
787  )
788  {
789  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
790  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
791  }
792  }// for i
793  }
794  else if (symmetry == pg_DNH )
795  {
796  for (size_t i = 0; i < sampling_points_angles.size(); i++)
797  {
798  if (XX(sampling_points_angles[i]) >= 0. &&
799  XX(sampling_points_angles[i]) <= 180. / sym_order &&
800  YY(sampling_points_angles[i]) <= 90.
801  )
802  {
803  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
804  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
805  }
806  }// for i
807  }
808  else if (symmetry == pg_T )
809  {//OK
810  Matrix1D<double> _3_fold_axis_1_by_3_fold_axis_2(3);
811  _3_fold_axis_1_by_3_fold_axis_2 = vectorR3(-0.942809, 0., 0.);
812  _3_fold_axis_1_by_3_fold_axis_2.selfNormalize();
813  Matrix1D<double> _3_fold_axis_2_by_3_fold_axis_3(3);
814  _3_fold_axis_2_by_3_fold_axis_3 = vectorR3(0.471405, 0.272165, 0.7698);
815  _3_fold_axis_2_by_3_fold_axis_3.selfNormalize();
816  Matrix1D<double> _3_fold_axis_3_by_3_fold_axis_1(3);
817  _3_fold_axis_3_by_3_fold_axis_1 = vectorR3(0.471404, 0.816497, 0.);
818  _3_fold_axis_3_by_3_fold_axis_1.selfNormalize();
819  for (size_t i = 0; i < sampling_points_angles.size(); i++)
820  {
821  if ((XX(sampling_points_angles[i]) >= 90. &&
822  XX(sampling_points_angles[i]) <= 150. )||
823  XX(sampling_points_angles[i]) == 0
824  )
825  if (
826  dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_3_fold_axis_2) >= 0 &&
827  dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_3_fold_axis_3) >= 0 &&
828  dotProduct(sampling_points_vector[i], _3_fold_axis_3_by_3_fold_axis_1) >= 0
829  )
830  {
831  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
832  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
833  }
834  }// for i
835  }
836  else if (symmetry == pg_TD )
837  {//OK
838  Matrix1D<double> _2_fold_axis_1_by_3_fold_axis_2(3);
839  _2_fold_axis_1_by_3_fold_axis_2 = vectorR3(-0.942809, 0., 0.);
840  _2_fold_axis_1_by_3_fold_axis_2.selfNormalize();
841  Matrix1D<double> _3_fold_axis_2_by_3_fold_axis_5(3);
842  _3_fold_axis_2_by_3_fold_axis_5 = vectorR3(0.471405, 0.272165, 0.7698);
843  _3_fold_axis_2_by_3_fold_axis_5.selfNormalize();
844  Matrix1D<double> _3_fold_axis_5_by_2_fold_axis_1(3);
845  _3_fold_axis_5_by_2_fold_axis_1 = vectorR3(0., 0.471405, -0.666667);
846  _3_fold_axis_5_by_2_fold_axis_1.selfNormalize();
847  for (size_t i = 0; i < sampling_points_angles.size(); i++)
848  {
849  // if ( XX(sampling_points_angles[i])>= 120. &&
850  // XX(sampling_points_angles[i])<= 150. ||
851  // XX(sampling_points_angles[i])== 0
852  // )
853  if (
854  dotProduct(sampling_points_vector[i], _2_fold_axis_1_by_3_fold_axis_2) >= 0 &&
855  dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_3_fold_axis_5) >= 0 &&
856  dotProduct(sampling_points_vector[i], _3_fold_axis_5_by_2_fold_axis_1) >= 0
857  )
858  {
859  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
860  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
861  }
862  }// for i
863  }
864  else if (symmetry == pg_TH )
865  {//OK
866  Matrix1D<double> _3_fold_axis_1_by_2_fold_axis_1(3);
867  _3_fold_axis_1_by_2_fold_axis_1 = vectorR3(-0.816496, 0., 0.);
868  _3_fold_axis_1_by_2_fold_axis_1.selfNormalize();
869  Matrix1D<double> _2_fold_axis_1_by_2_fold_axis_2(3);
870  _2_fold_axis_1_by_2_fold_axis_2 = vectorR3(0.707107, 0.408248, -0.57735);
871  _2_fold_axis_1_by_2_fold_axis_2.selfNormalize();
872  Matrix1D<double> _2_fold_axis_2_by_3_fold_axis_1(3);
873  _2_fold_axis_2_by_3_fold_axis_1 = vectorR3(-0.408248, -0.707107, 0.);
874  _2_fold_axis_2_by_3_fold_axis_1.selfNormalize();
875  for (size_t i = 0; i < sampling_points_angles.size(); i++)
876  {
877  // if ( XX(sampling_points_angles[i])>= 120. &&
878  // XX(sampling_points_angles[i])<= 150. ||
879  // XX(sampling_points_angles[i])== 0
880  // )
881  if (
882  dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_2_fold_axis_1) >= 0 &&
883  dotProduct(sampling_points_vector[i], _2_fold_axis_1_by_2_fold_axis_2) >= 0 &&
884  dotProduct(sampling_points_vector[i], _2_fold_axis_2_by_3_fold_axis_1) >= 0
885  )
886  {
887  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
888  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
889  }
890  }// for i
891  }
892  else if (symmetry == pg_O )
893  {//OK
894  Matrix1D<double> _3_fold_axis_1_by_3_fold_axis_2(3);
895  _3_fold_axis_1_by_3_fold_axis_2 = vectorR3(0., -1., 1.);
896  _3_fold_axis_1_by_3_fold_axis_2.selfNormalize();
897  Matrix1D<double> _3_fold_axis_2_by_4_fold_axis(3);
898  _3_fold_axis_2_by_4_fold_axis = vectorR3(1., 1., 0.);
899  _3_fold_axis_2_by_4_fold_axis.selfNormalize();
900  Matrix1D<double> _4_fold_axis_by_3_fold_axis_1(3);
901  _4_fold_axis_by_3_fold_axis_1 = vectorR3(-1., 1., 0.);
902  _4_fold_axis_by_3_fold_axis_1.selfNormalize();
903  for (size_t i = 0; i < sampling_points_angles.size(); i++)
904  {
905  if ((XX(sampling_points_angles[i]) >= 45. &&
906  XX(sampling_points_angles[i]) <= 135. &&
907  YY(sampling_points_angles[i]) <= 90.) ||
908  XX(sampling_points_angles[i]) == 0.
909  )
910  if (
911  dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_3_fold_axis_2) >= 0 &&
912  dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_4_fold_axis) >= 0 &&
913  dotProduct(sampling_points_vector[i], _4_fold_axis_by_3_fold_axis_1) >= 0
914  )
915  {
916  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
917  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
918  }
919  }// for i
920  }
921  else if (symmetry == pg_OH )
922  {//OK
923  Matrix1D<double> _3_fold_axis_1_by_3_fold_axis_2(3);
924  _3_fold_axis_1_by_3_fold_axis_2 = vectorR3(0., -1., 1.);
925  _3_fold_axis_1_by_3_fold_axis_2.selfNormalize();
926  Matrix1D<double> _3_fold_axis_2_by_4_fold_axis(3);
927  _3_fold_axis_2_by_4_fold_axis = vectorR3(1., 1., 0.);
928  _3_fold_axis_2_by_4_fold_axis.selfNormalize();
929  Matrix1D<double> _4_fold_axis_by_3_fold_axis_1(3);
930  _4_fold_axis_by_3_fold_axis_1 = vectorR3(-1., 1., 0.);
931  _4_fold_axis_by_3_fold_axis_1.selfNormalize();
932  for (size_t i = 0; i < sampling_points_angles.size(); i++)
933  {
934  if (XX(sampling_points_angles[i]) >= 90. &&
935  XX(sampling_points_angles[i]) <= 135. &&
936  YY(sampling_points_angles[i]) <= 90.)
937  if (
938  dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_3_fold_axis_2) >= 0 &&
939  dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_4_fold_axis) >= 0 &&
940  dotProduct(sampling_points_vector[i], _4_fold_axis_by_3_fold_axis_1) >= 0
941  )
942  {
943  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
944  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
945  }
946  }// for i
947  }
948  else if (symmetry == pg_I || symmetry == pg_I2)
949  {//OK
950  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
951  _5_fold_axis_1_by_5_fold_axis_2 = vectorR3(0., 1., 0.);
952  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
953  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
954  _5_fold_axis_2_by_3_fold_axis = vectorR3(-0.4999999839058737,
955  -0.8090170074556163,
956  0.3090169861701543);
957  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
958  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
959  _3_fold_axis_by_5_fold_axis_1 = vectorR3(0.4999999839058737,
960  -0.8090170074556163,
961  0.3090169861701543);
962  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
963  for (size_t i = 0; i < sampling_points_angles.size(); i++)
964  {
965  if (
966  dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
967  dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
968  dotProduct(sampling_points_vector[i], _3_fold_axis_by_5_fold_axis_1) >= 0
969  )
970  {
971  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
972  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
973  }
974  }// for i
975  }
976  else if (symmetry == pg_I1)
977  {//OK
978  Matrix2D<double> A(3, 3);
979  Euler_angles2matrix(0., 90., 0., A);
980  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
981  _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 1., 0.);
982  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
983  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
984  _5_fold_axis_2_by_3_fold_axis = A * vectorR3(-0.4999999839058737,
985  -0.8090170074556163,
986  0.3090169861701543);
987  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
988  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
989  _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.4999999839058737,
990  -0.8090170074556163,
991  0.3090169861701543);
992  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
993 
994  for (size_t i = 0; i < sampling_points_angles.size(); i++)
995  {
996  if (
997  dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
998  dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
999  dotProduct(sampling_points_vector[i], _3_fold_axis_by_5_fold_axis_1) >= 0
1000  )
1001  {
1002  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1003  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1004  }
1005  }// for i
1006  }
1007  else if (symmetry == pg_I3)
1008  {//OK
1009  Matrix2D<double> A(3, 3);
1010  Euler_angles2matrix(0.,31.7174745559,0., A);
1011  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
1012  _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 1., 0.);
1013  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
1014  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
1015  _5_fold_axis_2_by_3_fold_axis = A * vectorR3(-0.4999999839058737,
1016  -0.8090170074556163,
1017  0.3090169861701543);
1018  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
1019  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
1020  _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.4999999839058737,
1021  -0.8090170074556163,
1022  0.3090169861701543);
1023  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
1024 
1025  for (size_t i = 0; i < sampling_points_angles.size(); i++)
1026  {
1027  if (
1028  dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
1029  dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
1030  dotProduct(sampling_points_vector[i], _3_fold_axis_by_5_fold_axis_1) >= 0
1031  )
1032  {
1033  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1034  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1035  }
1036  }// for i
1037  }
1038  else if (symmetry == pg_I4)
1039  {//OK
1040  Matrix2D<double> A(3, 3);
1041  Euler_angles2matrix(0.,-31.7174745559,0., A);
1042  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
1043  _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 0., 1.);
1044  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
1045  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
1046  _5_fold_axis_2_by_3_fold_axis = A * vectorR3(0.187592467856686,
1047  -0.303530987314591,
1048  -0.491123477863004);
1049  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
1050  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
1051  _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.187592467856686,
1052  0.303530987314591,
1053  -0.491123477863004);
1054  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
1055 
1056  for (size_t i = 0; i < sampling_points_angles.size(); i++)
1057  {
1058  if (
1059  dotProduct(sampling_points_vector[i],
1060  _5_fold_axis_2_by_3_fold_axis) <= 0 &&
1061  dotProduct(sampling_points_vector[i],
1062  _3_fold_axis_by_5_fold_axis_1) <= 0 &&
1063  dotProduct(sampling_points_vector[i],
1064  _5_fold_axis_1_by_5_fold_axis_2) <= 0
1065  )
1066  {
1067  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1068  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1069  }
1070  }// for i
1071  }
1072  else if (symmetry == pg_I5)
1073  {//OK
1074  std::cerr << "ERROR: Symmetry pg_I5 not implemented" << std::endl;
1075  exit(0);
1076  }
1077  else if (symmetry == pg_IH || symmetry == pg_I2H)
1078  {//OK
1079  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
1080  _5_fold_axis_1_by_5_fold_axis_2 = vectorR3(0., 1., 0.);
1081  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
1082  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
1083  _5_fold_axis_2_by_3_fold_axis = vectorR3(-0.4999999839058737,
1084  -0.8090170074556163,
1085  0.3090169861701543);
1086  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
1087  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
1088  _3_fold_axis_by_5_fold_axis_1 = vectorR3(0.4999999839058737,
1089  -0.8090170074556163,
1090  0.3090169861701543);
1091  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
1092  Matrix1D<double> _3_fold_axis_by_2_fold_axis(3);
1093  _3_fold_axis_by_2_fold_axis = vectorR3(1.,0.,0.);
1094  _3_fold_axis_by_2_fold_axis.selfNormalize();
1095  for (size_t i = 0; i < sampling_points_angles.size(); i++)
1096  {
1097  if (
1098  dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
1099  dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
1100  dotProduct(sampling_points_vector[i], _3_fold_axis_by_2_fold_axis) >= 0
1101  )
1102  {
1103  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1104  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1105  }
1106  }// for i
1107  }
1108  else if (symmetry == pg_I1H)
1109  {//OK
1110  Matrix2D<double> A(3, 3);
1111  Euler_angles2matrix(0., 90., 0., A);
1112  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
1113  _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 1., 0.);
1114  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
1115  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
1116  _5_fold_axis_2_by_3_fold_axis = A * vectorR3(-0.4999999839058737,
1117  -0.8090170074556163,
1118  0.3090169861701543);
1119  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
1120  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
1121  _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.4999999839058737,
1122  -0.8090170074556163,
1123  0.3090169861701543);
1124  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
1125  Matrix1D<double> _3_fold_axis_by_2_fold_axis(3);
1126  _3_fold_axis_by_2_fold_axis = A * vectorR3(1.,0.,0.);
1127  _3_fold_axis_by_2_fold_axis.selfNormalize();
1128  for (size_t i = 0; i < sampling_points_angles.size(); i++)
1129  {
1130  if (
1131  dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
1132  dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
1133  dotProduct(sampling_points_vector[i], _3_fold_axis_by_2_fold_axis) >= 0
1134  )
1135  {
1136  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1137  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1138  }
1139  }// for i
1140  }
1141  else if (symmetry == pg_I3H)
1142  {//OK
1143  Matrix2D<double> A(3, 3);
1144  Euler_angles2matrix(0.,31.7174745559,0., A);
1145  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
1146  _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 0., 1.);
1147  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
1148  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
1149  _5_fold_axis_2_by_3_fold_axis = A * vectorR3(0.187592467856686,
1150  -0.303530987314591,
1151  -0.491123477863004);
1152  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
1153  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
1154  _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.187592467856686,
1155  0.303530987314591,
1156  -0.491123477863004);
1157  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
1158  Matrix1D<double> _3_fold_axis_by_2_fold_axis(3);
1159  _3_fold_axis_by_2_fold_axis = vectorR3(0.,1.,0.);
1160  _3_fold_axis_by_2_fold_axis.selfNormalize();
1161  for (size_t i = 0; i < sampling_points_angles.size(); i++)
1162  {
1163  if (
1164  dotProduct(sampling_points_vector[i],
1165  _5_fold_axis_2_by_3_fold_axis) >= 0 &&
1166  dotProduct(sampling_points_vector[i],
1167  _3_fold_axis_by_5_fold_axis_1) >= 0 &&
1168  dotProduct(sampling_points_vector[i],
1169  _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
1170  dotProduct(sampling_points_vector[i],
1171  _3_fold_axis_by_2_fold_axis) >= 0
1172  )
1173  {
1174  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1175  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1176  }
1177  }// for i
1178  }
1179  else if (symmetry == pg_I4H)
1180  {//OK
1181  Matrix2D<double> A(3, 3);
1182  Euler_angles2matrix(0.,-31.7174745559,0., A);
1183  Matrix1D<double> _5_fold_axis_1_by_5_fold_axis_2(3);
1184  _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 0., 1.);
1185  _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
1186  Matrix1D<double> _5_fold_axis_2_by_3_fold_axis(3);
1187  _5_fold_axis_2_by_3_fold_axis = A * vectorR3(0.187592467856686,
1188  -0.303530987314591,
1189  -0.491123477863004);
1190  _5_fold_axis_2_by_3_fold_axis.selfNormalize();
1191  Matrix1D<double> _3_fold_axis_by_5_fold_axis_1(3);
1192  _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.187592467856686,
1193  0.303530987314591,
1194  -0.491123477863004);
1195  _3_fold_axis_by_5_fold_axis_1.selfNormalize();
1196  Matrix1D<double> _3_fold_axis_by_2_fold_axis(3);
1197  _3_fold_axis_by_2_fold_axis = vectorR3(0.,1.,0.);
1198  _3_fold_axis_by_2_fold_axis.selfNormalize();
1199  for (size_t i = 0; i < sampling_points_angles.size(); i++)
1200  {
1201  if (
1202  dotProduct(sampling_points_vector[i],
1203  _5_fold_axis_2_by_3_fold_axis) <= 0 &&
1204  dotProduct(sampling_points_vector[i],
1205  _3_fold_axis_by_5_fold_axis_1) <= 0 &&
1206  dotProduct(sampling_points_vector[i],
1207  _5_fold_axis_1_by_5_fold_axis_2) <= 0 &&
1208  dotProduct(sampling_points_vector[i],
1209  _3_fold_axis_by_2_fold_axis) >= 0
1210  )
1211  {
1212  no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
1213  no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
1214  }
1215  }// for i
1216  }
1217  else if (symmetry == pg_I5H)
1218  {//OK
1219  std::cerr << "ERROR: pg_I5H Symmetry not implemented" << std::endl;
1220  exit(0);
1221  }
1222  else
1223  {
1224  std::cerr << "ERROR: Symmetry " << symmetry << "is not known" << std::endl;
1225  exit(0);
1226  }
1227 
1228  CREATE_INDEXES();
1229  numberSamplesAsymmetricUnit=no_redundant_sampling_points_vector.size();
1230  //#define DEBUG5
1231 #ifdef DEBUG5
1232 
1233  std::ofstream filestr5;
1234  filestr5.open ("removeRedundantPoints.bild");
1235  filestr5 << ".color 0 0 1" << std::endl;
1236  for (int i = 0;
1237  i < no_redundant_sampling_points_vector.size();
1238  i++)
1239  {
1240  filestr5 << ".color 1 0 0" << std::endl
1241  << ".sphere "
1242  << no_redundant_sampling_points_vector[i](0) << " "
1243  << no_redundant_sampling_points_vector[i](1) << " "
1244  << no_redundant_sampling_points_vector[i](2) << " "
1245  << " .03" << std::endl;
1246  }
1247  filestr5.close();
1248 #endif
1249 #undef DEBUG5
1250 
1251 }
1252 
1253 void Sampling::removeRedundantPointsExhaustive(const int symmetry,
1254  int sym_order,
1255  bool only_half_sphere,
1256  double max_ang)
1257 {
1258  // Maximum distance
1259  double cos_max_ang = cos(DEG2RAD(max_ang));
1260  double my_dotProduct;
1261  //int j_end=0;
1262  Matrix1D<double> direction(3), direction1(3);
1263 
1264  // First call to conventional removeRedundantPoints
1265  removeRedundantPoints(symmetry, sym_order);
1266  // for (int isym = 0; isym < no_redundant_sampling_points_vector.size(); isym++)
1267  // std::cout << "sampling:" << no_redundant_sampling_points_vector[isym];
1268  std::vector <Matrix1D<double> > old_vector = no_redundant_sampling_points_vector;
1269  std::vector <Matrix1D<double> > old_angles = no_redundant_sampling_points_angles;
1270 
1271  CLEAR_VECTORS();
1272 
1273  // Precalculate symmetry matrices
1274  fillLRRepository();
1275 
1276  // Then check all points versus each other
1277  for (size_t i = 0; i < old_angles.size(); i++)
1278  {
1279  //direction1=(old_vector[i]).transpose();
1280  direction1=old_vector[i];
1281  bool uniq = true;
1282  for (size_t j = 0; j < R_repository.size(); j++)
1283  {
1284  for (size_t k = 0; k < no_redundant_sampling_points_vector.size(); k++)
1285  {
1286  direction = L_repository[j] *
1287  (no_redundant_sampling_points_vector[k].transpose() *
1288  R_repository[j]).transpose();
1289  //Calculate distance
1290  my_dotProduct = dotProduct(direction,direction1);
1291  if (only_half_sphere)
1292  my_dotProduct = ABS(my_dotProduct);
1293 
1294  if (my_dotProduct > cos_max_ang)
1295  {
1296  uniq = false;
1297  break;
1298  }
1299  }// for k
1300  if (!uniq)
1301  break;
1302  } // for j
1303  if (uniq)
1304  {
1305  no_redundant_sampling_points_vector.push_back(old_vector[i]);
1306  no_redundant_sampling_points_angles.push_back(old_angles[i]);
1307  }
1308  } // for i
1309 
1310  CREATE_INDEXES();
1311 
1312 }
1313 
1314 
1315 //THIs FUNCTION IS NOW OBSOLETE
1316 //SINCE readSymmetryFile does not longer need a file
1317 //use symmetry functions instead
1318 /* Create symmetry file----------------------------------------------------- */
1319 void Sampling::createSymFile(const FileName &simFp,int symmetry, int sym_order)
1320 {
1321  symmetry_file = simFp + ".sym";
1322  std::ofstream SymFile;
1323  SymFile.open(symmetry_file.c_str(), std::ios::out);
1324  if (symmetry == pg_CN)
1325  {
1326  SymFile << "rot_axis " << sym_order << " 0 0 1";
1327  }
1328  else if (symmetry == pg_CI)
1329  {
1330  SymFile << "inversion ";
1331  }
1332  else if (symmetry == pg_CS)
1333  {
1334  SymFile << "mirror_plane 0 0 1";
1335  }
1336  else if (symmetry == pg_CNV)
1337  {
1338  SymFile << "rot_axis " << sym_order << " 0 0 1";
1339  SymFile << std::endl;
1340  SymFile << "mirror_plane 0 1 0";
1341  }
1342  else if (symmetry == pg_CNH)
1343  {
1344  SymFile << "rot_axis " << sym_order << " 0 0 1";
1345  SymFile << std::endl;
1346  SymFile << "mirror_plane 0 0 -1";
1347  }
1348  else if (symmetry == pg_SN)
1349  {
1350  int order = sym_order / 2;
1351  SymFile << "rot_axis " << order << " 0 0 1";
1352  SymFile << std::endl;
1353  SymFile << "inversion ";
1354  }
1355  else if (symmetry == pg_DN)
1356  {
1357  SymFile << "rot_axis " << sym_order << " 0 0 1";
1358  SymFile << std::endl;
1359  SymFile << "rot_axis " << "2" << " 0 1 0";
1360  }
1361  else if (symmetry == pg_DNV)
1362  {
1363  SymFile << "rot_axis " << sym_order << " 0 0 1";
1364  SymFile << std::endl;
1365  SymFile << "rot_axis " << "2" << " 0 1 0";
1366  SymFile << std::endl;
1367  SymFile << "mirror_plane 0 1 0";
1368  }
1369  else if (symmetry == pg_DNH)
1370  {
1371  SymFile << "rot_axis " << sym_order << " 0 0 1";
1372  SymFile << std::endl;
1373  SymFile << "rot_axis " << "2" << " 1 0 0";
1374  SymFile << std::endl;
1375  SymFile << "mirror_plane 0 0 1";
1376  }
1377  else if (symmetry == pg_T)
1378  {
1379  SymFile << "rot_axis " << "3" << " 0. 0. 1.";
1380  SymFile << std::endl;
1381  SymFile << "rot_axis " << "2" << " 0. 0.816496 0.577350";
1382  }
1383  else if (symmetry == pg_TD)
1384  {
1385  SymFile << "rot_axis " << "3" << " 0. 0. 1.";
1386  SymFile << std::endl;
1387  SymFile << "rot_axis " << "2" << " 0. 0.816496 0.577350";
1388  SymFile << std::endl;
1389  SymFile << "mirror_plane 1.4142136 2.4494897 0.0000000";
1390  }
1391  else if (symmetry == pg_TH)
1392  {
1393  SymFile << "rot_axis " << "3" << " 0. 0. 1.";
1394  SymFile << std::endl;
1395  SymFile << "rot_axis " << "2" << " 0. -0.816496 -0.577350";
1396  SymFile << std::endl;
1397  SymFile << "inversion";
1398  }
1399  else if (symmetry == pg_O)
1400  {
1401  SymFile << "rot_axis " << "3" << " .5773502 .5773502 .5773502";
1402  SymFile << std::endl;
1403  SymFile << "rot_axis " << "4" << " 0 0 1";
1404  }
1405  else if (symmetry == pg_OH)
1406  {
1407  SymFile << "rot_axis " << "3" << " .5773502 .5773502 .5773502";
1408  SymFile << std::endl;
1409  SymFile << "rot_axis " << "4" << " 0 0 1";
1410  SymFile << std::endl;
1411  SymFile << "mirror_plane 0 1 1";
1412  }
1413  else if (symmetry == pg_I)
1414  {
1415  SymFile << "rot_axis 2 0 0 1";
1416  SymFile << std::endl;
1417  SymFile << "rot_axis 5 -1.618033989 -1 0";
1418  SymFile << std::endl;
1419  SymFile << "rot_axis 3 -0.53934467 -1.4120227 0";
1420  }
1421  else if (symmetry == pg_IH)
1422  {
1423  SymFile << "rot_axis 2 0 0 1";
1424  SymFile << std::endl;
1425  SymFile << "rot_axis 5 -1.618033989 -1 0";
1426  SymFile << std::endl;
1427  SymFile << "rot_axis 3 -0.53934467 -1.4120227 0";
1428  SymFile << std::endl;
1429  SymFile << "mirror_plane 1 0 0";
1430  }
1431  else
1432  {
1433  std::cerr << "ERROR:: Symmetry " << symmetry << "is not known" << std::endl;
1434  exit(0);
1435  }
1436  SymFile.close();
1437  SL.readSymmetryFile(symmetry_file);
1438 
1439 }
1440 void Sampling::createAsymUnitFile(const FileName &docfilename)
1441 {
1442  MetaDataVec DF;
1443  FileName tmp_filename;
1444  //#define CHIMERA
1445 #ifdef CHIMERA
1446 
1447  std::ofstream filestr;
1448  filestr.open ("create_asym_unit_file.bild");
1449  filestr << ".color white"
1450  << std::endl
1451  << ".sphere 0 0 0 .95"
1452  << std::endl
1453  ;
1454  filestr << ".color green"
1455  << std::endl
1456  ;
1457 #endif
1458 
1459  MDRowVec row;
1460  for (size_t i = 0; i < no_redundant_sampling_points_vector.size(); i++)
1461  {
1462 #ifdef CHIMERA
1463  filestr << ".sphere "
1464  << no_redundant_sampling_points_vector[i]
1465  << " 0.018"
1466  << std::endl
1467  ;
1468 #endif
1469 
1470  row.setValue(MDL_REF, (int)i);
1471  row.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i]);
1472  row.setValue(MDL_ANGLE_ROT,XX(no_redundant_sampling_points_angles[i]));
1473  row.setValue(MDL_ANGLE_TILT,YY(no_redundant_sampling_points_angles[i]));
1474  row.setValue(MDL_ANGLE_PSI,ZZ(no_redundant_sampling_points_angles[i]));
1475  row.setValue(MDL_X,XX(no_redundant_sampling_points_vector[i]));
1476  row.setValue(MDL_Y,YY(no_redundant_sampling_points_vector[i]));
1477  row.setValue(MDL_Z,ZZ(no_redundant_sampling_points_vector[i]));
1478  DF.addRow(row);
1479  }
1480 #ifdef CHIMERA
1481  filestr.close();
1482 #endif
1483  #undef CHIMERA
1484 
1485  tmp_filename = docfilename + "_angles.doc";
1486  DF.setComment("REF refers to the projection directions BEFORE delete those not in a neighborhood, ");
1487  DF.write(tmp_filename);
1488 }
1489 
1490 #define FN_SAMPLING_NEI(base) formatString("neighbors@%s_sampling.xmd", base.c_str())
1491 #define FN_SAMPLING_PROJ(base) formatString("projectionDirections@%s_sampling.xmd", base.c_str())
1492 #define FN_SAMPLING_EXTRA(base) formatString("extra@%s_sampling.xmd", base.c_str())
1493 #define FN_SAMPLING_SPHERE(base) formatString("projectionDirectionsSphere@%s_sampling.xmd", base.c_str())
1494 
1495 void Sampling::saveSamplingFile(const FileName &fn_base, bool write_vectors, bool write_sampling_sphere)
1496 {
1497  {
1498  MetaDataVec md;
1499  MDRowVec row;
1500  row.setValue(MDL_SAMPLINGRATE, sampling_rate_rad);
1501  row.setValue(MDL_NEIGHBORHOOD_RADIUS, cos_neighborhood_radius);
1502  row.setValue(MDL_POINTSASYMETRICUNIT,numberSamplesAsymmetricUnit);
1503  md.setComment("data_extra -> sampling description;"\
1504  " data_neighbors --> List with order of each"\
1505  "experimental images and its neighbors"\
1506  );
1507  md.setColumnFormat(false);
1508  md.addRow(row);
1509  md.write(FN_SAMPLING_EXTRA(fn_base), MD_OVERWRITE);
1510  }
1511 
1512  {
1513  MetaDataVec md;
1514  md.setColumnFormat(true);
1515  size_t size = my_neighbors.size();
1516  //Write first block with experimental images order and its neighbors
1517  bool writeFileName = !exp_data_fileNames.empty();
1518  for(size_t i = 0; i < size; ++i)
1519  {
1520  MDRowVec row;
1521  row.setValue(MDL_ORDER,i+FIRST_IMAGE);
1522  if (writeFileName)
1523  row.setValue(MDL_IMAGE,exp_data_fileNames[i]);
1524  row.setValue(MDL_NEIGHBORS, my_neighbors[i]);
1525  md.addRow(row);
1526  }
1527 
1528  md.write(FN_SAMPLING_NEI(fn_base), MD_APPEND);
1529  }
1530 
1531  {
1532  //Write projection directions
1533  MetaDataVec md;
1534  size_t size = no_redundant_sampling_points_index.size();
1535 
1536  for (size_t i = 0; i < size; ++i)
1537  {
1538  MDRowVec row;
1539  Matrix1D<double> &angles = no_redundant_sampling_points_angles[i];
1540  row.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i]);
1541  row.setValue(MDL_ANGLE_ROT, XX(angles));
1542  row.setValue(MDL_ANGLE_TILT, YY(angles));
1543  row.setValue(MDL_ANGLE_PSI, ZZ(angles));
1544 
1545  if (write_vectors)
1546  {
1547  Matrix1D<double> &vectors = no_redundant_sampling_points_vector[i];
1548  row.setValue(MDL_X, XX(vectors));
1549  row.setValue(MDL_Y, YY(vectors));
1550  row.setValue(MDL_Z, ZZ(vectors));
1551  }
1552  md.addRow(row);
1553  }
1554  md.write(FN_SAMPLING_PROJ(fn_base), MD_APPEND);
1555  }
1556 
1557  if (write_sampling_sphere)
1558  {
1559  MetaDataVec md;
1560  size_t size = sampling_points_angles.size();
1561 
1562  for (size_t i = 0; i < size; ++i)
1563  {
1564  MDRowVec row;
1565  row.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i]);
1566  Matrix1D<double> &angles = sampling_points_angles[i];
1567  row.setValue(MDL_ANGLE_ROT, XX(angles));
1568  row.setValue(MDL_ANGLE_TILT, YY(angles));
1569  row.setValue(MDL_ANGLE_PSI, ZZ(angles));
1570 
1571 
1572  if (write_vectors)
1573  {
1574  Matrix1D<double> &vectors = sampling_points_vector[i];
1575  row.setValue(MDL_X, XX(vectors));
1576  row.setValue(MDL_Y, YY(vectors));
1577  row.setValue(MDL_Z, ZZ(vectors));
1578  }
1579  md.addRow(row);
1580  }
1581  md.setComment("NEIGHBOR index points to the slice in a stack file in which the projection projected in the direction defined by ref is stored");
1582  md.write(FN_SAMPLING_SPHERE(fn_base), MD_APPEND);
1583 
1584  }
1585 
1586 }
1587 
1588 //angle conventions changed by test not
1589 //label names changed but no test
1590 //readSamplingFile cannot work
1591 
1592 void Sampling::readSamplingFile(const FileName &fn_base,
1593  bool read_vectors,
1594  bool read_sampling_sphere)
1595 {
1596  //Read extra info
1597  MetaDataVec md(FN_SAMPLING_EXTRA(fn_base));
1598  size_t id = md.firstRowId();
1599  md.getValue(MDL_SAMPLINGRATE, sampling_rate_rad, id);
1600  md.getValue(MDL_NEIGHBORHOOD_RADIUS, cos_neighborhood_radius, id);
1601  md.getValue(MDL_POINTSASYMETRICUNIT,numberSamplesAsymmetricUnit,id);
1602 
1603  //Read neighbors
1604  md.read(FN_SAMPLING_NEI(fn_base));
1605  my_neighbors.resize(md.size());
1606  bool readFileName = md.containsLabel(MDL_IMAGE);
1607  if (readFileName)
1608  {
1609  exp_data_fileNames.clear();
1610  exp_data_fileNames.resize(md.size());
1611  }
1612 
1613  {
1614  size_t i = 0, ii = 0;
1615  for (size_t objId : md.ids())
1616  {
1617  if (readFileName)
1618  md.getValue(MDL_IMAGE,exp_data_fileNames[ii++], objId);
1619  md.getValue(MDL_NEIGHBORS, my_neighbors[i], objId);
1620  i++;
1621  }
1622  }
1623 
1624  //Read projection directions
1625  md.read(FN_SAMPLING_PROJ(fn_base));
1626  size_t size = md.size();
1627  no_redundant_sampling_points_index.resize(size);
1628  no_redundant_sampling_points_angles.resize(size);
1629  if (read_vectors)
1630  no_redundant_sampling_points_vector.resize(size);
1631 
1632  {
1633  size_t i = 0;
1634  for (size_t objId : md.ids())
1635  {
1637  //size_t objId;
1639  //size_t objIndex;
1640 
1641  Matrix1D<double> &angles = no_redundant_sampling_points_angles[i];
1642  angles.resizeNoCopy(3);
1643  md.getValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i], objId);
1644  md.getValue(MDL_ANGLE_ROT, XX(angles), objId);
1645  md.getValue(MDL_ANGLE_TILT, YY(angles), objId);
1646  md.getValue(MDL_ANGLE_PSI, ZZ(angles), objId);
1647 
1648  if (read_vectors)
1649  {
1650  Matrix1D<double> &vectors = no_redundant_sampling_points_vector[i];
1651  vectors.resizeNoCopy(3);
1652  md.getValue(MDL_X, XX(vectors), objId);
1653  md.getValue(MDL_Y, YY(vectors), objId);
1654  md.getValue(MDL_Z, ZZ(vectors), objId);
1655  }
1656 
1657  i++;
1658  }
1659  }
1660 
1661 //#define DEBUG5
1662 #ifdef DEBUG5
1663 //This bild file does not make sense since x,y,z are angles
1664  std::ofstream filestr5;
1665  filestr5.open ("/tmp/loadsamplingfile_removeRedundantPoints.bild");
1666  filestr5 << ".color 0 0 1" << std::endl;
1667  for (int i = 0;
1668  i < no_redundant_sampling_points_index.size();
1669  i++)
1670  {
1671  filestr5 << ".color 1 0 0" << std::endl
1672  << ".sphere "
1673  << no_redundant_sampling_points_vector[i](0) << " "
1674  << no_redundant_sampling_points_vector[i](1) << " "
1675  << no_redundant_sampling_points_vector[i](2) << " "
1676  << " .03" << std::endl;
1677  }
1678  filestr5.close();
1679 #endif
1680 #undef DEBUG5
1681 
1682  if (read_sampling_sphere)
1683  {
1684  //Read projection directions
1685  md.read(FN_SAMPLING_SPHERE(fn_base));
1686  size_t size = md.size();
1687  sampling_points_angles.resize(size);
1688  if (read_vectors)
1689  sampling_points_vector.resize(size);
1690 
1691  size_t i = 0;
1692  for (size_t objId : md.ids())
1693  {
1694  Matrix1D<double> &angles = sampling_points_angles[i];
1695  angles.resizeNoCopy(3);
1696  md.getValue(MDL_ANGLE_ROT, XX(angles), objId);
1697  md.getValue(MDL_ANGLE_TILT, YY(angles), objId);
1698  md.getValue(MDL_ANGLE_PSI, ZZ(angles), objId);
1699 
1700  if (read_vectors)
1701  {
1702  Matrix1D<double> &vectors = sampling_points_vector[i];
1703  vectors.resizeNoCopy(3);
1704  md.getValue(MDL_X, XX(vectors), objId);
1705  md.getValue(MDL_Y, YY(vectors), objId);
1706  md.getValue(MDL_Z, ZZ(vectors), objId);
1707  }
1708 
1709  i++;
1710  }
1711 
1712  }
1713 }
1714 
1715 void Sampling::computeNeighbors(bool only_winner)
1716 {
1717 
1718  double my_dotProduct;
1719  double winner_dotProduct;
1720  Matrix2D<double> L(4, 4), R(4, 4);
1721  std::vector<size_t> aux_neighbors;
1722  bool new_reference=true;
1723  my_neighbors.clear();
1724 #ifdef MYPSI
1725  std::vector<double> aux_neighbors_psi;
1726  my_neighbors_psi.clear();
1727 #endif
1728 
1729  // calculate some sizes only once
1730  size_t exp_data_projection_direction_by_L_R_size = exp_data_projection_direction_by_L_R.size();
1731  size_t no_redundant_sampling_points_vector_size = no_redundant_sampling_points_vector.size();
1732 
1733  if (verbose)
1734  {
1735  std::cout << "Find valid sampling points based on the neighborhood" <<std::endl;
1736  init_progress_bar(exp_data_projection_direction_by_L_R_size);
1737  }
1738  size_t ratio = exp_data_projection_direction_by_L_R_size / 60;
1739  ratio = XMIPP_MAX(ratio, 1);
1740 
1741  for(size_t j = 0; j < exp_data_projection_direction_by_L_R_size;)
1742  {
1743  if ((j%ratio) == 0 && verbose)
1744  progress_bar(j);
1745 
1746 #ifdef MYPSI
1747 
1748  aux_neighbors_psi.clear();
1749 #endif
1750 
1751  aux_neighbors.clear();
1752  if (cos_neighborhood_radius <= -1.0)
1753  {
1754  aux_neighbors=no_redundant_sampling_points_index;
1755  j+=R_repository.size();
1756  }
1757  else
1758  {
1759  size_t * aux_neighborsArray = nullptr;
1760  for (size_t k = 0; k < R_repository.size(); k++,j++)
1761  {
1762  winner_dotProduct = -1.;
1763  for (size_t i = 0; i < no_redundant_sampling_points_vector_size; ++i)
1764  {
1765  my_dotProduct = dotProduct(no_redundant_sampling_points_vector[i],
1766  exp_data_projection_direction_by_L_R[j]);
1767 
1768  if (my_dotProduct > cos_neighborhood_radius)
1769  {
1770  if(aux_neighbors.size()==0)
1771  {
1772  aux_neighbors.push_back(no_redundant_sampling_points_index[i]);
1773  winner_dotProduct=my_dotProduct;
1774 
1775  #ifdef MYPSI
1776 
1777  aux_neighbors_psi.push_back(exp_data_projection_direction_by_L_R_psi[j]);
1778  #endif
1779 
1780  }
1781  else
1782  {
1783  new_reference = true;
1784  if(only_winner)
1785  {
1786  if(winner_dotProduct<my_dotProduct)
1787  {
1788  if(winner_dotProduct!=-1)
1789  aux_neighbors.pop_back();
1790  #ifdef MYPSI
1791 
1792  if(winner_dotProduct!=-1)
1793  aux_neighbors_psi.pop_back();
1794  #endif
1795 
1796  winner_dotProduct=my_dotProduct;
1797  }
1798  else
1799  {
1800  new_reference=false;
1801  }
1802  }
1803  else
1804  {
1805  //precalculate size saves time here but
1806  //not in the whole loop
1807  aux_neighborsArray = &aux_neighbors[0];
1808  size_t _size = aux_neighbors.size();
1809  //std::cerr << "DEBUG_JM: no_redundant_sampling_points_index[i]: " << no_redundant_sampling_points_index[i] << std::endl;
1810  // for (size_t kkk = 0; kkk < aux_neighbors.size(); ++kkk)
1811  // std::cerr << aux_neighbors[kkk] << " ";
1812  for( size_t l=0;l< _size;l++)
1813  {
1814  //if (aux_neighbors[l]==i)
1815  if (aux_neighborsArray[l]==no_redundant_sampling_points_index[i])
1816  {
1817  new_reference=false;
1818  break;
1819  }
1820  }
1821  //std::cerr << "DEBUG_JM: new_reference: " << new_reference << std::endl;
1822  }
1823  if (new_reference)
1824  {
1825  //std::cerr << formatString("DEBUG_JM: j %lu k %lu i %lu ", j, k, i) << std::endl;
1826  aux_neighbors.push_back(no_redundant_sampling_points_index[i]);
1827 
1828  //for (size_t kkk = 0; kkk < aux_neighbors.size(); ++kkk)
1829  // std::cerr << aux_neighbors[kkk] << " ";
1830  // std::cerr << std::endl;
1831 
1832  #ifdef MYPSI
1833 
1834  aux_neighbors_psi.push_back(exp_data_projection_direction_by_L_R_psi[j]);
1835  #endif
1836 
1837  }
1838  }
1839  //same sampling point should appear only once
1840  //note that psi recorded here may be different from psi
1841  //recorded in _closest_sampling_points because
1842  //may refer to a different sampling point
1843  //in fact every point is degenerated
1844  }
1845  }//for i;
1846  }//for k
1847  }
1848  my_neighbors.push_back(aux_neighbors);
1849 #ifdef MYPSI
1850 
1851  my_neighbors_psi.push_back(aux_neighbors_psi);
1852 #endif
1853 
1854  }//for j
1855  if (verbose)
1856  progress_bar(exp_data_projection_direction_by_L_R_size);
1857 
1858  //#define DEBUG
1859 #ifdef DEBUG
1860 
1861  for (int i=0;i< my_neighbors.size();i++)
1862  for (int j=0;j< my_neighbors[i].size();j++)
1863  std::cerr << "image:" << i << " "<< my_neighbors[i][j]<<std::endl;
1864  exit(1);
1865 #endif
1866 #undef DEBUG
1867 
1868 //#define CHIMERA
1869 #ifdef CHIMERA
1870 
1871  std::ofstream filestr;
1872  filestr.open ("compute_neighbors.bild");
1873  filestr << ".color white"
1874  << std::endl << ".sphere 0 0 0 .95" << std::endl
1875  ;
1876  int exp_image=0;
1877  filestr << ".color yellow" << std::endl
1878  << ".sphere "
1879  << exp_data_projection_direction_by_L_R[exp_image*R_repository.size()](0) << " "
1880  << exp_data_projection_direction_by_L_R[exp_image*R_repository.size()](1) << " "
1881  << exp_data_projection_direction_by_L_R[exp_image*R_repository.size()](2) << " "
1882  << " .021"
1883  << std::endl;
1884  for(int i=(exp_image*R_repository.size());
1885  i< (exp_image+1)*R_repository.size();
1886  i++)
1887  {
1888  filestr << ".color red" << std::endl
1889  << ".sphere "
1890 // << exp_data_projection_direction_by_L_R[i]
1891  << exp_data_projection_direction_by_L_R[i](0) << " "
1892  << exp_data_projection_direction_by_L_R[i](1) << " "
1893  << exp_data_projection_direction_by_L_R[i](2) << " "
1894  << " .017" << std::endl;
1895  }
1896  double blue;
1897  for(int i=0;
1898  i< my_neighbors[exp_image].size();
1899  i++)
1900  {
1901 #ifdef MYPSI
1902  blue = (my_neighbors_psi[exp_image][i]+180.)/360.;
1903 #else
1904 
1905  blue = 1.;
1906 #endif
1907 
1908  filestr << ".color 0 0 " << blue << std::endl
1909  << ".sphere "
1910 // << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]]
1911  << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]](0) << " "
1912  << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]](1) << " "
1913  << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]](2) << " "
1914 
1915  << " .019" << std::endl;
1916  //std::cerr << my_neighbors_psi[exp_image][i] << std::endl;
1917  }
1918  filestr.close();
1919 
1920 #endif
1921  #undef CHIMERA
1922 
1923 }
1926 #define REMOVE_LAST(vector) vector[i] = vector[my_end]; vector.pop_back();
1927 
1928 void Sampling::removePointsFarAwayFromExperimentalData()
1929 {
1930  double my_dotProduct;
1931  Matrix1D<double> row(3),direction(3);
1932  Matrix2D<double> L(4, 4), R(4, 4);
1933 
1934  size_t my_end = no_redundant_sampling_points_vector.size() - 1;
1935 
1936  for (size_t i = 0; i <= my_end; i++)
1937  {
1938  bool my_delete=true;
1939  for (size_t j=0; my_delete && j< exp_data_projection_direction_by_L_R.size();j++)
1940  {
1941  my_dotProduct = dotProduct(no_redundant_sampling_points_vector[i],
1942  exp_data_projection_direction_by_L_R[j]);
1943  if (my_dotProduct > cos_neighborhood_radius)
1944  my_delete=false;
1945  }//for j
1946  if(my_delete)
1947  {
1948  REMOVE_LAST(no_redundant_sampling_points_vector);
1949  REMOVE_LAST(no_redundant_sampling_points_angles);
1950  REMOVE_LAST(no_redundant_sampling_points_index);
1951 
1952  --my_end;
1953  --i;//since a point has been swaped we should repeat the same index
1954  }// if(my_delete)
1955  }//for i end
1956  //#define CHIMERA
1957 #ifdef CHIMERA
1958  std::ofstream filestr;
1959  filestr.open ("remove_points_far_away_from_experimental_data.bild");
1960  filestr << ".color white"
1961  << std::endl
1962  << ".sphere 0 0 0 .95"
1963  << std::endl
1964  ;
1965  filestr << ".color green"
1966  << std::endl
1967  ;
1968  //green neighbours
1969  for (int i = 0;
1970  i < no_redundant_sampling_points_vector.size();
1971  i++)
1972  {
1973  filestr << ".color green" << std::endl
1974  << ".sphere " << no_redundant_sampling_points_vector[i].transpose() <<
1975  " .018" << std::endl;
1976  }
1977  filestr.close();
1978 
1979 #endif
1980  #undef CHIMERA
1981 }
1982 void Sampling::findClosestSamplingPoint(const FileName &FnexperimentalImages,
1983  const FileName &output_file_root)
1984 {
1985  //read input files
1986  MetaDataVec DFi;
1987  DFi.read(FnexperimentalImages);//experimental points
1988  findClosestSamplingPoint(DFi, output_file_root);
1989 
1990 }
1991 void Sampling::findClosestSamplingPoint(const MetaData &DFi,
1992  const FileName &output_file_root)
1993 {
1994  double my_dotProduct,my_dotProduct_aux;
1995  Matrix1D<double> row(3),direction(3);
1996  Matrix1D<double> docline;
1997  docline.initZeros(7);//three original angles, one winnir, new angles
1998  Matrix2D<double> L(4, 4), R(4, 4);
1999  int winner_sampling=-1;
2000 #if defined(CHIMERA) || defined(MYPSI)
2001 
2002  int winner_exp_L_R=-1;
2003 #endif
2004 
2005  MetaDataVec DFo;
2006  size_t id;
2007 
2008  DFo.setComment("Original rot, tilt, psi, Xoff, Yoff are stored as comments");
2009 
2010  //#define DEBUG3
2011 #ifdef DEBUG3
2012 
2013  std::ofstream filestr;
2014  filestr.open ("find_closest_sampling_point.bild");
2015  int exp_image=1;
2016 #endif
2017 
2018  auto idIter(DFi.ids().begin());
2019  for (size_t i = 0; i < exp_data_projection_direction_by_L_R.size(); )
2020  {
2021  my_dotProduct=-2;
2022  for (size_t k = 0; k < R_repository.size(); k++,i++)
2023  {
2024 #ifdef DEBUG3
2025  //experimental points plus symmetry
2026  if( i>(exp_image*R_repository.size()-1) && i< ((exp_image+1)*R_repository.size()))
2027  {
2028  filestr << ".color red" << std::endl
2029  << ".sphere " << exp_data_projection_direction_by_L_R[i]
2030  << " .019" << std::endl;
2031  }
2032 #endif
2033  for(size_t j=0;j< no_redundant_sampling_points_vector.size();j++)
2034  {
2035  my_dotProduct_aux =
2036  dotProduct(exp_data_projection_direction_by_L_R[i],
2037  no_redundant_sampling_points_vector[j]);
2038 
2039  if ( my_dotProduct_aux > my_dotProduct)
2040  {
2041  my_dotProduct = my_dotProduct_aux;
2042  winner_sampling = j;
2043 #if defined(CHIMERA) || defined(MYPSI)
2044 
2045  winner_exp_L_R = i;
2046 #endif
2047 
2048  }
2049  }//for j
2050  }//for k
2051 #ifdef DEBUG3
2052  if( i== ((exp_image+1)*R_repository.size()) )
2053  {
2054  filestr << ".color yellow" << std::endl
2055  << ".sphere " << no_redundant_sampling_points_vector[winner_sampling]
2056  << " .020" << std::endl;
2057  }
2058 #endif
2059  //add winner to the DOC fILE
2060  std::string fnImg, comment;
2061  double aux;
2062  size_t objId = *idIter;
2063  DFi.getValue(MDL_IMAGE, fnImg, objId);
2064  DFi.getValue(MDL_ANGLE_ROT,aux, objId);
2065  comment+=floatToString(aux)+" ";
2066  DFi.getValue(MDL_ANGLE_TILT,aux, objId);
2067  comment+=floatToString(aux)+" ";
2068  DFi.getValue(MDL_ANGLE_PSI,aux, objId);
2069  comment+=floatToString(aux)+" ";
2070  DFi.getValue(MDL_SHIFT_X,aux, objId);
2071  comment+=floatToString(aux)+" ";
2072  DFi.getValue(MDL_SHIFT_Y,aux, objId);
2073  comment+=floatToString(aux);
2074  id = DFo.addObject();
2075  DFo.setValue(MDL_STAR_COMMENT,comment, id);
2076  DFo.setValue(MDL_IMAGE,fnImg, id);
2077  DFo.setValue(MDL_REF, winner_sampling, id);
2078 #ifdef MYPSI
2079 
2080  DFo.set(6, exp_data_projection_direction_by_L_R_psi[winner_exp_L_R]);
2081 #endif
2082 
2083  DFo.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[winner_sampling], id);
2084  DFo.setValue(MDL_ANGLE_ROT,XX(no_redundant_sampling_points_angles[winner_sampling]), id);
2085  DFo.setValue(MDL_ANGLE_TILT,YY(no_redundant_sampling_points_angles[winner_sampling]), id);
2086  DFo.setValue(MDL_ANGLE_PSI,ZZ(no_redundant_sampling_points_angles[winner_sampling]), id);
2087 
2088  ++idIter;
2089  }//for i
2090 
2091  if (output_file_root.size() > 0)
2092  DFo.write(output_file_root+ "_closest_sampling_points.doc");
2093 #ifdef DEBUG3
2094 
2095  filestr.close();
2096 #endif
2097 #undef DEBUG3
2098 }
2099 
2100 void Sampling::findClosestExperimentalPoint()
2101 {
2102  double my_dotProduct,my_dotProduct_aux;
2103  Matrix1D<double> row(3),direction(3);
2104  int winner_sampling=-1;
2105  int winner_exp=-1;
2106  //#define CHIMERA
2107 #ifdef CHIMERA
2108 
2109  std::vector<std::vector<int> > aux_vec;
2110  aux_vec.resize(no_redundant_sampling_points_vector.size());
2111 #endif
2112 
2113  std::vector<std::vector<size_t> > aux_my_exp_img_per_sampling_point;
2114 
2115  //resize vector
2116  aux_my_exp_img_per_sampling_point.resize(
2117  no_redundant_sampling_points_vector.size());
2118 
2119  for(size_t i=0,l=0;i< exp_data_projection_direction_by_L_R.size();l++)
2120  {
2121  my_dotProduct=-2;
2122  for (size_t k = 0; k < R_repository.size(); k++,i++)
2123  {
2124  for(size_t j=0;j< no_redundant_sampling_points_vector.size();j++)
2125  {
2126  my_dotProduct_aux =
2127  dotProduct(exp_data_projection_direction_by_L_R[i],
2128  no_redundant_sampling_points_vector[j]);
2129 
2130  if ( my_dotProduct_aux > my_dotProduct)
2131  {
2132  my_dotProduct = my_dotProduct_aux;
2133  winner_sampling = j;
2134 #ifdef CHIMERA
2135 
2136  winner_exp_L_R = i;
2137 #endif
2138 
2139  winner_exp = l;
2140  }
2141  }//for j
2142  }//for k
2143  aux_my_exp_img_per_sampling_point[winner_sampling].push_back(winner_exp);
2144 #ifdef CHIMERA
2145 
2146  aux_vec[winner_sampling].push_back(winner_exp_L_R);
2147 #endif
2148 
2149  }//for i aux_my_exp_img_per_sampling_point
2150  for(size_t i=0;i< aux_my_exp_img_per_sampling_point.size();i++)
2151  if(aux_my_exp_img_per_sampling_point[i].size()!=0)
2152  my_exp_img_per_sampling_point.push_back(aux_my_exp_img_per_sampling_point[i]);
2153 #ifdef CHIMERA
2154 
2155  std::ofstream filestr;
2156  filestr.open ("find_closest_experimental_point.bild");
2157  filestr << ".color white"
2158  << std::endl
2159  << ".sphere 0 0 0 .95"
2160  << std::endl
2161  ;
2162  filestr << ".color red"
2163  << std::endl
2164  ;
2165  for (int i = 0;
2166  i < no_redundant_sampling_points_vector.size();
2167  i++)
2168  {
2169  filestr << ".sphere " << no_redundant_sampling_points_vector[i].transpose() <<
2170  " .018" << std::endl;
2171  }
2172  int my_sample_point=5;
2173  filestr << ".color green"
2174  << std::endl
2175  ;
2176  int ii;
2177  for (int i = 0;
2178  i < my_exp_img_per_sampling_point[my_sample_point].size();
2179  i++)
2180  {
2181  ii=aux_vec[my_sample_point][i];
2182  filestr << ".sphere "
2183  << exp_data_projection_direction_by_L_R[ii].transpose()
2184  << " .017" << std::endl;
2185  }
2186  filestr.close();
2187 
2188 #endif
2189  #undef CHIMERA
2190  //#define DEBUG4
2191 #ifdef DEBUG4
2192 
2193  std::ofstream filestr;
2194  filestr.open ("find_closest_experimental_point.txt");
2195 
2196  for (int i = 0;
2197  i < my_exp_img_per_sampling_point.size();
2198  i++)
2199  { //for each sampling point write its experimental images
2200  filestr << i << std::endl;
2201  for (int j = 0;
2202  j < my_exp_img_per_sampling_point[i].size();
2203  j++)
2204  {
2205  filestr << my_exp_img_per_sampling_point[i][j]
2206  << " " ;
2207  }
2208  filestr << std::endl;
2209  }
2210  filestr.close();
2211 
2212 #endif
2213  #undef DEBUG4
2214 }
2215 
2216 void Sampling::fillLRRepository(void)
2217 {
2218  Matrix2D<double> L(4, 4), R(4, 4);
2219  Matrix2D<double> Identity(3,3);
2220  Identity.initIdentity();
2221  //NEXT 2 ROB
2222  R_repository.clear();
2223  L_repository.clear();
2224  R_repository.push_back(Identity);
2225  L_repository.push_back(Identity);
2226  for (int isym = 0; isym < SL.symsNo(); isym++)
2227  {
2228  SL.getMatrices(isym, L, R);
2229  R.resize(3, 3);
2230  L.resize(3, 3);
2231  R_repository.push_back(R);
2232  L_repository.push_back(L);
2233  }
2234 //#define DEBUG3
2235 #ifdef DEBUG3
2236  std::cout << "==============================\n" ;
2237  for (int isym = 0; isym < R_repository.size(); isym++)
2238  {
2239  std::cout << R_repository[isym];
2240  //std::cout << L_repository[isym];
2241  }
2242  std::cout << "-------------------------------\n";
2243 #endif
2244 #undef DEBUG3
2245 }
2246 
2247 void Sampling::fillExpDataProjectionDirectionByLR(
2248  const FileName &FnexperimentalImages)
2249 {
2250  //read input files
2251  MetaDataVec DFi;
2252  DFi.read(FnexperimentalImages);//experimental points
2253  fillExpDataProjectionDirectionByLR(DFi);
2254 }
2255 
2256 void Sampling::fillExpDataProjectionDirectionByLR(const MetaData &DFi)
2257 {
2258  std::vector <Matrix1D<double> > exp_data_projection_direction;
2259  Matrix1D<double> direction(3);
2260  DFi.firstRowId();
2261  //#define CHIMERA
2262 #ifdef CHIMERA
2263 
2264  std::ofstream filestr;
2265  filestr.open ("exp_data_projection_direction_by_L_R.bild");
2266  filestr << ".color white"
2267  << std::endl
2268  << ".sphere 0 0 0 .95"
2269  << std::endl
2270  ;
2271  filestr << ".color green"
2272  << std::endl
2273  ;
2274 #endif
2275 
2276  double img_tilt,img_rot,img_psi;
2277  FileName imgName;
2278  exp_data_fileNames.clear();
2279  for (size_t objId : DFi.ids())
2280  {
2281  DFi.getValue(MDL_ANGLE_ROT,img_rot, objId);
2282  DFi.getValue(MDL_ANGLE_TILT,img_tilt, objId);
2283  DFi.getValue(MDL_ANGLE_PSI,img_psi, objId);
2284  Euler_direction(img_rot, img_tilt, img_psi, direction);
2285  exp_data_projection_direction.push_back(direction);
2286  DFi.getValue(MDL_IMAGE,imgName, objId);
2287  exp_data_fileNames.push_back(imgName);
2288  }
2289 
2290  exp_data_projection_direction_by_L_R.clear();
2291 #ifdef MYPSI
2292 
2293  exp_data_projection_direction_by_L_R_psi.clear();
2294 #endif
2295 
2296  for (size_t i = 0; i < exp_data_projection_direction.size(); i++)
2297  for (size_t j = 0; j < R_repository.size(); j++)
2298  {
2299  direction = L_repository[j] *
2300  (exp_data_projection_direction[i].transpose() *
2301  R_repository[j]).transpose();
2302  exp_data_projection_direction_by_L_R.push_back(direction);
2303 #ifdef MYPSI
2304 
2305  Euler_apply_transf(L_repository[j],
2306  R_repository[j], img_rot,
2307  img_tilt,
2308  img_psi,
2309  rotp,
2310  tiltp,
2311  psip);
2312  exp_data_projection_direction_by_L_R_psi.push_back(psip);
2313 #endif
2314  //#define CHIMERA
2315 #ifdef CHIMERA
2316 
2317  filestr << ".sphere " << direction.transpose()
2318  << " 0.02" << std::endl
2319  ;
2320 #endif
2321 
2322  }
2323  //#define CHIMERA
2324 #ifdef CHIMERA
2325  filestr.close();
2326 #endif
2327  #undef CHIMERA
2328 }
pg_I2H
#define pg_I2H
Definition: symmetries.h:98
pg_TD
#define pg_TD
Definition: symmetries.h:84
REMOVE_LAST
#define REMOVE_LAST(vector)
Definition: sampling.cpp:1926
init_progress_bar
void init_progress_bar(long total)
Definition: xmipp_funcs.cpp:782
MDL_NEIGHBOR
Particular neighbor (pointed myNEIGHBORS)
Definition: metadata_label.h:308
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
Sampling::my_exp_img_per_sampling_point
std::vector< std::vector< size_t > > my_exp_img_per_sampling_point
Definition: sampling.h:90
Sampling::setSampling
void setSampling(double sampling)
Definition: sampling.cpp:121
Sampling::removeRedundantPoints
void removeRedundantPoints(const int symmetry, int sym_order)
Definition: sampling.cpp:691
pg_DN
#define pg_DN
Definition: symmetries.h:80
XMIPP_MAX
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
Sampling::setNoise
void setNoise(double deviation, int my_seed=-1)
Definition: sampling.cpp:134
MDL_ORDER
Definition: metadata_label.h:323
Sampling::vertices_vectors
std::vector< Matrix1D< double > > vertices_vectors
Definition: sampling.h:66
pg_T
#define pg_T
Definition: symmetries.h:83
pg_DNH
#define pg_DNH
Definition: symmetries.h:82
SymList::getMatrices
void getMatrices(int i, Matrix2D< double > &L, Matrix2D< double > &R, bool homogeneous=true) const
Definition: symmetries.cpp:348
MetaDataVec::read
void read(const FileName &inFile, const std::vector< MDLabel > *desiredLabels=nullptr, bool decomposeStack=true) override
Definition: metadata_vec.cpp:120
pg_O
#define pg_O
Definition: symmetries.h:86
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
Sampling::symmetry_file
FileName symmetry_file
Definition: sampling.h:135
Euler_angles2matrix
void Euler_angles2matrix(T alpha, T beta, T gamma, Matrix2D< T > &A, bool homogeneous)
Definition: geometry.cpp:624
Euler_direction
void Euler_direction(double alpha, double beta, double gamma, Matrix1D< double > &v)
Definition: geometry.cpp:721
DEG2RAD
#define DEG2RAD(d)
Definition: xmipp_macros.h:312
Sampling::createAsymUnitFile
void createAsymUnitFile(const FileName &docfilename)
Definition: sampling.cpp:1440
pg_I3H
#define pg_I3H
Definition: symmetries.h:99
metadata_vec.h
MD_APPEND
Definition: metadata_writemode.h:34
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
MDL_Z
Z component (double)
Definition: metadata_label.h:503
Sampling::sampling_points_angles
std::vector< Matrix1D< double > > sampling_points_angles
Definition: sampling.h:103
Sampling::setNeighborhoodRadius
void setNeighborhoodRadius(double neighborhood)
Definition: sampling.cpp:140
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
MetaData::firstRowId
virtual size_t firstRowId() const =0
pg_I3
#define pg_I3
Definition: symmetries.h:93
Sampling::R_repository
std::vector< Matrix2D< double > > R_repository
Definition: sampling.h:105
Sampling::neighborhood_radius_rad
double neighborhood_radius_rad
Definition: sampling.h:81
Sampling::findClosestSamplingPoint
void findClosestSamplingPoint(const MetaData &DFi, const FileName &output_file_root)
Definition: sampling.cpp:1991
Sampling::verbose
int verbose
Definition: sampling.h:126
Sampling::number_of_samples
size_t number_of_samples
Definition: sampling.h:75
MDRowVec::setValue
void setValue(const MDObject &object) override
Definition: metadata_row_vec.cpp:199
SymList::readSymmetryFile
int readSymmetryFile(FileName fn_sym, double accuracy=SYM_ACCURACY)
Definition: symmetries.cpp:33
MDL_SHIFT_X
Shift for the image in the X axis (double)
Definition: metadata_label.h:431
vectorR3
Matrix1D< double > vectorR3(double x, double y, double z)
Definition: matrix1d.cpp:892
FN_SAMPLING_PROJ
#define FN_SAMPLING_PROJ(base)
Definition: sampling.cpp:1491
Sampling::findClosestExperimentalPoint
void findClosestExperimentalPoint()
Definition: sampling.cpp:2100
pg_I1H
#define pg_I1H
Definition: symmetries.h:97
Sampling::no_redundant_sampling_points_vector
std::vector< Matrix1D< double > > no_redundant_sampling_points_vector
Definition: sampling.h:118
pg_I4H
#define pg_I4H
Definition: symmetries.h:100
Sampling::computeNeighbors
void computeNeighbors(bool only_winner=false)
Definition: sampling.cpp:1715
MDL_ANGLE_PSI
Special label to be used when gathering MDs in MpiMetadataPrograms.
Definition: metadata_label.h:51
pg_CNV
#define pg_CNV
Definition: symmetries.h:77
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
Matrix1D::selfNormalize
void selfNormalize()
Definition: matrix1d.cpp:723
FN_SAMPLING_SPHERE
#define FN_SAMPLING_SPHERE(base)
Definition: sampling.cpp:1493
Sampling::removeRedundantPointsExhaustive
void removeRedundantPointsExhaustive(const int symmetry, int sym_order, bool only_half_sphere, double max_ang)
Definition: sampling.cpp:1253
SymList::symsNo
int symsNo() const
Definition: symmetries.h:268
pg_OH
#define pg_OH
Definition: symmetries.h:87
Sampling::saveSamplingFile
void saveSamplingFile(const FileName &fn_base, bool write_vectors=true, bool write_sampling_sphere=false)
Definition: sampling.cpp:1495
MetaData::ids
virtual IdIteratorProxy< false > ids()
Definition: metadata_base.h:754
floatToString
String floatToString(float F, int _width, int _prec)
Definition: xmipp_strings.cpp:204
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
Euler_apply_transf
void Euler_apply_transf(const Matrix2D< double > &L, const Matrix2D< double > &R, double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
Definition: geometry.cpp:1038
MetaData::getValue
virtual bool getValue(MDObject &mdValueOut, size_t id) const =0
MetaDataVec::size
size_t size() const override
Definition: metadata_vec.cpp:393
i
#define i
Definition: numerical_recipes.cpp:2493
init_random_generator
void init_random_generator(int seed)
Definition: xmipp_funcs.cpp:394
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
MetaData::setComment
virtual void setComment(const String &newComment="No comment")
Definition: metadata_base.h:315
MetaDataVec::addRow
size_t addRow(const MDRow &row) override
Definition: metadata_vec.cpp:225
Matrix1D::transpose
Matrix1D< T > transpose() const
Definition: matrix1d.cpp:644
XMIPP_EQUAL_REAL
#define XMIPP_EQUAL_REAL(x, y)
Definition: xmipp_macros.h:122
Sampling::sampling_rate_rad
double sampling_rate_rad
Definition: sampling.h:69
pg_CS
#define pg_CS
Definition: symmetries.h:75
Sampling::exp_data_projection_direction_by_L_R
std::vector< Matrix1D< double > > exp_data_projection_direction_by_L_R
Definition: sampling.h:108
Sampling::operator==
bool operator==(const Sampling &op) const
Definition: sampling.cpp:107
FOR_ALL_ELEMENTS_IN_MATRIX1D
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
pg_I5H
#define pg_I5H
Definition: symmetries.h:101
Sampling::no_redundant_sampling_points_index
std::vector< size_t > no_redundant_sampling_points_index
Definition: sampling.h:123
Sampling::computeSamplingPoints
void computeSamplingPoints(bool only_half_sphere=true, double max_tilt=180, double min_tilt=0)
Definition: sampling.cpp:155
pg_I4
#define pg_I4
Definition: symmetries.h:94
XX
#define XX(v)
Definition: matrix1d.h:85
Sampling::numberSamplesAsymmetricUnit
size_t numberSamplesAsymmetricUnit
Definition: sampling.h:78
transpose
void transpose(double *array, int size)
Definition: primitive_FST.cpp:121
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
ERR_ARG_INCORRECT
Incorrect argument received.
Definition: xmipp_error.h:113
MetaDataVec::firstRowId
size_t firstRowId() const override
Definition: metadata_vec.cpp:508
pg_I1
#define pg_I1
Definition: symmetries.h:91
pg_DNV
#define pg_DNV
Definition: symmetries.h:81
Sampling::vertices_angles
Vect_angles vertices_angles
Definition: sampling.h:62
MD_OVERWRITE
Definition: metadata_writemode.h:33
sampling
#define sampling
Definition: spherical_harmonics.h:36
progress_bar
void progress_bar(long rlen)
Definition: xmipp_funcs.cpp:791
Sampling::cos_neighborhood_radius
double cos_neighborhood_radius
Definition: sampling.h:84
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
ROUND
#define ROUND(x)
Definition: xmipp_macros.h:210
Sampling::Sampling
Sampling()
Definition: sampling.cpp:32
direction
void direction(const MultidimArray< double > &orMap, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &dirMap, int x, int y)
Definition: fringe_processing.cpp:436
FN_SAMPLING_NEI
#define FN_SAMPLING_NEI(base)
Definition: sampling.cpp:1490
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
Sampling::fillLRRepository
void fillLRRepository(void)
Definition: sampling.cpp:2216
Sampling::fillDistance
void fillDistance(const Matrix1D< double > &starting_point, const Matrix1D< double > &ending_point, std::vector< Matrix1D< double > > &edge_vector, int number, bool only_half_spheree, double min_z=-10., double max_z=+10.)
Definition: sampling.cpp:631
FN_SAMPLING_EXTRA
#define FN_SAMPLING_EXTRA(base)
Definition: sampling.cpp:1492
Sampling::my_neighbors
std::vector< std::vector< size_t > > my_neighbors
Definition: sampling.h:87
pg_CNH
#define pg_CNH
Definition: symmetries.h:78
j
#define j
Definition: numerical_recipes.cpp:2493
MDL_POINTSASYMETRICUNIT
Definition: metadata_label.h:353
YY
#define YY(v)
Definition: matrix1d.h:93
MetaDataVec::getValue
bool getValue(MDObject &mdValueOut, size_t id) const override
Definition: metadata_vec.cpp:293
pg_TH
#define pg_TH
Definition: symmetries.h:85
MDL_REF
Class to which the image belongs (int)
Definition: metadata_label.h:377
MDL_SAMPLINGRATE
rfactor
Definition: metadata_label.h:398
pg_CN
#define pg_CN
Definition: symmetries.h:76
Sampling::readSamplingFile
void readSamplingFile(const FileName &infilename, bool read_vectors=true, bool read_sampling_sphere=false)
Definition: sampling.cpp:1592
Sampling::fillEdge
void fillEdge(const Matrix1D< double > &starting_point, const Matrix1D< double > &ending_point, std::vector< Matrix1D< double > > &edge_vector, bool FLAG_END)
Definition: sampling.cpp:606
Sampling::no_redundant_sampling_points_angles
std::vector< Matrix1D< double > > no_redundant_sampling_points_angles
Definition: sampling.h:121
MetaData::setColumnFormat
virtual void setColumnFormat(bool column)
Definition: metadata_base.h:293
MDL_NEIGHBORS
Vector of indexes to points some "neighbors".
Definition: metadata_label.h:307
Sampling::L_repository
std::vector< Matrix2D< double > > L_repository
Definition: sampling.h:106
dotProduct
T dotProduct(const Matrix1D< T > &v1, const Matrix1D< T > &v2)
Definition: matrix1d.h:1140
MDL_Y
Y component (double)
Definition: metadata_label.h:499
rnd_gaus
double rnd_gaus()
Definition: xmipp_funcs.cpp:466
cte_w
#define cte_w
Definition: sampling.h:35
pg_I2
#define pg_I2
Definition: symmetries.h:92
Sampling::sortFunc
int sortFunc(const Matrix1D< double > &a, const Matrix1D< double > &b)
Definition: sampling.cpp:586
MDL_X
X component (double)
Definition: metadata_label.h:495
Matrix1D::resizeNoCopy
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
pg_SN
#define pg_SN
Definition: symmetries.h:79
pg_CI
#define pg_CI
Definition: symmetries.h:74
pg_IH
#define pg_IH
Definition: symmetries.h:89
FIRST_IMAGE
#define FIRST_IMAGE
Definition: xmipp_image_macros.h:38
Sampling::sampling_points_vector
std::vector< Matrix1D< double > > sampling_points_vector
Definition: sampling.h:101
MDL_SHIFT_Y
Shift for the image in the Y axis (double)
Definition: metadata_label.h:435
MDL_NEIGHBORHOOD_RADIUS
Radius of the neighborhood (radians)
Definition: metadata_label.h:309
Sampling::exp_data_fileNames
std::vector< FileName > exp_data_fileNames
Definition: sampling.h:110
Sampling::fillExpDataProjectionDirectionByLR
void fillExpDataProjectionDirectionByLR(const MetaData &DFi)
Definition: sampling.cpp:2256
MetaDataVec::containsLabel
bool containsLabel(const MDLabel label) const override
Definition: metadata_vec.cpp:397
Sampling::removePointsFarAwayFromExperimentalData
void removePointsFarAwayFromExperimentalData()
Definition: sampling.cpp:1928
PI
#define PI
Definition: tools.h:43
Sampling::sampling_noise
double sampling_noise
Definition: sampling.h:72
Sampling::createSymFile
void createSymFile(const FileName &simFp, int symmetry, int sym_order)
Definition: sampling.cpp:1319
n
int * n
Definition: numerical_recipes.cpp:2229
MDL_IMAGE
Name of an image (std::string)
Definition: metadata_label.h:246
Matrix2D::resize
void resize(size_t Ydim, size_t Xdim, bool noCopy=false)
Definition: matrix2d.cpp:1022
ZZ
#define ZZ(v)
Definition: matrix1d.h:101
a
doublereal * a
Definition: numerical_recipes.cpp:2230
Sampling::SL
SymList SL
Definition: sampling.h:138
CLEAR_VECTORS
#define CLEAR_VECTORS()
Definition: sampling.cpp:679
Matrix2D::initIdentity
void initIdentity()
Definition: matrix2d.h:673
pg_I
#define pg_I
Definition: symmetries.h:88
xmipp_image_macros.h
CREATE_INDEXES
#define CREATE_INDEXES()
Definition: sampling.cpp:684
MDL_STAR_COMMENT
A comment for this object /*** NOTE THIS IS A SPECIAL CASE AND SO IS TREATED ***/.
Definition: metadata_label.h:453
pg_I5
#define pg_I5
Definition: symmetries.h:95
Generated by   doxygen 1.8.13