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