aimsalgo 6.0.0
Neuroimaging image processing
primalSketch.h
Go to the documentation of this file.
1/* This software and supporting documentation are distributed by
2 * Institut Federatif de Recherche 49
3 * CEA/NeuroSpin, Batiment 145,
4 * 91191 Gif-sur-Yvette cedex
5 * France
6 *
7 * This software is governed by the CeCILL-B license under
8 * French law and abiding by the rules of distribution of free software.
9 * You can use, modify and/or redistribute the software under the
10 * terms of the CeCILL-B license as circulated by CEA, CNRS
11 * and INRIA at the following URL "http://www.cecill.info".
12 *
13 * As a counterpart to the access to the source code and rights to copy,
14 * modify and redistribute granted by the license, users are provided only
15 * with a limited warranty and the software's author, the holder of the
16 * economic rights, and the successive licensors have only limited
17 * liability.
18 *
19 * In this respect, the user's attention is drawn to the risks associated
20 * with loading, using, modifying and/or developing or reproducing the
21 * software by the user in light of its specific status of free software,
22 * that may mean that it is complicated to manipulate, and that also
23 * therefore means that it is reserved for developers and experienced
24 * professionals having in-depth computer knowledge. Users are therefore
25 * encouraged to load and test the software's suitability as regards their
26 * requirements in conditions enabling the security of their systems and/or
27 * data to be ensured and, more generally, to use and operate it in the
28 * same conditions as regards security.
29 *
30 * The fact that you are presently reading this means that you have had
31 * knowledge of the CeCILL-B license and that you accept its terms.
32 */
33
34#ifndef AIMS_PRIMALSKETCH_PRIMALSKETCH_H
35#define AIMS_PRIMALSKETCH_PRIMALSKETCH_H
36
37#include <cstdlib>
42#include <set>
43#include <map>
44#include <list>
45#include <aims/mesh/surfaceOperation.h>
46#include <aims/mesh/surfacegen.h>
47
48
49namespace aims
50{
51
57
58 template<typename Geom, typename Text> class PrimalSketch
59 {
60
61 protected:
62 typedef typename SiteType<Geom>::type Site;
63 typedef typename TexType<Text>::type Val;
64
65 std::string _subject;
66
68 std::list<Bifurcation<Site>*> bifurcationList;
69 std::list<ScaleSpaceBlob<Site>*> blobList;
70 int labelMax; // ssblob label max (=final numbre of blobs); incremented during detection
71
72 float t_min; // top and bottom scales
73 float t_max;
74
75 TexturedData<Geom, Text> *_mask; // masking : outside the mask, blob are removed
76
77 float min_delta_t; // minimal scale interval for refinements
78
79 int _type; // is it a surface or image-based primal sketch ?
80 // default is image
81
82 public:
83
84 PrimalSketch() : _scaleSpace(NULL), _mask(NULL) {}
85 PrimalSketch( const std::string & subject, int type )
86 : _subject(subject) , _scaleSpace(NULL), labelMax(0),_mask(NULL), _type(type) {}
87 PrimalSketch( const std::string & subject,
89 : _subject(subject) , _scaleSpace(scaleSpace), labelMax(0), _mask(NULL), _type(type)
90 { SetMinDt(); }
91 PrimalSketch( const std::string & subject,
95 { SetMinDt(); }
96
97 std::string Subject() { return _subject; }
98 int Type() { return _type; }
99 void setType( int t ) { _type = t; }
100
101 void SetMinDt()
102 {
103 if (!_scaleSpace->smoother()->optimal())
104 min_delta_t=0.5;
105 else
106 min_delta_t=_scaleSpace->smoother()->dt();
107 }
110 std::list<Bifurcation<Site>*> BifurcationList()
111 { return bifurcationList; }
112 std::list<ScaleSpaceBlob<Site>*> BlobSet() { return blobList; }
114
116 { blobList.push_back(blob); }
118 { bifurcationList.push_back(bifurcation); }
119
120 ScaleSpaceBlob<Site> *Blob( int label );
121
122 void ComputePrimalSketch( float tmin, float tMax,
123 const std::string & statFile="",
124 uint intersection_param=10 );
125
126 void MatchScaleLevels( float t_up, float t_down,
127 uint intersection_param );
128
129 void ComputeBlobMeasurements( const std::string & statName="" );
130
132 GreyLevelBlob<Site> *glBlob );
133
134 };
135
136
137 //--------------------------------------------------------------
138 // IMPLEMENTATION
139 //--------------------------------------------------------------
140
141 template<typename Geom, typename Text>
142 ScaleSpaceBlob<typename SiteType<Geom>::type>
144 {
145 typename std::list<ScaleSpaceBlob<Site>*>::iterator itBlobList = blobList.begin();
146 for (; itBlobList!=blobList.end(); ++itBlobList)
147 {
148 if ((*itBlobList).Label() == label)
149 return (*itBlobList);
150 }
151 }
152
153 //---------------------------------------------------------------------
154
155 template<typename Geom, typename Text> // Pas super optimal mais bon...
157 {
158 typename std::list<ScaleSpaceBlob<Site>*>::iterator itBlobList=blobList.begin();
159 ScaleSpaceBlob<Site> *ssBlob;
160 std::list<GreyLevelBlob<Site>*> glBlobList;
161 typename std::list<GreyLevelBlob<Site>*>::iterator
162 itGLBlobList;
163 for (; itBlobList!=blobList.end(); ++itBlobList)
164 {
165 ssBlob=*itBlobList;
166 glBlobList=ssBlob->glBlobs;
167 itGLBlobList=glBlobList.begin();
168 for (; itGLBlobList!=glBlobList.end(); ++itGLBlobList)
169 if ( ((*itGLBlobList)->Label()==glBlob->Label()) && ((*itGLBlobList)->GetScale()==glBlob->GetScale()) )
170 return ssBlob;
171 }
172 // blob not found, return NULL
173 std::cout << "SS BLOB label=" << glBlob->Label() << " and scale " << glBlob->GetScale() << " cannot be found...." << std::endl;
174 return NULL;
175 }
176
177 //---------------------------------------------------------------------
178
179 template<typename Geom, typename Text>
181 float tmin, float tmax, const std::string & statFile,
182 uint intersection_param )
183 {
184 // Algo : On part del'echelle la plus haute et on descend vers l'image originale
185 // A chaque niveau : on regarde vers le niveau en dessous et on applique
186 // MatchScaleLevels();
187 // Si ca marche on descend d'un cran et on recommence
188 // Si ca ne marche pas on calcule un niveau intermediaire et on recommence (~recursif)
189
190 std::cout << "Computing primal sketch between tmin=" << tmin << " and tmax=" << tmax << std::endl;
191
192 t_min = tmin;
193 t_max = tmax;
194 float t_bif = tmax*2;
195
196 if ( _scaleSpace == NULL ) {
197 std::cerr << "Trying to compute a primal sketch without a scale-space" << std::endl;
198 std::cerr << "Subject : " << _subject;
199 exit(EXIT_FAILURE);
200 }
201
202 std::set<float> scaleList=_scaleSpace->GetScaleList();
203 std::set<float>::reverse_iterator itScaleUp=scaleList.rbegin();
204 std::set<float>::reverse_iterator itScaleDown;
205 float t_up, t_down;
206
207 for ( ; (*itScaleUp) != tmax ; ++itScaleUp ) {
208 if ( itScaleUp == scaleList.rend() ) {
209 std::cerr << "Scale max " << tmax << " does not exist" << std::endl;
210 exit(EXIT_FAILURE);
211 }
212 }
213
214 ScaleLevel<Geom,Text> *levelUp;
215 levelUp = _scaleSpace->Scale(*itScaleUp);
216 ScaleLevel<Geom,Text> *levelDown;
217 levelDown=_scaleSpace->Scale(tmin);
218
219 levelUp->DetectBlobs( _mask );
220
221 GreyLevelBlob<Site> *glBlob;
222 ScaleSpaceBlob<Site> *ssBlob;
223 Bifurcation<Site> *bifurc;
224
225 // Scale-space blob initialisation
226 std::map<int, GreyLevelBlob<Site> *>
227 listeGLB=levelUp->BlobList();
228
229 typename std::map<int, GreyLevelBlob<Site> *>::iterator
230 itGL = listeGLB.begin();
231
232 for ( ; itGL != listeGLB.end() ; ++itGL ) {
233 glBlob=(*itGL).second;
235 ssBlob->SetScaleMax(*itScaleUp);
236 ssBlob->SetScaleMin(t_min);
237 ssBlob->AddGreyLevelBlob(glBlob);
238 bifurc = new Bifurcation<Site>(DISAPPEAR, t_bif, t_max);
239 bifurc->AddBottomBlob(ssBlob);
240 ssBlob->SetTopBifurcation(bifurc);
241 AddBlob(ssBlob);
242 AddBifurcation(bifurc);
243 labelMax++;
244 }
245
246 // Primal Sketch Computation, Top To Bottom
247
248 for ( ; (*itScaleUp) > tmin ; ++itScaleUp ) {
249 itScaleDown=itScaleUp;
250 itScaleDown++;
251 t_up = *itScaleUp;
252 t_down = *itScaleDown;
253 MatchScaleLevels(t_up, t_down, intersection_param);
254 }
255
256 // "closing" the primal sketch at bottom
257
258 listeGLB = levelDown->BlobList();
259 itGL = listeGLB.begin();
260 t_bif = t_min / 2.0;
261 for ( ; itGL != listeGLB.end() ; ++itGL ) {
262 glBlob = (*itGL).second;
263 ssBlob = GetSSBlobFromGLBlob(glBlob);
264 if ( ssBlob != NULL) {
265 bifurc = new Bifurcation<Site>(APPEAR, t_min, t_bif);
266 ssBlob->SetBottomBifurcation(bifurc);
267 bifurc->AddTopBlob(ssBlob);
268 AddBifurcation(bifurc);
269 }
270 }
271
272 // blob updates
273 typename std::list<ScaleSpaceBlob<Site>*>::iterator
274 itBlob = blobList.begin();
275 for ( ; itBlob != blobList.end() ; ++itBlob ) {
276 (*itBlob)->ComputeLifeTime();
277 (*itBlob)->ComputeScaleRep();
278 }
279 // blob (multiscale) measurements are computed here.
280
281 ComputeBlobMeasurements(statFile);
282 std::cout << "ssblobs avant : " << blobList.size() << std::endl;
283
284 // ELAGAGE
285 typename std::list<Bifurcation<Site>*>::iterator bifit = bifurcationList.begin();
286 //cout << "bifurcation number:" << bifurcationList.size() << endl;
287 uint merge_count = 0;
288
289 for ( bifit = bifurcationList.begin() ; bifit != bifurcationList.end() ; bifit++ )
290 if ((*bifit)->Type() == MERGE)
291 merge_count ++;
292 //cout << "MERGECOUNT:" << merge_count << endl;
293
294
295 for ( bifit = bifurcationList.begin() ; bifit != bifurcationList.end() ; bifit++ ) {
296 if ( (*bifit)->Type() == APPEAR ) {
297 // std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->TopBlobs();
298 // typename std::list<ScaleSpaceBlob<Site>*>::iterator topit=blist.begin();
299 // ScaleSpaceBlob<Site> *ssblob1;
300 // ssblob1 = *topit;
301 // if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1) {
302 // int label1 = ssblob1->Label();
303 // for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
304 // if ((*itBlob)->Label() == label1) {
305 // blobList.erase(itBlob);
306 // break;
307 // }
308 //
309 // }
310
311 }
312 else if ( (*bifit)->Type() == DISAPPEAR ) {
313 // std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->BottomBlobs();
314 // typename std::list<ScaleSpaceBlob<Site>*>::iterator bottomit=blist.begin();
315 // ScaleSpaceBlob<Site> *ssblob1;
316 // ssblob1 = *bottomit;
317 // if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1) {
318 // int label1 = ssblob1->Label();
319 // for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
320 // if ((*itBlob)->Label() == label1) {
321 // blobList.erase(itBlob);
322 // break;
323 // }
324 // }
325 }
326 else if ( (*bifit)->Type() == MERGE ) {
327 //cout << "merge" << (*bifit)->BottomBlobs().size()<< ";" << (*bifit)->TopBlobs().size() << endl;
328 std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->BottomBlobs();
329 typename std::list<ScaleSpaceBlob<Site>*>::iterator bottomit=blist.begin();
330 ScaleSpaceBlob<Site> *ssblob1, *ssblob2;
331 ssblob1 = *bottomit;
332 bottomit++;
333 ssblob2 = *bottomit;
334 /* float area1 = (*ssblob1).GetMeasurements().area;
335 float area2 = (*ssblob2).GetMeasurements().area; */
336 //cout << /*area1 << ";" << area2 << ";" <<*/ max(area1,area2)/min(area1,area2) << endl;
337 if ((*ssblob2).GlBlobRep()->GetListePoints().size() == 1) {
338 //cout << "ELAG" << endl;
339 // LE BLOB2 EST PETIT : ON FUSIONNE BLOB1 AVEC TOPBLOB ET TOPBIF DE BLOB2 DEVIENT DISPARITION
340
341 ssblob2->TopBifurcation()->setType(DISAPPEAR);
342 ssblob2->TopBifurcation()->TopBlobs().clear();
343 ssblob2->TopBifurcation()->BottomBlobs().clear();
344 ssblob2->TopBifurcation()->AddBottomBlob(ssblob2);
345
346 blist = (*bifit)->TopBlobs();
347 ScaleSpaceBlob<Site> *topblob = *(blist.begin());
348 // RAJOUT DES GLB DE SSBLOB1
349 std::list<GreyLevelBlob<Site>*> glb = ssblob1->glBlobs;
350 typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
351 for (; glbit != glb.end() ; glbit++)
352 topblob->AddGreyLevelBlob(*glbit);
353
354 topblob->SetBottomBifurcation(ssblob1->BottomBifurcation());
355 topblob->SetScaleMin(ssblob1->ScaleMin());
356
357 int label1 = ssblob1->Label();
358 for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
359 if ((*itBlob)->Label() == label1) {
360 blobList.erase(itBlob);
361 break;
362 }
363 }
364 else if ( (*ssblob1).GlBlobRep()->GetListePoints().size() == 1 ) {
365 // LE BLOB1 EST PETIT : ON FUSIONNE BLOB2 AVEC TOPBLOB ET TOPBIF DE BLOB1 DEVIENT DISPARITION
366 //cout << "ELAG" << endl;
367 ssblob1->TopBifurcation()->setType(DISAPPEAR);
368 ssblob1->TopBifurcation()->TopBlobs().clear();
369 ssblob1->TopBifurcation()->BottomBlobs().clear();
370 ssblob1->TopBifurcation()->AddBottomBlob(ssblob1);
371
372 blist = (*bifit)->TopBlobs();
373 ScaleSpaceBlob<Site> *topblob = *(blist.begin());
374 // RAJOUT DES GLB DE SSBLOB2
375 std::list<GreyLevelBlob<Site>*> glb = ssblob2->glBlobs;
376 typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
377 for (; glbit != glb.end() ; glbit++)
378 topblob->AddGreyLevelBlob(*glbit);
379
380 topblob->SetBottomBifurcation(ssblob2->BottomBifurcation());
381 topblob->SetScaleMin(ssblob2->ScaleMin());
382
383 int label2 = ssblob2->Label();
384 for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
385 if ((*itBlob)->Label() == label2) {
386 blobList.erase(itBlob);
387 break;
388 }
389 }
390
391 }
392 else if ((*bifit)->Type() == SPLIT) {
393
394 //cout << "split" << (*bifit)->TopBlobs().size()<< ";" << (*bifit)->BottomBlobs().size() << endl;
395 std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->TopBlobs();
396 typename std::list<ScaleSpaceBlob<Site>*>::iterator topit=blist.begin();
397 ScaleSpaceBlob<Site> *ssblob1, *ssblob2;
398 ssblob1 = *topit;
399 topit++;
400 ssblob2 = *topit;
401 /* float area1 = (*ssblob1).GetMeasurements().area;
402 float area2 = (*ssblob2).GetMeasurements().area; */
403 // cout << /*area1 << ";" << area2 << ";" <<*/ max(area1,area2)/min(area1,area2) << endl;
404 if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1) {
405 // LE BLOB2 EST PETIT : ON FUSIONNE BLOB1 AVEC BOTTOMBLOB ET BOTTOMBIF DE BLOB2 DEVIENT DISPARITION
406
407 ssblob2->BottomBifurcation()->setType(DISAPPEAR);
408 ssblob2->BottomBifurcation()->BottomBlobs().clear();
409 ssblob2->BottomBifurcation()->TopBlobs().clear();
410 ssblob2->BottomBifurcation()->AddTopBlob(ssblob2);
411
412 blist = (*bifit)->BottomBlobs();
413 ScaleSpaceBlob<Site> *topblob = *(blist.begin());
414 // RAJOUT DES GLB DE SSBLOB1
415 std::list<GreyLevelBlob<Site>*> glb = ssblob1->glBlobs;
416 typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
417 for (; glbit != glb.end() ; glbit++)
418 topblob->AddGreyLevelBlob(*glbit);
419
420 topblob->SetTopBifurcation(ssblob1->TopBifurcation());
421 topblob->SetScaleMin(ssblob1->ScaleMin());
422
423 int label1 = ssblob1->Label();
424 for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
425 if ((*itBlob)->Label() == label1) {
426 blobList.erase(itBlob);
427 break;
428 }
429 }
430 else if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1){
431 // LE BLOB1 EST PETIT : ON FUSIONNE BLOB2 AVEC BOTTOMBLOB ET BOTTOMBIF DE BLOB1 DEVIENT DISPARITION
432
433 ssblob1->BottomBifurcation()->setType(DISAPPEAR);
434 ssblob1->BottomBifurcation()->BottomBlobs().clear();
435 ssblob1->BottomBifurcation()->TopBlobs().clear();
436 ssblob1->BottomBifurcation()->AddTopBlob(ssblob1);
437
438 blist = (*bifit)->BottomBlobs();
439 ScaleSpaceBlob<Site> *topblob = * ( blist.begin() );
440 // RAJOUT DES GLB DE SSBLOB2
441 std::list<GreyLevelBlob<Site>*> glb = ssblob2->glBlobs;
442 typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
443 for ( ; glbit != glb.end() ; glbit++ )
444 topblob->AddGreyLevelBlob( *glbit );
445
446 topblob->SetTopBifurcation ( ssblob2->TopBifurcation() );
447 topblob->SetScaleMin ( ssblob2->ScaleMin() );
448
449 int label2 = ssblob2->Label();
450 for ( itBlob = blobList.begin() ; itBlob != blobList.end() ; ++itBlob )
451 if ( (*itBlob)->Label() == label2 ) {
452 blobList.erase ( itBlob );
453 break;
454 }
455 }
456 }
457 }
458 merge_count = 0;
459 bifit = bifurcationList.begin();
460 for (bifit = bifurcationList.begin() ; bifit != bifurcationList.end() ; bifit++)
461 if ((*bifit)->Type() == MERGE)
462 merge_count ++;
463 //cout << "MERGECOUNT:" << merge_count << endl;
464
465 // blob updates
466 std::cout << "ssblobs après : " << blobList.size() << std::endl;
467
468 for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob) {
469 (*itBlob)->ComputeLifeTime();
470 (*itBlob)->ComputeScaleRep();
471 }
472
473 ComputeBlobMeasurements(statFile);
474 }
475
476 //---------------------------------------------------------------------
477
478 template<typename Geom, typename Text>
479 void PrimalSketch<Geom,Text>::MatchScaleLevels(float t_up, float t_down, uint intersection_param)
480 {
481 float t, t_new;
482 int stat1=0, stat2=0, stat3=0;
483
484 std::cout << "Matching " << t_up << " <-> " << t_down << std::endl;
485
486 ScaleLevel<Geom,Text> *levelUp=_scaleSpace->Scale(t_up);
487 ScaleLevel<Geom,Text> *levelDown=_scaleSpace->Scale(t_down);
488 ScaleSpaceBlob<Site> *ssBlobUp, *ssBlobDown;
489 Bifurcation<Site> *bifurc;
490
491 std::cout << "Blob detection at scale " << t_up << std::endl;
492 levelUp->DetectBlobs(_mask);
493 std::cout << "Blob detection at scale " << t_down << std::endl;
494 levelDown->DetectBlobs(_mask);
495
496 std::cout << "\tt=" << t_up << " : " << levelUp->nbBlobs() << " blobs" << std::endl;
497 std::cout << "\tt=" << t_down << " : " << levelDown->nbBlobs() << " blobs" << std::endl;
498
499 std::map<int, GreyLevelBlob<Site> *> listeGLB_up=levelUp->BlobList();
500 std::map<int, GreyLevelBlob<Site> *> listeGLB_down=levelDown->BlobList();
501
502 // On essaye la maniere (simple mais) brutale
503
504 std::cout << "\tComputing potential links" << std::endl;
505
506 std::map<int, std::set<int> > matchUp;
507 std::map<int, std::set<int> > matchDown;
508 typename std::map<int, GreyLevelBlob<Site> *>::iterator
509 iteratorUp;
510 typename std::map<int, GreyLevelBlob<Site> *>::iterator
511 iteratorDown;
512 GreyLevelBlob<Site> *blobUp, *blobDown;
513 int labelUp, labelDown;
514
515 iteratorUp=listeGLB_up.begin();
516 for (; iteratorUp!=listeGLB_up.end(); ++iteratorUp)
517 {
518 blobUp=(*iteratorUp).second; labelUp=(*iteratorUp).first;
519 matchUp[labelUp]=std::set<int>();
520 std::set<Site,ltstr_p3d<Site> > pointsUp=blobUp->GetListePoints();
521 iteratorDown=listeGLB_down.begin();
522 for (; iteratorDown!=listeGLB_down.end(); ++iteratorDown)
523 {
524 blobDown=(*iteratorDown).second; labelDown=(*iteratorDown).first;
525 if (matchDown.find(labelDown) == matchDown.end())
526 matchDown[labelDown]=std::set<int>();
527 // blobs intersect eachother
528 std::set<Site,ltstr_p3d<Site> > pointsInter;
529 std::set<Site,ltstr_p3d<Site> > pointsDown=blobDown->GetListePoints();
530 std::set_intersection(pointsDown.begin(), pointsDown.end(),
531 pointsUp.begin(), pointsUp.end(),
532 std::inserter(pointsInter, pointsInter.begin()),
533 ltstr_p3d<Site>() );
534// if (2*pointsInter.size()/(pointsUp.size() + pointsDown.size()) >0.5)
535 if (pointsInter.size()>intersection_param)
536 {
537 matchUp[labelUp].insert(labelDown);
538 matchDown[labelDown].insert(labelUp);
539 }
540 }
541 }
542
543 //-------------------------------------------------------------------
544 std::cout << "\tChecking potential refinements" << std::endl;
545 std::map<int, std::set<int> >::iterator itMatchUp=matchUp.begin();
546 std::map<int, std::set<int> >::iterator itMatchDown=matchDown.begin();
547 int refinement=0;
548 for (; itMatchUp!=matchUp.end(); ++itMatchUp)
549 {
550 int nUp=(*itMatchUp).second.size();
551 if (nUp > 2)
552 {
553 refinement=2;
554 } //---------------------------------------------------------------------
555 else if (nUp == 2)
556 {
557 std::set<int> setUp=(*itMatchUp).second;
558 std::set<int>::iterator itUp=setUp.begin();
559 for (; itUp!=setUp.end(); ++itUp)
560 if (matchDown[*itUp].size()>1)
561 refinement=1;
562 }
563 if (refinement>0) break;
564 }
565 if (refinement==0)
566 {
567 for (; itMatchDown!=matchDown.end(); ++itMatchDown)
568 {
569 if ((*itMatchDown).second.size() > 2)
570 refinement=3;
571 if (refinement>0) break;
572 }
573 }
574 t=exp((log(t_up) + log(t_down))/2.0); t_new=0;
575 int i=0;
576 while ((i*0.1) <= t) i++;
577 t_new=(i-1)*0.1;
578 if ((refinement>0) && ((t_up-t_new)>min_delta_t) && ((t_new-t_down)>min_delta_t))
579 {
580 std::cout << "\tRefinement needed (type " ;
581 switch(refinement)
582 {
583 case 1 :
584 std::cout << " 2 <-> 2)" << std::endl;
585 stat1++;
586 break;
587 case 2 :
588 std::cout << "split 1->3)" << std::endl;
589 stat2++;
590 break;
591 case 3 :
592 std::cout << "merge 3->1)" << std::endl;
593 stat3++;
594 break;
595 default :
596 std::cout << "unknown)" << std::endl;
597 break;
598 }
599 std::cout << "\tNew scale : " << t_new << std::endl;
600 MatchScaleLevels(t_up, t_new, intersection_param);
601 MatchScaleLevels(t_new, t_down, intersection_param);
602 std::cout << "\trefinement successful" << std::endl;
603 }
604 else
605 {
606 if (refinement >0) std::cout << "\tAccepting matching but refinement needed (type " << refinement << ")" << std::endl;
607 else std::cout << "\tMatching OK" << std::endl;
608 itMatchUp=matchUp.begin();
609 itMatchDown=matchDown.begin();
610
611 //----- new (sort of) optimised version
612
613 std::cout << "\tResolving matching - version 2" << std::endl;
614 std::vector<std::pair<std::set<int>, std::set<int> > > bifurcTable;
615 std::set<int> set1, set2;
616 for (; itMatchUp!=matchUp.end(); ++itMatchUp)
617 {
618 labelUp=(*itMatchUp).first;
619 std::set<int> setUp=(*itMatchUp).second;
620 if (setUp.size()==0)
621 {
622 //cr�ation d'une paire (apparition)
623 set1=std::set<int>();set1.insert(labelUp);
624 set2=setUp;
625 bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
626 }
627 else
628 {
629 int checkRepeat=0;
630 std::set<int>::iterator itSetUp;
631 for (itSetUp=setUp.begin(); (itSetUp!=setUp.end()&&(checkRepeat==0)); ++itSetUp)
632 {
633 int blobLab=(*itSetUp);
634 if (matchDown[blobLab].size() > 1)
635 checkRepeat=1;
636 }
637 if (checkRepeat==0)
638 {
639 // cr�ation d'une paire
640 set1=std::set<int>();set1.insert(labelUp);
641 set2=setUp;
642 bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
643 }
644 else
645 {
646 // doit on cr�er une nouvelle paire
647 // ou ins�rer dans une existante
648 std::vector<std::pair<std::set<int>, std::set<int> > >::iterator
649 itBifurc=bifurcTable.begin();
650 int flagMerge=0;
651 for (; (itBifurc!=bifurcTable.end())&&(flagMerge==0); ++itBifurc)
652 {
653 std::set<int> setTmp=(*itBifurc).second;
654 std::set<int> setInter;
655 std::set_intersection(setUp.begin(), setUp.end(),
656 setTmp.begin(), setTmp.end(),
657 std::inserter(setInter, setInter.begin()));
658 if (setInter.size()>0)
659 {
660 //on merge ici
661 (*itBifurc).first.insert(labelUp);
662 itSetUp=setUp.begin();
663 for (; itSetUp!=setUp.end(); ++itSetUp)
664 (*itBifurc).second.insert(*itSetUp);
665 flagMerge=1;
666 }
667 }
668 if (flagMerge==0)
669 {
670 //creation d'une paire
671 set1=std::set<int>();set1.insert(labelUp);
672 set2=setUp;
673 bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
674 }
675 }
676 }
677 }
678 for (; itMatchDown!=matchDown.end(); ++itMatchDown)
679 {
680 // cas des disparitions non trait� par la boucle pr�c�dente.
681 labelDown=(*itMatchDown).first;
682 std::set<int> setDown=(*itMatchDown).second;
683 if (setDown.size()==0)
684 {
685 set2=std::set<int>();set2.insert(labelDown);
686 set1=setDown;
687 bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
688 }
689 }
690
691 std::vector< std::pair<std::set<int>, std::set<int> > >::iterator itTable=bifurcTable.begin();
692 std::pair<std::set<int>, std::set<int> > pairSet;
693 //------ end of new version
694
695
696 std::cout << "\t Table -> Blobs and bifurcations" << std::endl;
697 // transformation table -> bifurcations et blobs
698
699 for (; itTable!=bifurcTable.end(); ++itTable) {
700 pairSet=(*itTable);
701 set1=pairSet.first;
702 set2=pairSet.second;
703 if (set2.size()==0) // creation
704 {
705// cout << "creation :" << flush;
706 labelUp = *(set1.begin());
707 blobUp = levelUp->Blob(labelUp);
708 ssBlobUp = GetSSBlobFromGLBlob(blobUp);
709 if ( ssBlobUp != NULL) {
710 bifurc = new Bifurcation<Site>(APPEAR, t_up, t_down);
711 ssBlobUp->SetScaleMin(t_down);
712 bifurc->AddTopBlob(ssBlobUp);
713
714 ssBlobUp->SetBottomBifurcation(bifurc);
715 AddBifurcation(bifurc);
716 // Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
717 // bif = bifurc;
718 // cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
719 }
720 }
721 else if (set1.size()==0) // annihilation
722 {
723// cout << "annihilation :" << flush;
724 labelDown=*(set2.begin());
725 blobDown=levelDown->Blob(labelDown);
726 ssBlobDown=new ScaleSpaceBlob<Site>( _subject, labelMax );
727 labelMax++;
728 ssBlobDown->SetScaleMax(t_down);
729 ssBlobDown->SetScaleMin(t_min);
730 ssBlobDown->AddGreyLevelBlob(blobDown);
731 bifurc = new Bifurcation<Site>(DISAPPEAR, t_up, t_down);
732 bifurc->AddBottomBlob(ssBlobDown);
733 ssBlobDown->SetTopBifurcation(bifurc);
734 AddBifurcation(bifurc);
735 AddBlob(ssBlobDown);
736// Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
737// bif = bifurc;
738// cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
739 }
740 else
741 {
742 if ((set1.size()==1) && (set2.size()==1)) // blob growth
743 {
744 labelUp = *(set1.begin());
745 labelDown = *(set2.begin());
746 blobUp = levelUp->Blob ( labelUp );
747 ssBlobUp = GetSSBlobFromGLBlob ( blobUp );
748 if ( ssBlobUp != NULL) {
749 blobDown = levelDown->Blob ( labelDown );
750 ssBlobUp->AddGreyLevelBlob ( blobDown );
751 }
752 }
753 else if ((set1.size()==1) && (set2.size()==2)) // merge
754 {
755// cout << "merge " << flush ;
756 bifurc = new Bifurcation<Site>(MERGE, t_up, t_down);
757 labelUp = *(set1.begin());
758 blobUp = levelUp->Blob(labelUp);
759 ssBlobUp = GetSSBlobFromGLBlob(blobUp);
760 if ( ssBlobUp != NULL) {
761 ssBlobUp->SetScaleMin(t_down);
762 bifurc->AddTopBlob(ssBlobUp);
763 ssBlobUp->SetBottomBifurcation(bifurc);
764 std::set<int>::iterator itSet2=set2.begin();
765 for (; itSet2!=set2.end(); ++itSet2)
766 {
767 labelDown=*itSet2;
768 blobDown=levelDown->Blob(labelDown);
769 ssBlobDown=new ScaleSpaceBlob<Site>(_subject, labelMax);
770 labelMax++;
771 ssBlobDown->SetScaleMax(t_down);
772 ssBlobDown->SetScaleMin(t_min);
773 ssBlobDown->AddGreyLevelBlob(blobDown);
774 bifurc->AddBottomBlob(ssBlobDown);
775
776 AddBlob(ssBlobDown);
777 ssBlobDown->SetTopBifurcation(bifurc);
778 }
779 AddBifurcation(bifurc);
780 // Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
781 // bif = bifurc;
782 // cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
783 }
784 }
785 else if ( (set1.size()==2) && (set2.size()==1) ) // split
786 {
787// cout << "split " << flush;
788 bifurc = new Bifurcation<Site>(SPLIT, t_up, t_down);
789 labelDown = *(set2.begin());
790 blobDown = levelDown->Blob(labelDown);
791 ssBlobDown = new ScaleSpaceBlob<Site>(_subject, labelMax);
792 labelMax++;
793 ssBlobDown->SetScaleMax(t_down);
794 ssBlobDown->SetScaleMin(t_min);
795 ssBlobDown->AddGreyLevelBlob(blobDown);
796 AddBlob(ssBlobDown);
797 ssBlobDown->SetTopBifurcation(bifurc);
798 bifurc->AddBottomBlob(ssBlobDown);
799 std::set<int>::iterator itSet1=set1.begin();
800 for (; itSet1!=set1.end(); ++itSet1)
801 {
802 labelUp = *itSet1;
803 blobUp = levelUp->Blob(labelUp);
804 ssBlobUp = GetSSBlobFromGLBlob(blobUp);
805 if ( ssBlobUp != NULL) {
806 ssBlobUp->SetScaleMin(t_down);
807 bifurc->AddTopBlob(ssBlobUp);
808 ssBlobUp->SetBottomBifurcation(bifurc);
809 }
810 }
811 AddBifurcation(bifurc);
812// Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
813// bif = bifurc;
814// cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
815
816 }
817 else // complex
818 {
819 bifurc=new Bifurcation<Site>(COMPLEX, t_up, t_down);
820 std::set<int>::iterator itSet2=set2.begin();
821 for (; itSet2!=set2.end(); ++itSet2)
822 {
823 labelDown = *itSet2;
824 blobDown = levelDown->Blob(labelDown);
825 ssBlobDown = new ScaleSpaceBlob<Site>(_subject, labelMax);
826 labelMax++;
827 ssBlobDown->SetScaleMax(t_down);
828 ssBlobDown->SetScaleMin(t_min);
829 ssBlobDown->AddGreyLevelBlob(blobDown);
830 bifurc->AddBottomBlob(ssBlobDown);
831 AddBlob(ssBlobDown);
832 ssBlobDown->SetTopBifurcation(bifurc);
833 }
834 std::set<int>::iterator itSet1=set1.begin();
835 for (; itSet1!=set1.end(); ++itSet1)
836 {
837 labelUp = *itSet1;
838 blobUp = levelUp->Blob(labelUp);
839 ssBlobUp = GetSSBlobFromGLBlob(blobUp);
840 if ( ssBlobUp != NULL) {
841 ssBlobUp->SetScaleMin(t_down);
842 bifurc->AddTopBlob(ssBlobUp);
843 ssBlobUp->SetBottomBifurcation(bifurc);
844 }
845 }
846 AddBifurcation(bifurc);
847 }
848 }
849 }
850 std::cout << "\t OK" << std::endl;
851
852 }
853 }
854
855
856 //---------------------------------------------------------------------
857
858 template<typename Geom, typename Text>
860 const std::string & statName )
861 {
862 ScaleSpaceBlob<Site> *ssBlob;
863 std::map<float, BlobMeasurements > statsM, statsSD;
864 FILE *fileStats;
865 float t;
866 float max_int, mean_int, max_cont, mean_cont, area, areamoy, areavar, tvalue;
867 BlobMeasurements *measurements;
868 GreyLevelBlob<Site> *glBlob, *glBlob1, *glBlob2;
869 int res;
870
871 // Computing SSBlob measurements
872 // statM contains all the stats for normalizing GLBlobs measurements
873 // before integration across scales.
874
875 std::cout << "Computing multiscale measurements for SSBlobs"
876 << std::endl;
877
878 if (statName.length()>0)
879 {
880 if ((fileStats=fopen(statName.c_str(), "r"))==NULL)
881 {
882 std::cerr << "Problem : cannot open stat file " << statName
883 << std::endl;
884 exit(EXIT_FAILURE);
885 }
886 std::cout << "Found stat file, reading it" << std::endl;
887 while (!feof(fileStats))
888 {
889 res = fscanf(fileStats, "%f\n", &t); // !!! this defines the stat file format (simple)
890 if ( res != 1 )
891 {
892 std::cerr << "file corrupted" << std::endl;
893 }
894 res = fscanf(fileStats, "%f %f %f %f %f", &max_int, &mean_int, &max_cont, &mean_cont, &area);
895 if ( res != 5 )
896 {
897 std::cerr << "file corrupted" << std::endl;
898 }
899 measurements=new BlobMeasurements(max_int, mean_int, max_cont, mean_cont, area);
900 statsM.insert(std::pair<float, BlobMeasurements>(t, *measurements));
901 res = fscanf(fileStats, "%f %f %f %f %f", &max_int, &mean_int, &max_cont, &mean_cont, &area);
902 if ( res != 5 )
903 {
904 std::cerr << "file corrupted" << std::endl;
905 }
906 measurements=new BlobMeasurements(max_int, mean_int, max_cont, mean_cont, area);
907 statsSD.insert(std::pair<float, BlobMeasurements>(t, *measurements));
908 }
909 fclose(fileStats);
910 typename std::list<ScaleSpaceBlob<Site>*>::iterator itSSBlobs=blobList.begin();
911 std::cout << "processing GLblobs for normalisation" << std::endl;
912 for (; itSSBlobs!=blobList.end(); ++itSSBlobs)
913 {
914 ssBlob=*itSSBlobs;
915 typename std::list<GreyLevelBlob<Site>*>::iterator itglb;
916
917 for (itglb=ssBlob->glBlobs.begin();itglb!=ssBlob->glBlobs.end();++itglb)
918 {
919 glBlob=*itglb;
920 t=glBlob->GetScale();
921 if (statsM.find(t)!=statsM.end())
922 {
923 glBlob->measurements.maxIntensity=(glBlob->measurements.maxIntensity-statsM[t].maxIntensity)/statsSD[t].maxIntensity;
924 glBlob->measurements.meanIntensity=(glBlob->measurements.meanIntensity-statsM[t].meanIntensity)/statsSD[t].meanIntensity;
925 glBlob->measurements.maxContrast=(glBlob->measurements.maxContrast-statsM[t].maxContrast)/statsSD[t].maxContrast;
926 glBlob->measurements.meanContrast=(glBlob->measurements.meanContrast-statsM[t].meanContrast)/statsSD[t].meanContrast;
927 glBlob->measurements.area=(glBlob->measurements.area-statsM[t].area)/statsSD[t].area;
928 if (glBlob->measurements.maxIntensity > 0)
930 else
932 if (glBlob->measurements.meanIntensity > 0)
934 else
936 if (glBlob->measurements.maxContrast > 0)
938 else
940 if (glBlob->measurements.meanContrast > 0)
942 else
944 if (glBlob->measurements.area > 0)
945 glBlob->measurements.area=1+glBlob->measurements.area;
946 else
947 glBlob->measurements.area=exp(glBlob->measurements.area);
948 }
949 else
950 {
951 std::cerr << "Problem during measurements normalisation !" << std::endl;
952 std::cerr << "Stat file " << statName << " incomplete, can't find scale " << t << std::endl;
953 exit(EXIT_FAILURE);
954 }
955 }
956 }
957 }
958 typename std::list< ScaleSpaceBlob<Site>* >::iterator itSSBlobs=blobList.begin();
959 std::cout << "Computing multi-scale measurements" << std::endl;
960
961 for (; itSSBlobs!=blobList.end(); ++itSSBlobs)
962 {
963 ssBlob=*itSSBlobs;
964 typename std::list<GreyLevelBlob<Site>*>::iterator itGLBlobs=ssBlob->glBlobs.begin();
965 BlobMeasurements measure;
966 float t2,t1, dt;
967 max_int=0; mean_int=0; max_cont=0; mean_cont=0; area=0; areamoy=0; areavar=0; tvalue=0;
968 glBlob1=*itGLBlobs;
969 t1=glBlob1->GetScale();
970 dt=log(ssBlob->TopBifurcation()->tMid())-log(t1);
971 // DEBUG
972 if (dt< 0)
973 std::cout << "PROBLEME DE DT<0" << std::endl;
974 // FIN DEBUG
975 max_int += dt * (glBlob1->measurements.maxIntensity);
976 mean_int += dt * (glBlob1->measurements.meanIntensity);
977 max_cont += dt * (glBlob1->measurements.maxContrast);
978 mean_cont += dt * (glBlob1->measurements.meanContrast);
979 area += dt * (glBlob1->measurements.area);
980 tvalue += dt * (glBlob1->measurements.t);
981
982
983 ++itGLBlobs;
984 for (; itGLBlobs!=ssBlob->glBlobs.end(); ++itGLBlobs)
985 {
986 glBlob2=*itGLBlobs;
987 t2=glBlob2->GetScale();
988 dt=log(t1)-log(t2);
989 max_int += dt * (glBlob2->measurements.maxIntensity + glBlob1->measurements.maxIntensity)/2.0;
990 mean_int += dt * (glBlob2->measurements.meanIntensity + glBlob1->measurements.meanIntensity)/2.0;
991 max_cont += dt * (glBlob2->measurements.maxContrast + glBlob1->measurements.maxContrast)/2.0;
992 mean_cont += dt * (glBlob2->measurements.meanContrast + glBlob1->measurements.meanContrast)/2.0;
993 area += dt * (glBlob2->measurements.area + glBlob1->measurements.area)/2.0;
994 tvalue += dt * (glBlob2->measurements.t + glBlob1->measurements.t)/2.0;
995
996
997 glBlob1=glBlob2;
998 t1=t2;
999 }
1000 t1=glBlob1->GetScale();
1001 dt=log(t1) - log(ssBlob->BottomBifurcation()->tMid());
1002 // DEBUG
1003 if (dt< 0)
1004 std::cout << "PROBLEME DE DT<0" << std::endl;
1005 // FIN DEBUG
1006 max_int += dt * (glBlob1->measurements.maxIntensity);
1007 mean_int += dt * (glBlob1->measurements.meanIntensity);
1008 max_cont += dt * (glBlob1->measurements.maxContrast);
1009 mean_cont += dt * (glBlob1->measurements.meanContrast);
1010 area += dt * (glBlob1->measurements.area);
1011 tvalue += dt * (glBlob1->measurements.t);
1012
1013
1014 // Getting T-Value From Original Map...
1015
1016 TexturedData<Geom, Text> ima=scaleSpace()->Scale(0.0)->Data();
1017
1018 glBlob1=ssBlob->GlBlobRep();
1019 std::set<Site,ltstr_p3d<Site> > pixels;
1020 pixels=glBlob1->GetListePoints();
1021
1022
1023 typename std::set<Site, ltstr_p3d<Site> >::iterator itPix;
1024 float tvmax=-100.0;
1025 for (itPix=pixels.begin();itPix!=pixels.end();itPix++){
1026 if (ima.intensity(*itPix) > tvmax )
1027 tvmax = ima.intensity(*itPix);
1028 }
1029
1030
1031 measure=BlobMeasurements(max_int, mean_int, max_cont,mean_cont, area, tvalue, tvmax);
1032
1033 areamoy +=area;
1034 areavar +=area*area;
1035 ssBlob->SetMeasurements(measure);
1036 }
1037 }
1038
1039
1040}
1041
1042#endif
void AddTopBlob(ScaleSpaceBlob< T > *blob)
void AddBottomBlob(ScaleSpaceBlob< T > *blob)
Class for grey-level blobs Templated with respect to the type of points: TypeSite<carto::VolumeRef<T>...
std::set< T, ltstr_p3d< T > > & GetListePoints()
BlobMeasurements measurements
ScaleSpaceBlob< Site > * GetSSBlobFromGLBlob(GreyLevelBlob< Site > *glBlob)
TexturedData< Geom, Text > * _mask
void AddBlob(ScaleSpaceBlob< Site > *blob)
std::list< Bifurcation< Site > * > BifurcationList()
std::string Subject()
void MatchScaleLevels(float t_up, float t_down, uint intersection_param)
void setType(int t)
ScaleSpaceBlob< Site > * Blob(int label)
std::list< ScaleSpaceBlob< Site > * > blobList
void SetScaleSpace(ScaleSpace< Geom, Text > *scaleSpace)
ScaleSpace< Geom, Text > * _scaleSpace
void ComputePrimalSketch(float tmin, float tMax, const std::string &statFile="", uint intersection_param=10)
std::list< Bifurcation< Site > * > bifurcationList
std::string _subject
void ComputeBlobMeasurements(const std::string &statName="")
SiteType< Geom >::type Site
ScaleSpace< Geom, Text > * scaleSpace()
void AddBifurcation(Bifurcation< Site > *bifurcation)
PrimalSketch(const std::string &subject, ScaleSpace< Geom, Text > *scaleSpace, int type)
PrimalSketch(const std::string &subject, int type)
PrimalSketch(const std::string &subject, ScaleSpace< Geom, Text > *scaleSpace, TexturedData< Geom, Text > *mask, int type)
TexType< Text >::type Val
std::list< ScaleSpaceBlob< Site > * > BlobSet()
void DetectBlobs(TexturedData< Geom, Text > *mask, char *stats=0)
GreyLevelBlob< Site > * Blob(int label)
Definition scaleLevel.h:73
std::map< int, GreyLevelBlob< Site > * > BlobList()
Definition scaleLevel.h:69
Bifurcation< T > * TopBifurcation()
void AddGreyLevelBlob(GreyLevelBlob< T > *blob)
GreyLevelBlob< T > * GlBlobRep()
void SetScaleMin(float t)
void SetBottomBifurcation(Bifurcation< T > *bif)
Bifurcation< T > * BottomBifurcation()
void SetMeasurements(BlobMeasurements meas)
std::list< GreyLevelBlob< T > * > glBlobs
void SetTopBifurcation(Bifurcation< T > *bif)
void SetScaleMax(float t)
BucketMap< Void > * mask(const BucketMap< Void > &src, const BucketMap< Void > &m, bool intersect=true)
PrimalsketchDataType
@ SURFACE
unsigned int uint