A.I.M.S algorithms


displacement_field_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_DISPLACEMENT_FIELD_D_H
13 #define AIMS_REGISTRATION_DISPLACEMENT_FIELD_D_H
14 
15 
17 #include <aims/data/data_g.h>
18 #include <aims/io/reader.h>
19 #include <aims/io/writer.h>
20 
21 
22 
24 #include <cartobase/type/limits.h>
25 #include <iostream>
27 #include <aims/io/reader.h>
28 #include <aims/io/writer.h>
30 
31 #ifdef AIMS_USE_MPI
32 #include <mpi.h>
33 
34 using namespace MPI;
35 #endif
36 
37 
38 // fonction calculant la moyenne et variance de l'imagette(i,j) en parcourant l'offset
39 template<class T> inline
40 static void meanVar(const AimsData<T>& im, const AimsData<int>& offset, int i, int j, int k, double &m, double &s, const T &lth, const T &hth)
41 {
43  AimsData<int>::const_iterator eoit = offset.end(); //important optimisation not performed by g++ -O2
44  // pointeur first sur le niveau im(i,j) du premier point de l'imagette de coin sup. gauche (i,j)
45  T *first = (T*) (im.begin() + im.oFirstPoint() + j*im.oLine() + i + k*im.oSlice() );
46  T * pt ; //pointeur courant parcourant l'imagette
47  int N = offset.dimX() * offset.dimY() * offset.dimZ(); // nombre total de points dans l'imagette
48 
49  m = 0. ; s = 0. ; //initialisation de moyenne et variance a 0
50 
51  for(oit =offset.begin(); oit != eoit; oit++ )
52  {
53  pt = first + *oit; // absolu + relatif
54 
55  // Ca bogue ici parfois (sortie de l'image)
56  if( (*pt<lth) || (*pt>hth) ) {s=0; m=0; break;}
57 
58  m += *pt;
59  s += *pt * *pt;
60  }
61 
62  m /= N;
63 
64  s /= N; // E(X^2) ...
65  s -= (m * m); // ... - (E(X))^2
66 
67  // On decide de renvoyer la variance qui est le carre de l'ecart type
68  // s = sqrt( s ); // passage de variance a ecart type en prenant la racine carree
69 
70 
71 }
72 
73 template<class T>
74 void DisplacementField<T>::selectBlock( const AimsData<T>& image, double cV)
75 {
76  double m = 0.,s = 0.;
77 
78  int i, j, k;
79 
80  // Ecriture sur disque de pyrRef
81  //aims::Writer< AimsData<T> > ecr_T("pyrRef" + carto::toString(_level) );
82  //ecr_T.write(image);
83 
84  _sel.clear();
85  _moy.clear();
86 
87 
88  ForEach3d(image, i, j, k)
89  {
90  //std::cout << "selectBlock - [" << i << ", " << j << ", " << k << " ], [" << (image.dimX()- _nx) << ", " << (image.dimY() - _ny) << ", " << (image.dimZ() - _nz23d) << "]"<< std::endl << std::flush;
91  if ( ( i < image.dimX() - _nx) && (j < image.dimY() - _ny) && (k < image.dimZ() - _nz23d) )
92  {
93  meanVar( image, _offset, i, j, k, m, s, _lthr, _hthr);
94  //std::cout << "selectBlock - [" << i << ", " << j << ", " << k << " ], mean : " << m << ", var : " << s << std::endl << std::flush;
95  if( s > 0. )
96  {
97  _sel.insert(std::pair<float, Point3d>(s, Point3d(i,j, k) ) ); // on remplit 2 map en classant par variance croissante
98  _moy.insert(std::pair<float, double>(s, m ) ); // l'une avec les coord, l'autre avec la moyenne correspondante
99  }
100  }
101  }
102 
103  _selThresh = _sel.size() * cV;
104 
105 
106 }
107 
108 
109 template<class T>
110 void DisplacementField<T>::init(AimsData<T>& ref, ScaleControl& scaleControl, T* seuils)
111 {
112  // Initialisation des caracteristiques de ref
113  _dimRef = Point3d(ref.dimX(), ref.dimY(), ref.dimZ());
114  _sizeRef = Point3df(ref.sizeX(), ref.sizeY(), ref.sizeZ());
115 
116  // Initialisation de l'echelle
117  _level = scaleControl.getScale();
118 
119  // taille des blocs
120  _nx=(scaleControl.getTailleBloc()[0]);
121  _ny=(scaleControl.getTailleBloc()[1]);
122  _nz=(scaleControl.getTailleBloc()[2]);
123  _nz23d=( _nz == 1 ? 0 : _nz );
124 
125 
126  // On fixe le seuil sur le coeff de simi
127  _seuilCorrel = scaleControl.getSeuilCorrel();
128 
129  // On rentre les seuils sur les niveaux de gris
130  _lthr = seuils[0];
131  _hthr = seuils[1];
132  _ltht = seuils[2];
133  _htht = seuils[3];
134 
135 
136  // Construction de la pyramide pour l'image de reference
137  Pyramid<T> *Pref = NULL ;
138  PyramidFunc<T> *pyramidfunc = NULL;
139 
140  pyramidfunc = new MedianPyramidFunc<T>();
141  Pref = new Pyramid<T>(*pyramidfunc);
142  Pref->setRef( ref );
143  Pref->setLevel( _level );
144 
145  _pyrRef = Pref->item( _level ).clone(); // _pyrRef est l'image au niveau _level
146 // aims::Writer<AimsData<T> > w("pyref_debug.ima");
147 // w.write(_pyrRef);
148 
149  // menage...
150  delete Pref ;
151  delete pyramidfunc ;
152 
153 
154  // ATTENTION PROBLEME PYRAMIDE, SOLUTION TEMPORAIRE!!!!!!!!!!!!!!!!!!!
155  // _pyrRef = ref;
156 
157 
158  // Creation des offsets comme pre-calcul des coordonnees
159 // std::cerr<<"offsets"<<std::endl;
160  AimsData<int> tmp(_nx,_ny,_nz);
161  _offset = tmp;
162  int i, j, k;
163  ForEach3d(_offset, i, j, k)
164  {
165  _offset(i, j, k) = i + j * _pyrRef.oLine() + k * _pyrRef.oSlice() ;
166  }
167 
168  // Selection des blocks
169  //std::cout << "cur var : " << scaleControl.getcutVar() << std::endl << std::flush;
170  selectBlock( _pyrRef, scaleControl.getcutVar());
171 }
172 
173 
174 
175 template <class T>
177 {
178 
179  // Initialisation des barycentres
180  _baryref[0] = _baryref[1] = _baryref[2] = 0.;
181  _barytest[0] = _barytest[1] =_barytest[2] = 0.;
182 
183  // Initialisation de moyenne et variance courante pour calcul correlation
184  double mLocalTest = 0, sLocalTest = 0;
185 
186 
188 
189  // Declaration pyramide pour l'image testtrans (transformee courante)
190  Pyramid<T> *Ptest = NULL ;
191  PyramidFunc<T> *pyramidfunc = NULL;
192 
193  pyramidfunc = new MedianPyramidFunc<T>();
194  Ptest = new Pyramid<T>(*pyramidfunc);
195  Ptest->setRef(test);
196  Ptest->setLevel( _level );
197 
198  _pyrTest = Ptest->item( _level ).clone(); // _pyrTest est l'image testtrans degrade niveau _level
199 
200  // menage...
201  delete Ptest ;
202  delete pyramidfunc ;
203 
204 
205  // Creation des offsets comme pre-calcul des coordonnees pour Test si dim differentes
206  AimsData<int> tmp(_nx,_ny,_nz);
207  _offset2 = tmp;
208  int i, j, k;
209  ForEach3d(_offset2, i, j, k)
210  {
211  _offset2(i, j, k) = i + j * _pyrTest.oLine() + k * _pyrTest.oSlice() ;
212  }
213 
214 
215 
216 
218 
219 
220  // Declaration de l'image du champ niveau _level, qui sera renvoyee en sortie
221  AimsData<Point3d> field( _pyrRef.dimX(), _pyrRef.dimY(), _pyrRef.dimZ());
222  field.setSizeXYZT( _pyrRef.sizeX(), _pyrRef.sizeY(), _pyrRef.sizeZ(), _pyrRef.sizeT() );
223  field = Point3d(0, 0, 0) ;
224 
225 
226  std::multimap<float, Point3d>::reverse_iterator it = _sel.rbegin();
227  std::multimap<float, double>::reverse_iterator itm = _moy.rbegin();
228  AimsData<int>::iterator oit, oit2;
229  std::vector<Point3df> pointstest, pointsref; /* PAR: useless temporary vectors? */
230 
231  // Calcul de n puissance 4, pour prendre ro carre, pour positivite entre 0 et 1
232  double norm_fctr = _offset.dimX() * _offset.dimY() * _offset.dimZ() * _offset.dimX() * _offset.dimY() * _offset.dimZ();
233 // norm_fctr *= norm_fctr; // mise au carre
234 
235  float puis_level = float(pow((double) 2, (double) _level));
236 // std::cout << "_level = " << _level << "\tpuiss_level = " << puis_level
237 // << "\t _nx = " << _nx << "\t _ny = " << _ny <<std::endl ;
238  int realPointNumber = 0 ;
239  int entLimit = int(_selThresh + .5); //after this number of blocks, all significant (in termes of enthropy) blocks have been dealt with
240  // Declaration d'une image des blocs selectionnes par variance uniquement sur Ref
241 
242  /*AimsData<double> selected_blocks(_pyrRef.dimX(),
243  _pyrRef.dimY(),
244  _pyrRef.dimZ() );
245  selected_blocks.setSizeXYZT(_pyrRef.sizeX(),
246  _pyrRef.sizeY(),
247  _pyrRef.sizeZ(), 1.0);
248  selected_blocks=0;
249 
250 
251  AimsData<double> test_match(_pyrTest.dimX(),
252  _pyrTest.dimY(),
253  _pyrTest.dimZ() );
254  test_match.setSizeXYZT(_pyrTest.sizeX(),
255  _pyrTest.sizeY(),
256  _pyrTest.sizeZ(), 1.0);
257  test_match=0;*/
258 
259 #ifdef AIMS_USE_MPI
260  //Calculate index info. Homogenous mapping of blocks on procs, "first" procs recieve one more block than last.
261  int rank = COMM_WORLD.Get_rank(); //we recalculate rank and size for each getField (~9-15 times), so as to not modify the prototype
262  int size = COMM_WORLD.Get_size();
263  int nbBlock = entLimit/size; //entLimit may change so it's ok.
264  int nbBlockLeft = entLimit%size;
265  if(rank < nbBlockLeft) {
266  for(j = 0;j < rank * (nbBlock + 1);j++) { //can't iterate non-sequentially on reverse_iterators, so should change _sel to descending multimap and use a normal iterator..
267  it++;
268  itm++;
269  }
270  entLimit = nbBlock + 1;
271  }
272  else {
273  for(j = 0;j < nbBlockLeft + rank * nbBlock;j++) {
274  it++;
275  itm++;
276  }
277  entLimit = nbBlock;
278  }
279  //Buffer declarations
280  float* pointsRecv; //for recieving pointrefs and pointtests
281  int recvCounts[size]; //as in MPI
282  int dspls[size]; //as in MPI
283 #endif
284  /* PAR: Important optimisation that g++ does not perform */
285  AimsData<int>::const_iterator oitend = _offset.end();
286 
287  // boucle sur tous les blocs selectionnes dans pyrRef
288  for(int n = 0; n < entLimit; ++n, ++it, ++itm)
289  {
290  Point3d p = it->second; // p contient les coordonnees des points de ref selectionnes
291  double mo = itm->second; // moyenne du bloc de coin sup. gauche p
292  double sd = itm->first; // variance du meme bloc
293  int k, l, m; // parcours voisinage bloc...
294  double n4_by_sd = norm_fctr * sd; // pre-calcul facteur normalisation
295  bool badPoint = true; // booleen pour savoir si le correl est correct
296  bool one_max = true; // booleen pour assurer unicite correl max
297  double correlMax = 0; // on a normalement 0 <= correl_au_carre <= 1
298 
299 
300  // boucle recherche bloc correspondant dans voisinage
301  for(k=-3;k<4;++k)
302  {
303  for(l=-3;l<4;++l)
304  {
305  for(m=-3;m<4;++m)
306  {
307  if( (p[0]+k >= 0) &&
308  (p[1]+l >= 0) &&
309  (p[2]+m >= 0) &&
310  (p[0]+k+_nx < _pyrTest.dimX() ) &&
311  (p[1]+l+_ny < _pyrTest.dimY() ) &&
312  (p[2]+m+ _nz23d < _pyrTest.dimZ() ) )
313  {
314  T *pRef = _pyrRef.begin() + _pyrRef.oFirstPoint() + p[0] + p[1]*_pyrRef.oLine() + p[2]*_pyrRef.oSlice() ;
315  T *pTest = _pyrTest.begin() + _pyrTest.oFirstPoint() + (p[0]+k) + (p[1]+l)*_pyrTest.oLine() + (p[2]+m)*_pyrTest.oSlice() ;
316 
317  meanVar( _pyrTest, _offset2, p[0]+k, p[1]+l, p[2]+m, mLocalTest, sLocalTest, _ltht, _htht);
318  //std::cout << "mLocalTest : " << mLocalTest << ", sLocalTest : " << sLocalTest << std::endl << std::flush;
319 
320  //test_match(p[0]+k + _nx/2, p[1]+l + _ny/2, p[2]+m + _nz/2) = mLocalTest;
321 
322 
323  // calcul coefficient de correlation entre les deux blocs courants
324  double correl = 0. ;
325  for(oit = _offset.begin(), oit2 = _offset2.begin(); oit != oitend; oit++, oit2++ )
326  {
327  correl += ( *(pRef + *oit) - mo ) * ( *(pTest + *oit2 ) - mLocalTest);
328  }
329  correl *= correl; // mise au carre
330 
331  if(sLocalTest)
332  {
333  correl /= (n4_by_sd * sLocalTest); // normalisation
334  }
335  else correl = 0;
336 
337  if(correl) badPoint = false;
338 
339  // selection du correlMax
340  if ( (correl == correlMax ) && (!badPoint) )
341  {
342  one_max = false;
343  field( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2) = Point3d(0,0,0);
344  }
345  else if ( (correl > correlMax ) && (!badPoint) )
346  {
347  one_max = true;
348  correlMax = correl;
349  field( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2) = Point3d(k*(int16_t)puis_level, l*(int16_t)puis_level, m*(int16_t)puis_level);
350  }
351  }
352  }
353  }
354  }
355 
356  // Calcul des coordonnes REELLES des CENTRES des blocs
357  Point3df curRef, curTest;
358 
359  curRef[0] = ( (float)p[0] + float(_nx)/2.0 ) * puis_level;
360  curRef[1] = ( (float)p[1] + float(_ny)/2.0 ) * puis_level;
361  curRef[2] = ( (float)p[2] + float(_nz23d)/2.0 ) * puis_level;
362 
363  curTest[0] = curRef[0] + float( field( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2)[0] );
364  curTest[1] = curRef[1] + float( field( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2)[1] );
365  curTest[2] = curRef[2] + float( field( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2)[2] );
366 
367 /* curTest[0] = curRef[0] + float( (field(p))[0] ) * puis_level;
368  curTest[1] = curRef[1] + float( (field(p))[1] ) * puis_level;
369  curTest[2] = curRef[2] + float( (field(p))[2] ) * puis_level;*/
370 
371 
372  // Ecriture de l'image des blocs selectionnes uniquement par variance sur Ref
373  //selected_blocks( unsigned( (curRef[0])/10.0 + .5 ), unsigned( (curRef[1])/10.0 + .5) )= correlMax;
374 
375  //selected_blocks( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2) = correlMax;
376 
377  if ( (!badPoint) && ( correlMax > _seuilCorrel ) && (one_max) )
378  {
379  // On ajoute le centre REEL des 2 blocs apparies dans 2 listes
380  pointsref.push_back( curRef );
381  pointstest.push_back( curTest );
382 
383  // Calcul du barycentre REEL des 2 listes de points
384  _baryref[0] += curRef[0]; _baryref[1] += curRef[1]; _baryref[2] += curRef[2];
385  _barytest[0] += curTest[0]; _barytest[1] += curTest[1]; _barytest[2] += curTest[2];
386  ++realPointNumber;
387  }
388  else
389  {
390  //selected_blocks( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2 ) = 255;
391  field( p[0] + _nx/2, p[1] + _ny/2, p[2] + _nz23d/2) = Point3d(0,0,0);
392  }
393  } //fin boucle sur chaque bloc
394 
395  _pointstest = pointstest;
396  _pointsref = pointsref;
397 #ifdef AIMS_USE_MPI
398  /* Mutualize calculated fieldpoints.. (order of points is important, see yx in AffineLeastSquareEstimation::computeMotion()) */
399  //Tech:Can't place directly in pointRecv instead of pushing back cause don't know where. So copy points into oversized buffer, than to pointRecv, or
400  //get from _pointref, almost as fast..
401  recvCounts[rank] = 6*realPointNumber;
402  COMM_WORLD.Allgather(IN_PLACE, 0, INT, recvCounts, 1, INT);
403  COMM_WORLD.Barrier();
404 
405  dspls[0] = 0;
406  for(i = 1; i < size; i++) { //faster than MPI allgather
407  dspls[i] = dspls[i-1] + recvCounts[i-1];
408  }
409 
410  int totalNumFloats = dspls[size-1] + recvCounts[size-1]; //total number of points across all processes * 6
411  pointsRecv = new float[totalNumFloats];
412  for(i = 0; i < realPointNumber; i++) {
413  for(int j = 0; j < 3; j++) {
414  pointsRecv[dspls[rank] + i*6 + j] = (_pointsref[i])[j];
415  pointsRecv[dspls[rank] + i*6 + 3 + j] = (_pointstest[i])[j];
416  }
417  }
418  _pointsref.clear();_pointstest.clear();
419  COMM_WORLD.Allgatherv(IN_PLACE, 0, INT, pointsRecv, recvCounts, dspls, FLOAT);
420  COMM_WORLD.Barrier();
421  for(int i = 0; i < totalNumFloats; i += 6) {
422  _pointsref.push_back(Point3df(pointsRecv[i],pointsRecv[i+1],pointsRecv[i+2]));
423  _pointstest.push_back(Point3df(pointsRecv[i+3],pointsRecv[i+4],pointsRecv[i+5]));
424  }
425  delete[] pointsRecv;
426 
427  /*..and centers of mass */
428  realPointNumber = totalNumFloats/6;
429  COMM_WORLD.Allreduce(IN_PLACE, _baryref, 3, DOUBLE, SUM);
430  COMM_WORLD.Allreduce(IN_PLACE, _barytest, 3, DOUBLE, SUM);
431 #endif
432 
433 // Writer< AimsData<T> > testWriter("pyrTest_" + toString(_level) );
434 // testWriter.write(_pyrTest);
435 
436 
437  // Ecriture sur disque de selected_bloks
438 // aims::Writer< AimsData<double> > ecr_T("selected_blocks" + carto::toString(_level) );
439 // ecr_T.write(selected_blocks);
440 
441 // aims::Writer< AimsData<double> > ecr_T2("test_match" + carto::toString(_level) );
442 // ecr_T2.write(test_match);
443 
444  _baryref[0] /= realPointNumber;
445  _baryref[1] /= realPointNumber;
446  _baryref[2] /= realPointNumber;
447  _barytest[0] /= realPointNumber;
448  _barytest[1] /= realPointNumber;
449  _barytest[2] /= realPointNumber;
450 
451  return ( field );
452 }
453 
454 
455 template <class T>
457 {
458  // Declaration pyramide pour l'image testtrans (transformee finale au niveau _level)
459  Pyramid<T> *Ptesttrans = NULL ;
460  PyramidFunc<T> *pyramidfunc = NULL;
461 
462  pyramidfunc = new MedianPyramidFunc<T>();
463  Ptesttrans = new Pyramid<T>(*pyramidfunc);
464  Ptesttrans->setRef(testtrans);
465  Ptesttrans->setLevel( level );
466 
467  AimsData<T> pyrTrans = Ptesttrans->item( level ).clone(); // pyrTrans est l'image testtrans degradee niveau level
468 
469  // menage...
470  delete Ptesttrans ;
471  delete pyramidfunc ;
472 
473  // Declaration pyramide pour l'image ref (transformee en niveau _level)
474  Pyramid<T> *Pref = NULL ;
475  PyramidFunc<T> *pyramidfunc2 = NULL;
476 
477  pyramidfunc2 = new MedianPyramidFunc<T>();
478  Pref = new Pyramid<T>(*pyramidfunc2);
479  Pref->setRef(ref);
480  Pref->setLevel( level );
481 
482  AimsData<T> pyrRef = Pref->item( level ).clone(); // pyrRef est l'image ref degradee niveau level
483 
484  // menage...
485  delete Pref ;
486  delete pyramidfunc2 ;
487 
488 
489 
490  // Creation des offsets comme pre-calcul des coordonnees
491  AimsData<int> tmp(_nx,_ny,_nz);
492  _offset = tmp.clone();
493  _offset2 = tmp.clone();
494  int i, j, k;
495  ForEach3d(_offset2, i, j, k)
496  {
497  _offset(i, j, k) = i + j * pyrRef.oLine() + k * pyrRef.oSlice() ;
498  _offset2(i, j, k) = i + j * pyrTrans.oLine() + k * pyrTrans.oSlice() ;
499  }
500 
501 
502 
503 
504 
505  // Initialisation de moyenne et variance courante pour calcul correlation
506  double mRef = 0, sRef = 0, mTest = 0, sTest = 0;
507 
508  // Creation d'une image affichant la correlation en chaque point e un niveau de pyramide
509  AimsData<T> pyrQuality = pyrTrans.clone();
510 
512  pyrQuality = 0;
513 
514  AimsData<int>::iterator oit, oit2;
515  AimsData<int>::const_iterator eoit = _offset.end(); //important optimisation that g++ does not perform
516 
517  // Calcul de n puissance 4, pour prendre ro carre, pour positivite entre 0 et 1
518  double norm = _offset.dimX() * _offset.dimY() * _offset.dimZ() * _offset.dimX() * _offset.dimY() * _offset.dimZ();
519 
521 // std::cout<<"thresh : "<<thresh<<std::endl;
522 // std::cout<<"_nz23d : "<<_nz23d<<std::endl;
523 
524 // aims::Writer< AimsData<T> > pyrTransWriter("pyrTrans");
525 // pyrTransWriter.write(pyrTrans);
526 
527 
528 
529  ForEach3d(pyrTrans, i, j, k)
530  {
531 /*
532  std::cout<<"i, j, k : "<<i<<", "<<j<<", "<<k<<std::endl;
533  std::cout<<"i + _nx, j + _ny, k + _nz : "<<i + _nx<<", "<<j + _ny<<", "<<k + _nz23d*_nz<<std::endl;
534  std::cout<<"pyrTrans.dimX(), pyrTrans.dimY(), pyrTrans.dimZ() : "<<pyrTrans.dimX()<<", "<<pyrTrans.dimY()<<", "<<pyrTrans.dimZ()<<std::endl;
535 */
536 
537 
538  if( (i + _nx < pyrTrans.dimX() ) && (j + _ny < pyrTrans.dimY() ) && (k + _nz23d*_nz < pyrTrans.dimZ() ) )
539  {
541  // std::cout<<"boucle 1"<<std::endl;
542  if( (pyrTrans(i + _nx/2,j + _ny/2, k + _nz23d*_nz/2) < thresh) /*&& (pyrRef(i,j) < thresh) */)
543  {
544  T *pRef = pyrRef.begin() + pyrRef.oFirstPoint() + i + j*pyrRef.oLine() + k*pyrRef.oSlice();
545  T *pTest = pyrTrans.begin() + pyrTrans.oFirstPoint() + i + j*pyrTrans.oLine() + k*pyrRef.oSlice();
546 
547 
548  meanVar( pyrTrans, _offset2, i, j , k, mTest, sTest, _ltht, _htht);
549  meanVar( pyrRef, _offset, i, j , k, mRef, sRef, _lthr, _hthr);
550  double normS = norm * (sRef*sTest);
552 // std::cout<<"boucle 2"<<std::endl;
553 
554 
555  if(normS!=0)
556  {
557  // calcul coefficient de correlation entre les deux blocs courants
558  double correl = 0. ;
559  for(oit = _offset.begin(), oit2 = _offset2.begin(); oit != eoit; oit++, oit2++ )
560  {
561  correl += ( *(pRef + *oit) - mRef ) * ( *(pTest + *oit2 ) - mTest);
562  }
563  correl *= correl; // mise au carre
564  correl /= normS; // normalisation
565  correl = sqrt(correl);
566 
568 // std::cout<<"correl : "<<correl<<std::endl;
569 // std::cout<<"correl rescale : "<<T((std::numeric_limits<T>::max() -1)*correl + 0.5)<<std::endl;
570 
571  pyrQuality(i + _nx/2,j + _ny/2, k + _nz23d*_nz/2) = T((std::numeric_limits<T>::max() -1)*correl + 0.5);
572  }
573 
575 // else pyrQuality(i + _nx/2,j + _ny/2, k + _nz23d*_nz/2) = 1000;
576  }
577  else pyrQuality(i + _nx/2,j + _ny/2, k + _nz23d*_nz/2) = 0;
578  }
579  }
580 
581 
582  return(pyrQuality);
583 }
584 
585 
586 #endif
587 
AimsData< Point3d > getField(AimsData< T > &test)
int dimX() const
float sizeY() const
iterator begin()
#define ForEach3d(thing, x, y, z)
Point3d getTailleBloc()
Definition: scale_control.h:37
AimsData< T > getQuality(AimsData< T > &testtrans, AimsData< T > &ref, int level, T thresh=std::numeric_limits< T >::max())
double getcutVar()
Definition: scale_control.h:34
AimsData< T > clone() const
double getSeuilCorrel()
Definition: scale_control.h:36
#define FLOAT
float sizeZ() const
int oLine() const
float sizeX() const
int dimZ() const
void setSizeXYZT(float sizex=1.0f, float sizey=1.0f, float sizez=1.0f, float sizet=1.0f)
int dimY() const
AIMSDATA_API float norm(const Tensor &thing)
iterator begin()
iterator end()
int oSlice() const
#define DOUBLE
void init(AimsData< T > &ref, ScaleControl &scaleControl, T *seuils)
int oFirstPoint() const