# A.I.M.S algorithms

minimisation_d.h
Go to the documentation of this file.
1 /* Copyright (C) 2000-2013 CEA
2  *
3  * This software and supporting documentation were developed by
4  * bioPICSEL
5  * CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
6  * 18, route du Panorama
7  * 92265 Fontenay-aux-Roses
8  * France
9  */
10
11
12 #ifndef AIMS_REGISTRATION_MINIMISATION_D_H
13 #define AIMS_REGISTRATION_MINIMISATION_D_H
14
15
18 #include <iostream>
19
20 #include <aims/io/motionW.h>
23
24
25 //methode principale:
26
27 template <class T>
29 {
30
31  //on charge l'ensemble des points a apparier
32  std::vector<Point3df> ptstest = displacementField.getpointstest();
33  std::vector<Point3df> ptsref = displacementField.getpointsref();
34  int size = ptstest.size();
35  unsigned sizekept = unsigned ( double(ptstest.size())*_Pkept );
36
37  Motion m;
38
39  if ((sizekept > 0) && (ptsref.size() > 0))
40  {
41  // Calcul de la transformation LS choisie
42  switch(_transfo)
43  {
44  case 1:
45  {
46  // std::cout<<"Calcul transfo rigide..."<<std::endl;
47  aims::RigidLeastSquareEstimation leastSquareEstimation(ptstest, ptsref);
48  m = *(leastSquareEstimation.motion());
49  break;
50  }
51
52  case 2:
53  {
54  // std::cout<<"Calcul transfo similitude..."<<std::endl;
55  aims::SimiLeastSquareEstimation leastSquareEstimation(ptstest, ptsref);
56  m = *(leastSquareEstimation.motion());
57  break;
58  }
59
60  case 3:
61  {
62  // std::cout<<"Calcul transfo affine..."<<std::endl;
63  aims::AffineLeastSquareEstimation leastSquareEstimation(ptstest, ptsref);
64  m = *(leastSquareEstimation.motion());
65  break;
66  }
67
68  default:
69  {
70  m.setToIdentity();
71  throw std::runtime_error("The transformation chosen is not valid !") ;
72  break;
73  }
74  }
75
76  //ROBUSTESSE DEBUT
77
78  // Le tri des residus pour ref et test
79  std::multimap<float, Point3df> mt;
80  std::multimap<float, Point3df> mr;
81
82  std::multimap<float, Point3df>::iterator itt;
83  std::multimap<float, Point3df>::iterator itr;
84
85  // Diverses variables utiles
86  Point3df r; // point residu courant
87  float r2; // norme carree de r
88  unsigned indice; // indice de boucle multiple (variable locale boucle impossible)
89  Motion rigide_kept;
90  Motion similitude_kept;
91
92  // La liste des points ref et test ayant les plus petits residus
93  std::vector<Point3df> ptstestkept;
94  std::vector<Point3df> ptsrefkept;
95
96
97  for(int boucle=0; boucle < 3; boucle++)
98  {
99
100  //std::cout << "Minimisation - s : " << size << ", m : " << m << std::endl << std::flush;
101
102  mt.clear(); mr.clear();
103  ptstestkept.clear(); ptsrefkept.clear();
104
105  for(int l = 0; l < size; l++)
106  {
107  Point3df x = ptstest[l];
108  Point3df y = ptsref[l];
109
110  r = y - m.transform(x); //calcul du residu pour chaque point
111  r2 = r.norm(); // sa norme carree
112  mt.insert(std::pair<float, Point3df>(r2, x) ); //insertion du residu et du point test dans map
113  mr.insert(std::pair<float, Point3df>(r2, y) ); //insertion du residu et du point test dans map
114  }
115
116  // On remplit les listes avec une partie seulement des meilleurs residus
117  for (itt = mt.begin(),itr = mr.begin(), indice = 0 ;indice<sizekept;++itt,++itr,++indice)
118  {
119  ptstestkept.push_back( (*itt).second );
120  ptsrefkept.push_back( (*itr).second );
121  }
122
123  // Calcul de la transformation LS choisie
124  switch(_transfo)
125  {
126  case 1:
127  {
128  // cout<<"Calcul transfo rigide..."<<cout;
129  aims::RigidLeastSquareEstimation leastSquareEstimation_kept(ptstestkept, ptsrefkept);
130  m = *(leastSquareEstimation_kept.motion());
131  break;
132  }
133
134  case 2:
135  {
136  // cout<<"Calcul transfo similitude..."<<cout;
137  aims::SimiLeastSquareEstimation leastSquareEstimation_kept(ptstestkept, ptsrefkept);
138  m = *(leastSquareEstimation_kept.motion());
139  break;
140  }
141
142  case 3:
143  {
144  // cout<<"Calcul transfo affine..."<<cout;
145  aims::AffineLeastSquareEstimation leastSquareEstimation_kept(ptstestkept, ptsrefkept);
146  m = *(leastSquareEstimation_kept.motion());
147  break;
148  }
149
150  default:
151  {
152  m.setToIdentity();
153  throw std::runtime_error("The transformation chosen is not valid !") ;
154  break;
155  }
156  }
157  } // fin boucle robustesse
158  }
159  else
160  {
161  cartoDbgMsg( 1, "Not enough information in reference or test image to estimate transformation" );
162  }
163
164  return (m);
165 }
166
167
168 #endif
169
std::vector< Point3df > getpointstest()
std::vector< Point3df > getpointsref()
Motion quaternion(DisplacementField< T > &displacementField)
Point3dd transform(double x, double y, double z) const
float norm() const
#define cartoDbgMsg(level, message)