A.I.M.S algorithms


block_matching_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_BLOCK_MATCHING_D_H
13 #define AIMS_REGISTRATION_BLOCK_MATCHING_D_H
14 
15 #include <aims/data/data_g.h>
17 #include <aims/math/math_g.h>
23 #include <cartobase/type/limits.h>
24 #include <vector>
25 
26 
27 template<typename T>
29 {
30  _done=false;
31  _initialisation=Motion();
32  _level_start=-1;
33  _level_stop=-1;
34  _transfo=1;
35  _tailleBloc=Point3d(4,4,4);
36  _cutVar=-1;
37  _stopVar=-1;
38  _Pkept=0.5;
39  _seuilCorrel=-1;
40  _seuils[0]=std::numeric_limits<T>::min();
41  _seuils[1]=std::numeric_limits<T>::max();
42  _seuils[2]=_seuils[0];
43  _seuils[3]=_seuils[1];
44  _itermax=3;
45  _info=false;
46 }
47 
48 
49 template<typename T>
51 {
52 
53 
54  AimsData<T> test=test_orig.clone();
55 
56  // Mise de l'image test a la resolution de l'image de reference si besoin
57  if( (ref.sizeX()!=test.sizeX()) || (ref.sizeY()!=test.sizeY()) || (ref.sizeZ()!=test.sizeZ()) )
58  {
59  Motion identity;
60  identity.setToIdentity();
61  LinearResampler<T> reech;
62  reech.setRef( test );
63 
64  test = reech.doit( identity,
65  unsigned (test.dimX()*test.sizeX()/(1.0*ref.sizeX()) + .5),
66  unsigned (test.dimY()*test.sizeY()/(1.0*ref.sizeY()) + .5),
67  unsigned (test.dimZ()*test.sizeZ()/(1.0*ref.sizeZ()) + .5),
68  Point3df(ref.sizeX(),ref.sizeY(),ref.sizeZ()));
69  }
70 
71 
72  // Instanciation des differentes classes externes
73  ScaleControl scaleControl; // Gere la pyramide
74  DisplacementField<T> displacementField; // Calcule le champ d'appariemment
75  Transformation transformation; // Gere le critere de changement d'echelle delta
76  Minimisation minimisation(_Pkept, _transfo); // Calcule la transformation a partir du champ
77 
78 
79  //declaration de l'image transformee a chaque etape
80  AimsData<T> testtrans = test.clone();
81 
82 
83 
84  // Declaration des differents motions necessaires
85  // p transfo totale precedente, q transfo ponctuelle calculee, r nouvelle transfo totale composee
86  Motion p, q, r ;
87  p=_initialisation;
88  q.setToIdentity();
89  r.setToIdentity();
90 
91  {
92  LinearResampler<T> reech;
93  reech.setRef( test );
94 
95  testtrans = reech.doit( p,
96  testtrans.dimX(),
97  testtrans.dimY(),
98  testtrans.dimZ(),
99  Point3df(testtrans.sizeX(),testtrans.sizeY(),testtrans.sizeZ()));
100  }
101 
102 
103 
104  // Quelques initialisations en millimetres
105  transformation.setX(int(test.dimX()*test.sizeX()));
106  transformation.setY(int(test.dimY()*test.sizeY()));
107  transformation.setZ(int(test.dimZ()*test.sizeZ()));
108  transformation.setiterMax(_itermax); //on fixe le nombre maxi d'iter par niveau
109  transformation.setcx(0);
110  transformation.setcy(0);
111  transformation.setcz(0);
112  if(ref.dimZ()==1) _tailleBloc[2] = 1;
113  scaleControl.init( ref, _level_start, _level_stop, _cutVar, _stopVar, _seuilCorrel, _tailleBloc);
114  int count = 0 ;
115  // Initialisation de delta au carre de 4 fois un majorant de 2 fois la diagonale en mm...
116  float delta_init = 8*pow(test.dimX()*test.sizeX() + test.dimY()*test.sizeY() + test.dimZ()*test.sizeZ(), 2);
117 
118 
119  // Affichage de tous les parametres avant debut :
120  if(_info)
121  {
122  std::cout<<std::endl<<"PARAMETRES INITIAUX :"<<std::endl;
123  std::cout<<"Transfo initiale : "<<std::flush;
124  std::cout << p;
125 
126  std::cout<<"Transformation recherchee : ";
127  switch(_transfo){
128  case 1: std::cout<<"RIGIDE"<<std::endl; break;
129  case 2: std::cout<<"SIMILITUDE"<<std::endl; break;
130  case 3: std::cout<<"AFFINE"<<std::endl; break;}
131 
132  std::cout<<"Pyramide debut : "<<scaleControl.getScale()<<" fin : "<<scaleControl.getlevel_stop()<<std::endl;
133  std::cout<<"Taille des blocs : x = "<<_tailleBloc[0]<<" y = "<<_tailleBloc[1]<<" z = "<<_tailleBloc[2]<<std::endl;
134  std::cout<<"Pourcentage de blocs conserves par variance : %initial = "<<scaleControl.getcutVar()<<" %final ="<<scaleControl.getstopVar()<<std::endl;
135  std::cout<<"Seuil sur la mesure de similarite = "<<_seuilCorrel<<std::endl;
136  std::cout<<"Seuils sur les niveaux de gris : "<<std::endl;
137  std::cout<<" seuilBasRef = "<<_seuils[0]<<" seuilHautRef = "<<_seuils[1]<<std::endl;
138  std::cout<<" seuilBasTest = "<<_seuils[2]<<" seuilHautTest = "<<_seuils[3]<<std::endl;
139  std::cout<<"Nombre d'iterations par niveau de pyramide : "<<_itermax<<std::endl;
140  }
141 
142  // MISE DE p EN VOXELS !
143  p.translation()[0] /= test.sizeX();
144  p.translation()[1] /= test.sizeY();
145  p.translation()[2] /= test.sizeZ();
146 
147 
148  do
149  {
150  // boucle echelle
151  std::cout << std::endl << std::endl << "PYRAMIDE NIVEAU : " << scaleControl.getScale() << std::endl;
152 
153  transformation.setdeltaprev(delta_init);
154  displacementField.init( ref, scaleControl, _seuils);
155 
156  count = 1; // count compte le nbre d'iterations sans changer d'echelle
157 
158 
159  do
160  {
161  // boucle delta
162  std::cout<<" ITERATION No "<<count<<std::endl;
163 
164  // mise a jour transfo ancienne a ancienne nouvelle...
165  p = q * p ;
166 
167  // Calcul effectif du champ
168  AimsData< Point3d > df = displacementField.getField( testtrans);
169 
170  // Attention, le champ va de Ref a Test et non l'inverse !!!!!!!!!!!!!!!!!!!!!!!!
171 
172 
173  // Calcul de la transformation d'etape a partir du champ estime, rappel : q est un motion
174  q = minimisation.quaternion(displacementField); //transfo de test vers ref !!!!
175 
176 
177  // Rappel : r est la transformation totale actualisee, q la calculee courante et p l'ancienne totale
178  // On actualise effectivement r en composant q avec r
179  r = q * p ;
180 
181  // MISE DES TRANSLATIONS EN MILLIMETRES !
182  r.translation()[0] *= test.sizeX();
183  r.translation()[1] *= test.sizeY();
184  r.translation()[2] *= test.sizeZ();
185 
186 
187  // Resampling de test en testtrans par r la transfo totale la plus recente
188  LinearResampler<T> resampler;
189  resampler.setRef( test );
190  testtrans = resampler.doit( r, test.dimX(), test.dimY(), test.dimZ(),
191  Point3df(test.sizeX(),test.sizeY(),test.sizeZ()));
192 
193 
194  // Compteur du nombre d'iterations a cette echelle
195  ++count ;
196 
197  } while ( count < (_itermax + 1) ); // On fait itermax iterations exactement a chaque niveau de pyramide
198 
199 
200  // On a fini d'iterer a cette echelle, on monte un niveau de pyramide
201  scaleControl.nextScale();
202 
203 
204  } while (scaleControl.goOn()); // test fin pyramide
205 
206 
207  _result = test_orig.clone();
208  LinearResampler<T> linres;
209  linres.setRef( test_orig );
210  _result = linres.doit( r,
211  test_orig.dimX(),
212  test_orig.dimY(),
213  test_orig.dimZ(),
214  Point3df(test_orig.sizeX(),test_orig.sizeY(),test_orig.sizeZ()));
215 
216  _done=true;
217 
218 // std::cout<<"C'est fini!"<<std::endl<<"Le motion trouve est : "<<std::endl;
219 // std::cout << r << std::endl;
220  return(r);
221 }
222 
223 #endif
224 
void setZ(int Z)
AimsData< Point3d > getField(AimsData< T > &test)
void setX(int X)
int dimX() const
float sizeY() const
void init(const AimsData< T > &ref, int level_start, int level_stop, double cutVar, double stopVar, double seuilCorrel, const Point3d &tailleBloc=Point3d(4, 4, 4))
void setcx(float cx)
Motion quaternion(DisplacementField< T > &displacementField)
void setRef(const AimsData< T > &ref)
set the input data to be resampled
Definition: resampler.h:125
void setiterMax(int max)
void setdeltaprev(float d)
void setcz(float cz)
void setY(int Y)
double getcutVar()
Definition: scale_control.h:34
void nextScale()
static _Tp min()
AimsData< T > clone() const
virtual void doit(const aims::AffineTransformation3d &transform, AimsData< T > &output_data)
Resample the reference input data (set via setRef()) into an existing output data.
Definition: resampler_d.h:169
float sizeZ() const
double getstopVar()
Definition: scale_control.h:35
void setcy(float cy)
float sizeX() const
int dimZ() const
int dimY() const
int getlevel_stop()
Definition: scale_control.h:31
void init(AimsData< T > &ref, ScaleControl &scaleControl, T *seuils)
aims::AffineTransformation3d doit(AimsData< T > &ref, const AimsData< T > &test_orig)
static _Tp max()