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)