Xmipp  v3.23.11-Nereus
Public Member Functions | Public Attributes | List of all members
EulerSolver Class Reference
Libraries » Reconstruction library » angular_common_lines

#include <angular_commonline.h>

Inheritance diagram for EulerSolver:
Inheritance graph
[legend]
Collaboration diagram for EulerSolver:
Collaboration graph
[legend]

Public Member Functions

 EulerSolver (int dim, int pop, const Matrix1D< int > &newAlreadyOptimized, const Matrix1D< double > &newCurrentSolution, const MultidimArray< int > &newImgIdx, const Prog_Angular_CommonLine *newParent)
 
double EnergyFunction (double trial[], bool &bAtSolution)
 
double similarityBetweenTwoLines (int imgi, int imgj)
 
void setShow (bool newShow)
 
- Public Member Functions inherited from DESolver
 DESolver (int dim, int popSize)
 Empty constructor. More...
 
 DESolver (const DESolver &)=delete
 
DESolver & operator= (const DESolver &)=delete
 
virtual ~DESolver (void)
 Destructor. More...
 
void Setup (double min[], double max[], int deStrategy, double diffScale, double crossoverProb)
 Setup() must be called before solve to set min, max, strategy etc. More...
 
virtual bool Solve (int maxGenerations)
 
int Dimension () const
 Return dimension. More...
 
int Population () const
 Return population. More...
 
double Energy () const
 Call these functions after Solve() to get results. More...
 
double * Solution (void)
 Return best solution. More...
 
int Generations () const
 Return the number of generations. More...
 

Public Attributes

bool show
 
int NToSolve
 
int Ndim
 
int Nimg
 
double roti
 
double tilti
 
double psii
 
double rotj
 
double tiltj
 
double psij
 
Matrix1D< double > normali
 
Matrix1D< double > normalj
 
Matrix1D< double > commonline
 
Matrix1D< double > commonlinei
 
Matrix1D< double > commonlinej
 
const Prog_Angular_CommonLine * parent
 
const Matrix1D< int > * alreadyOptimized
 
const Matrix1D< double > * currentSolution
 
const MultidimArray< int > * imgIdx
 
Matrix1D< double > imgAvgCorrelation
 
Matrix1D< double > imgMinCorrelation
 
MultidimArray< double > correlationMatrix
 

Additional Inherited Members

- Protected Member Functions inherited from DESolver
void SelectSamples (int candidate, int *r1, int *r2=nullptr, int *r3=nullptr, int *r4=nullptr, int *r5=nullptr)
 
- Protected Attributes inherited from DESolver
int nDim
 
int nPop
 
int generations
 
int strategy
 
StrategyFunction calcTrialSolution
 
double scale
 
double probability
 
double trialEnergy
 
double bestEnergy
 
double * trialSolution
 
double * bestSolution
 
double * popEnergy
 
double * population
 

Detailed Description

Definition at line 45 of file angular_commonline.h.

Constructor & Destructor Documentation

◆ EulerSolver()

EulerSolver::EulerSolver ( int  dim,
int  pop,
const Matrix1D< int > &  newAlreadyOptimized,
const Matrix1D< double > &  newCurrentSolution,
const MultidimArray< int > &  newImgIdx,
const Prog_Angular_CommonLine *  newParent 
)

Definition at line 33 of file angular_commonline.cpp.

37  : DESolver(dim, pop)
38 {
39  parent = newParent;
40  Ndim = dim;
41  NToSolve=dim/3;
42  Nimg=VEC_XSIZE(newAlreadyOptimized);
43  imgAvgCorrelation.initZeros(Nimg);
44  imgMinCorrelation.initZeros(Nimg);
45  correlationMatrix.initZeros(Nimg,Nimg);
46  commonline.initZeros(3);
47  alreadyOptimized=&newAlreadyOptimized;
48  currentSolution=&newCurrentSolution;
49  imgIdx=&newImgIdx;
50  show=false;
51 }
EulerSolver::parent
const Prog_Angular_CommonLine * parent
Definition: angular_commonline.h:52
EulerSolver::imgAvgCorrelation
Matrix1D< double > imgAvgCorrelation
Definition: angular_commonline.h:57
VEC_XSIZE
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
EulerSolver::imgIdx
const MultidimArray< int > * imgIdx
Definition: angular_commonline.h:55
EulerSolver::currentSolution
const Matrix1D< double > * currentSolution
Definition: angular_commonline.h:54
DESolver::DESolver
DESolver(int dim, int popSize)
Empty constructor.
Definition: numerical_tools.cpp:387
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
EulerSolver::commonline
Matrix1D< double > commonline
Definition: angular_commonline.h:51
EulerSolver::imgMinCorrelation
Matrix1D< double > imgMinCorrelation
Definition: angular_commonline.h:58
EulerSolver::alreadyOptimized
const Matrix1D< int > * alreadyOptimized
Definition: angular_commonline.h:53
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
EulerSolver::correlationMatrix
MultidimArray< double > correlationMatrix
Definition: angular_commonline.h:59
pop
void pop(struct stack_T *stack)
Definition: stack.cpp:59

Member Function Documentation

◆ EnergyFunction()

double EulerSolver::EnergyFunction ( double  testSolution[],
bool &  bAtSolution 
)
virtual

EnergyFunction must be overridden for problem to solve testSolution[] is nDim array for a candidate solution setting bAtSolution = true indicates solution is found and Solve() immediately returns true.

Implements DESolver.

Definition at line 59 of file angular_commonline.cpp.

60 {
61  // Check limits
62  int i,j;
63 
64  for (i=0, j=0; i<Ndim; i++,j=(j+1)%3)
65  {
66  if (j==1)
67  trial[i]=realWRAP(trial[i],0,180);
68  else
69  trial[i]=realWRAP(trial[i],0,360);
70  }
71 
72  // Get number of symmetries
73  int Nsym=parent->L.size();
74 
75  // Evaluate goal function for this trial
76  double retval=0;
77  double Ncomparisons=0;
78  double worseSimilarity=2;
79  int idx;
80  static Matrix1D<int> imgCorrelationN;
81  static double minval=0;
82  imgAvgCorrelation.initZeros(Nimg);
83  imgMinCorrelation.initZeros(Nimg);
84  imgMinCorrelation.initConstant(2);
85  correlationMatrix.resize(Nimg,Nimg);
86  imgCorrelationN.initZeros(Nimg);
87  for (int imgi=0; imgi<Nimg; imgi++)
88  {
89  // Get the right angles for this image
90  if ((*alreadyOptimized)(imgi)==0)
91  continue;
92  else if ((*alreadyOptimized)(imgi)==1)
93  {
94  for (idx=0; idx<NToSolve; idx++)
95  if ((*imgIdx)(idx)==imgi)
96  break;
97  idx*=3;
98  /*
99  trial[idx] = parent->initialSolution(3*imgi);
100  trial[idx+1] = parent->initialSolution(3*imgi+1);
101  trial[idx+2] = parent->initialSolution(3*imgi+2);
102  */
103  roti = trial[idx++];
104  tilti = trial[idx++];
105  psii = trial[idx];
106  }
107  else
108  {
109  idx=3*imgi;
110  roti = (*currentSolution)(idx++);
111  tilti = (*currentSolution)(idx++);
112  psii = (*currentSolution)(idx);
113  }
114 
115  // Loop for the second image
116  for (int imgj=imgi+1; imgj<Nimg; imgj++)
117  {
118  // Get the right angles for this image
119  if ((*alreadyOptimized)(imgj)==0)
120  continue;
121  else if ((*alreadyOptimized)(imgj)==1)
122  {
123  for (idx=0; idx<NToSolve; idx++)
124  if ((*imgIdx)(idx)==imgj)
125  break;
126  idx*=3;
127  /*
128  trial[idx] = parent->initialSolution(3*imgj);
129  trial[idx+1] = parent->initialSolution(3*imgj+1);
130  trial[idx+2] = parent->initialSolution(3*imgj+2);
131  */
132  rotj = trial[idx++];
133  tiltj = trial[idx++];
134  psij = trial[idx];
135  }
136  else
137  {
138  idx=3*imgj;
139  rotj = (*currentSolution)(idx++);
140  tiltj = (*currentSolution)(idx++);
141  psij = (*currentSolution)(idx);
142  }
143 
144  // Check that at least one of the two images is new
145  if ((*alreadyOptimized)(imgi)==2 && (*alreadyOptimized)(imgj)==2)
146  continue;
147 
148  double similarity=similarityBetweenTwoLines(imgi,imgj);
149  if (similarity>0)
150  {
151  retval -= similarity;
152  Ncomparisons++;
153  worseSimilarity=XMIPP_MIN(worseSimilarity,similarity);
154  imgAvgCorrelation(imgi)+=similarity;
155  imgAvgCorrelation(imgj)+=similarity;
156  imgMinCorrelation(imgi)=XMIPP_MIN(similarity,
157  imgMinCorrelation(imgi));
158  imgMinCorrelation(imgj)=XMIPP_MIN(similarity,
159  imgMinCorrelation(imgj));
160  correlationMatrix(imgi,imgj)=correlationMatrix(imgj,imgi)=
161  similarity;
162  imgCorrelationN(imgi)++;
163  imgCorrelationN(imgj)++;
164  }
165 
166  double original_rotj=rotj;
167  double original_tiltj=tiltj;
168  double original_psij=psij;
169 
170  for (int sym=0; sym<Nsym; sym++)
171  {
172  Euler_apply_transf(parent->L[sym],parent->R[sym],
173  original_rotj,original_tiltj,original_psij,
174  rotj, tiltj, psij);
175  similarity=similarityBetweenTwoLines(imgi,imgj);
176  if (similarity>0)
177  {
178  retval -= similarity;
179  Ncomparisons++;
180  worseSimilarity=XMIPP_MIN(worseSimilarity,similarity);
181  imgAvgCorrelation(imgi)+=similarity;
182  imgAvgCorrelation(imgj)+=similarity;
183  imgMinCorrelation(imgi)=XMIPP_MIN(similarity,
184  imgMinCorrelation(imgi));
185  imgMinCorrelation(imgj)=XMIPP_MIN(similarity,
186  imgMinCorrelation(imgj));
187  correlationMatrix(imgi,imgj)=correlationMatrix(imgj,imgi)=
188  similarity;
189  imgCorrelationN(imgi)++;
190  imgCorrelationN(imgj)++;
191  }
192  }
193  }
194  }
195  if (Ncomparisons>0)
196  retval/=Ncomparisons;
197  FOR_ALL_ELEMENTS_IN_MATRIX1D(imgAvgCorrelation)
198  if (imgCorrelationN(i)!=0)
199  imgAvgCorrelation(i)/=imgCorrelationN(i);
200  if (show)
201  std::cout
202  << "Average Distance = " << retval << std::endl
203  << "Worse Distance = " << -worseSimilarity << std::endl
204  << "imgAvgCorrelation = " << imgAvgCorrelation.transpose() << std::endl
205  << "imgMinCorrelation = " << imgMinCorrelation.transpose() << std::endl
206  << std::endl
207  << std::endl
208  ;
209  if (retval<minval)
210  {
211  minval=retval;
212  std::cout << "MinEnergy=" << minval << std::endl;
213  }
214  return retval;
215  // return 0.5*(retval-worseSimilarity);
216 }
Prog_Angular_CommonLine::R
std::vector< Matrix2D< double > > R
Definition: angular_commonline.h:181
EulerSolver::parent
const Prog_Angular_CommonLine * parent
Definition: angular_commonline.h:52
MultidimArray::resize
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
Definition: multidim_array.cpp:983
EulerSolver::imgAvgCorrelation
Matrix1D< double > imgAvgCorrelation
Definition: angular_commonline.h:57
Prog_Angular_CommonLine::L
std::vector< Matrix2D< double > > L
Definition: angular_commonline.h:178
EulerSolver::imgIdx
const MultidimArray< int > * imgIdx
Definition: angular_commonline.h:55
EulerSolver::similarityBetweenTwoLines
double similarityBetweenTwoLines(int imgi, int imgj)
Definition: angular_commonline.cpp:220
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
i
#define i
Definition: numerical_recipes.cpp:2493
Matrix1D::transpose
Matrix1D< T > transpose() const
Definition: matrix1d.cpp:644
FOR_ALL_ELEMENTS_IN_MATRIX1D
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
j
#define j
Definition: numerical_recipes.cpp:2493
EulerSolver::imgMinCorrelation
Matrix1D< double > imgMinCorrelation
Definition: angular_commonline.h:58
EulerSolver::alreadyOptimized
const Matrix1D< int > * alreadyOptimized
Definition: angular_commonline.h:53
realWRAP
#define realWRAP(x, x0, xF)
Definition: xmipp_macros.h:304
Matrix1D::initConstant
void initConstant(T val)
Definition: matrix1d.cpp:83
EulerSolver::correlationMatrix
MultidimArray< double > correlationMatrix
Definition: angular_commonline.h:59

◆ setShow()

void EulerSolver::setShow ( bool  newShow)

Definition at line 53 of file angular_commonline.cpp.

54 {
55  show=newShow;
56 }

◆ similarityBetweenTwoLines()

double EulerSolver::similarityBetweenTwoLines ( int  imgi,
int  imgj 
)

Definition at line 220 of file angular_commonline.cpp.

221 {
222  // Compute the direction of the common line in the
223  // universal coordinate system
224  Euler_direction(roti, tilti, psii, normali);
225  Euler_direction(rotj, tiltj, psij, normalj);
226  vectorProduct(normali,normalj,commonline);
227  if (commonline.module()<XMIPP_EQUAL_ACCURACY)
228  {
229  // They do not share a common line but a common plane
230  if (show)
231  {
232  std::cout
233  << "imgi=" << imgi << " (rot,tilt,psi)=("
234  << roti << "," << tilti << "," << psii << ") normali="
235  << normali.transpose() << std::endl
236  << "imgj=" << imgj << " (rot,tilt,psi)=("
237  << rotj << "," << tiltj << "," << psij << ") normalj="
238  << normalj.transpose() << std::endl
239  ;
240  std::cout << "These two images are taken from the same direction\n";
241  std::cout << "Press any key\n";
242  char c;
243  std::cin >> c;
244  }
245  return -1;
246  }
247  commonline.selfNormalize();
248 
249  // Compute the direction of the common line in each
250  // of the projection coordinate systems
251  Uproject_to_plane(commonline, roti, tilti, psii, commonlinei);
252  Uproject_to_plane(commonline, rotj, tiltj, psij, commonlinej);
253 
254  // Compute the angle of the common line in i and j images
255  double angi=RAD2DEG(atan2(YY(commonlinei),XX(commonlinei)));
256  double angj=RAD2DEG(atan2(YY(commonlinej),XX(commonlinej)));
257 
258  auto idxAngi = (int)intWRAP(-((int)angi),0,359);
259  auto idxAngj = (int)intWRAP(-((int)angj),0,359);
260 
261  double retval1=
262  correlationIndex(parent->radon[imgi][idxAngi],
263  parent->radon[imgj][idxAngj]);
264  int idxBesti=parent->bestLine(imgi,imgj);
265  int idxBestj=parent->bestLine(imgj,imgi);
266  int retvali=ABS(idxBesti-idxAngi);
267  if (retvali>180)
268  retvali=ABS(retvali-360);
269  int retvalj=ABS(idxBestj-idxAngj);
270  if (retvalj>180)
271  retvalj=ABS(retvalj-360);
272  retvali=XMIPP_MIN(retvali,10);
273  retvalj=XMIPP_MIN(retvalj,10);
274 
275  const double i180=0.5*1.0/180.0;
276  double retval2=1-(retvali+retvalj)*i180;
277 
278  if (show)
279  {
280  std::cout
281  << "imgi=" << imgi << " (rot,tilt,psi)=("
282  << roti << "," << tilti << "," << psii << ") normali="
283  << normali.transpose() << std::endl
284  << "imgj=" << imgj << " (rot,tilt,psi)=("
285  << rotj << "," << tiltj << "," << psij << ") normalj="
286  << normalj.transpose() << std::endl
287  << "commonline= " << commonline.transpose() << std::endl
288  << "in imgi=" << commonlinei.transpose() << " anglei=" << angi
289  << " (" << idxAngi << ")\n"
290  << "in imgj=" << commonlinej.transpose() << " anglej=" << angj
291  << " (" << idxAngj << ")\n"
292  << "Distance between lines = " << retval1 << " " << retval2
293  << std::endl
294  ;
295  parent->radon[imgi][idxAngi].write("PPPradoni1.txt");
296  parent->radon[imgj][idxAngj].write("PPPradonj1.txt");
297  std::cout << "Press any key\n";
298  char c;
299  std::cin >> c;
300  }
301 
302 #ifdef DEBUG
303  for (idxAngi=0; idxAngi<360; idxAngi++)
304  {
305  double retval3=
306  correlationIndex(parent->radon[imgi][idxAngi],
307  parent->radon[imgj][idxAngj]);
308  ;
309 
310  if (show)
311  {
312  std::cout
313  << " (" << idxAngi << ")\n"
314  << "Distance between lines = " << retval3 << std::endl
315  << std::endl
316  ;
317  parent->radon[imgi][idxAngi].write("PPPradoni2.txt");
318  parent->radon[imgj][idxAngj].write("PPPradonj2.txt");
319  parent->radonDerivative[imgi][idxAngi].write("PPPradonDerivativei2.txt");
320  parent->radonDerivative[imgj][idxAngj].write("PPPradonDerivativej2.txt");
321  std::cout << "Press any key\n";
322  char c;
323  std::cin >> c;
324  if (c=='q')
325  break;
326  }
327  }
328 #endif
329  return retval1*retval2;
330 }
EulerSolver::parent
const Prog_Angular_CommonLine * parent
Definition: angular_commonline.h:52
Matrix1D::module
double module() const
Definition: matrix1d.h:983
c
doublereal * c
Definition: numerical_recipes.cpp:2230
Euler_direction
void Euler_direction(double alpha, double beta, double gamma, Matrix1D< double > &v)
Definition: geometry.cpp:721
correlationIndex
double correlationIndex(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)
Definition: multidim_array.h:3974
EulerSolver::commonlinei
Matrix1D< double > commonlinei
Definition: angular_commonline.h:51
EulerSolver::normali
Matrix1D< double > normali
Definition: angular_commonline.h:51
Matrix1D::selfNormalize
void selfNormalize()
Definition: matrix1d.cpp:723
vectorProduct
Matrix1D< T > vectorProduct(const Matrix1D< T > &v1, const Matrix1D< T > &v2)
Definition: matrix1d.h:1165
EulerSolver::normalj
Matrix1D< double > normalj
Definition: angular_commonline.h:51
Matrix1D::transpose
Matrix1D< T > transpose() const
Definition: matrix1d.cpp:644
XX
#define XX(v)
Definition: matrix1d.h:85
XMIPP_EQUAL_ACCURACY
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
Prog_Angular_CommonLine::radon
std::vector< std::vector< MultidimArray< double > > > radon
Definition: angular_commonline.h:175
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
EulerSolver::commonline
Matrix1D< double > commonline
Definition: angular_commonline.h:51
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
YY
#define YY(v)
Definition: matrix1d.h:93
Prog_Angular_CommonLine::bestLine
Matrix2D< int > bestLine
Definition: angular_commonline.h:166
RAD2DEG
#define RAD2DEG(r)
Definition: xmipp_macros.h:320
EulerSolver::commonlinej
Matrix1D< double > commonlinej
Definition: angular_commonline.h:51
Uproject_to_plane
void Uproject_to_plane(const Matrix1D< double > &point, const Matrix1D< double > &direction, double distance, Matrix1D< double > &result)
Definition: geometry.cpp:39
intWRAP
#define intWRAP(x, x0, xF)
Definition: xmipp_macros.h:272

Member Data Documentation

◆ alreadyOptimized

const Matrix1D<int>* EulerSolver::alreadyOptimized

Definition at line 53 of file angular_commonline.h.

◆ commonline

Matrix1D<double> EulerSolver::commonline

Definition at line 51 of file angular_commonline.h.

◆ commonlinei

Matrix1D<double> EulerSolver::commonlinei

Definition at line 51 of file angular_commonline.h.

◆ commonlinej

Matrix1D<double> EulerSolver::commonlinej

Definition at line 51 of file angular_commonline.h.

◆ correlationMatrix

MultidimArray<double> EulerSolver::correlationMatrix

Definition at line 59 of file angular_commonline.h.

◆ currentSolution

const Matrix1D<double>* EulerSolver::currentSolution

Definition at line 54 of file angular_commonline.h.

◆ imgAvgCorrelation

Matrix1D<double> EulerSolver::imgAvgCorrelation

Definition at line 57 of file angular_commonline.h.

◆ imgIdx

const MultidimArray<int>* EulerSolver::imgIdx

Definition at line 55 of file angular_commonline.h.

◆ imgMinCorrelation

Matrix1D<double> EulerSolver::imgMinCorrelation

Definition at line 58 of file angular_commonline.h.

◆ Ndim

int EulerSolver::Ndim

Definition at line 48 of file angular_commonline.h.

◆ Nimg

int EulerSolver::Nimg

Definition at line 48 of file angular_commonline.h.

◆ normali

Matrix1D<double> EulerSolver::normali

Definition at line 51 of file angular_commonline.h.

◆ normalj

Matrix1D<double> EulerSolver::normalj

Definition at line 51 of file angular_commonline.h.

◆ NToSolve

int EulerSolver::NToSolve

Definition at line 48 of file angular_commonline.h.

◆ parent

const Prog_Angular_CommonLine* EulerSolver::parent

Definition at line 52 of file angular_commonline.h.

◆ psii

double EulerSolver::psii

Definition at line 49 of file angular_commonline.h.

◆ psij

double EulerSolver::psij

Definition at line 50 of file angular_commonline.h.

◆ roti

double EulerSolver::roti

Definition at line 49 of file angular_commonline.h.

◆ rotj

double EulerSolver::rotj

Definition at line 50 of file angular_commonline.h.

◆ show

bool EulerSolver::show

Definition at line 47 of file angular_commonline.h.

◆ tilti

double EulerSolver::tilti

Definition at line 49 of file angular_commonline.h.

◆ tiltj

double EulerSolver::tiltj

Definition at line 50 of file angular_commonline.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.13