A.I.M.S algorithms


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>
46 #include <aims/mesh/surfacegen.h>
47 
48 
49 namespace aims
50 {
51 
52  enum PrimalsketchDataType
53  {
54  SURFACE,
55  IMAGE
56  };
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,
93  TexturedData<Geom, Text> *mask, int type)
94  : _subject(subject) , _scaleSpace(scaleSpace), labelMax(0), _mask(mask), _type(type)
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); }
117  void AddBifurcation( Bifurcation<Site> *bifurcation )
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  int nbBlobs;
220 
221  levelUp->DetectBlobs( _mask );
222 
223  nbBlobs = levelUp->nbBlobs();
224 
225  GreyLevelBlob<Site> *glBlob;
226  ScaleSpaceBlob<Site> *ssBlob;
227  Bifurcation<Site> *bifurc;
228 
229  // Scale-space blob initialisation
230  std::map<int, GreyLevelBlob<Site> *>
231  listeGLB=levelUp->BlobList();
232 
233  typename std::map<int, GreyLevelBlob<Site> *>::iterator
234  itGL = listeGLB.begin();
235 
236  for ( ; itGL != listeGLB.end() ; ++itGL ) {
237  glBlob=(*itGL).second;
238  ssBlob=new ScaleSpaceBlob<Site>(_subject, labelMax);
239  ssBlob->SetScaleMax(*itScaleUp);
240  ssBlob->SetScaleMin(t_min);
241  ssBlob->AddGreyLevelBlob(glBlob);
242  bifurc = new Bifurcation<Site>(DISAPPEAR, t_bif, t_max);
243  bifurc->AddBottomBlob(ssBlob);
244  ssBlob->SetTopBifurcation(bifurc);
245  AddBlob(ssBlob);
246  AddBifurcation(bifurc);
247  labelMax++;
248  }
249 
250  // Primal Sketch Computation, Top To Bottom
251 
252  for ( ; (*itScaleUp) > tmin ; ++itScaleUp ) {
253  itScaleDown=itScaleUp;
254  itScaleDown++;
255  t_up = *itScaleUp;
256  t_down = *itScaleDown;
257  MatchScaleLevels(t_up, t_down, intersection_param);
258  }
259 
260  // "closing" the primal sketch at bottom
261 
262  listeGLB = levelDown->BlobList();
263  itGL = listeGLB.begin();
264  t_bif = t_min / 2.0;
265  for ( ; itGL != listeGLB.end() ; ++itGL ) {
266  glBlob = (*itGL).second;
267  ssBlob = GetSSBlobFromGLBlob(glBlob);
268  if ( ssBlob != NULL) {
269  bifurc = new Bifurcation<Site>(APPEAR, t_min, t_bif);
270  ssBlob->SetBottomBifurcation(bifurc);
271  bifurc->AddTopBlob(ssBlob);
272  AddBifurcation(bifurc);
273  }
274  }
275 
276  // blob updates
277  typename std::list<ScaleSpaceBlob<Site>*>::iterator
278  itBlob = blobList.begin();
279  for ( ; itBlob != blobList.end() ; ++itBlob ) {
280  (*itBlob)->ComputeLifeTime();
281  (*itBlob)->ComputeScaleRep();
282  }
283  // blob (multiscale) measurements are computed here.
284 
285  ComputeBlobMeasurements(statFile);
286  std::cout << "ssblobs avant : " << blobList.size() << std::endl;
287 
288  // ELAGAGE
289  typename std::list<Bifurcation<Site>*>::iterator bifit = bifurcationList.begin();
290  //cout << "bifurcation number:" << bifurcationList.size() << endl;
291  uint merge_count = 0;
292 
293  for ( bifit = bifurcationList.begin() ; bifit != bifurcationList.end() ; bifit++ )
294  if ((*bifit)->Type() == MERGE)
295  merge_count ++;
296  //cout << "MERGECOUNT:" << merge_count << endl;
297 
298 
299  for ( bifit = bifurcationList.begin() ; bifit != bifurcationList.end() ; bifit++ ) {
300  if ( (*bifit)->Type() == APPEAR ) {
301  // std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->TopBlobs();
302  // typename std::list<ScaleSpaceBlob<Site>*>::iterator topit=blist.begin();
303  // ScaleSpaceBlob<Site> *ssblob1;
304  // ssblob1 = *topit;
305  // if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1) {
306  // int label1 = ssblob1->Label();
307  // for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
308  // if ((*itBlob)->Label() == label1) {
309  // blobList.erase(itBlob);
310  // break;
311  // }
312  //
313  // }
314 
315  }
316  else if ( (*bifit)->Type() == DISAPPEAR ) {
317  // std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->BottomBlobs();
318  // typename std::list<ScaleSpaceBlob<Site>*>::iterator bottomit=blist.begin();
319  // ScaleSpaceBlob<Site> *ssblob1;
320  // ssblob1 = *bottomit;
321  // if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1) {
322  // int label1 = ssblob1->Label();
323  // for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
324  // if ((*itBlob)->Label() == label1) {
325  // blobList.erase(itBlob);
326  // break;
327  // }
328  // }
329  }
330  else if ( (*bifit)->Type() == MERGE ) {
331  //cout << "merge" << (*bifit)->BottomBlobs().size()<< ";" << (*bifit)->TopBlobs().size() << endl;
332  std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->BottomBlobs();
333  typename std::list<ScaleSpaceBlob<Site>*>::iterator bottomit=blist.begin();
334  ScaleSpaceBlob<Site> *ssblob1, *ssblob2;
335  ssblob1 = *bottomit;
336  bottomit++;
337  ssblob2 = *bottomit;
338  /* float area1 = (*ssblob1).GetMeasurements().area;
339  float area2 = (*ssblob2).GetMeasurements().area; */
340  //cout << /*area1 << ";" << area2 << ";" <<*/ max(area1,area2)/min(area1,area2) << endl;
341  if ((*ssblob2).GlBlobRep()->GetListePoints().size() == 1) {
342  //cout << "ELAG" << endl;
343  // LE BLOB2 EST PETIT : ON FUSIONNE BLOB1 AVEC TOPBLOB ET TOPBIF DE BLOB2 DEVIENT DISPARITION
344 
345  ssblob2->TopBifurcation()->setType(DISAPPEAR);
346  ssblob2->TopBifurcation()->TopBlobs().clear();
347  ssblob2->TopBifurcation()->BottomBlobs().clear();
348  ssblob2->TopBifurcation()->AddBottomBlob(ssblob2);
349 
350  blist = (*bifit)->TopBlobs();
351  ScaleSpaceBlob<Site> *topblob = *(blist.begin());
352  // RAJOUT DES GLB DE SSBLOB1
353  std::list<GreyLevelBlob<Site>*> glb = ssblob1->glBlobs;
354  typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
355  for (; glbit != glb.end() ; glbit++)
356  topblob->AddGreyLevelBlob(*glbit);
357 
358  topblob->SetBottomBifurcation(ssblob1->BottomBifurcation());
359  topblob->SetScaleMin(ssblob1->ScaleMin());
360 
361  int label1 = ssblob1->Label();
362  for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
363  if ((*itBlob)->Label() == label1) {
364  blobList.erase(itBlob);
365  break;
366  }
367  }
368  else if ( (*ssblob1).GlBlobRep()->GetListePoints().size() == 1 ) {
369  // LE BLOB1 EST PETIT : ON FUSIONNE BLOB2 AVEC TOPBLOB ET TOPBIF DE BLOB1 DEVIENT DISPARITION
370  //cout << "ELAG" << endl;
371  ssblob1->TopBifurcation()->setType(DISAPPEAR);
372  ssblob1->TopBifurcation()->TopBlobs().clear();
373  ssblob1->TopBifurcation()->BottomBlobs().clear();
374  ssblob1->TopBifurcation()->AddBottomBlob(ssblob1);
375 
376  blist = (*bifit)->TopBlobs();
377  ScaleSpaceBlob<Site> *topblob = *(blist.begin());
378  // RAJOUT DES GLB DE SSBLOB2
379  std::list<GreyLevelBlob<Site>*> glb = ssblob2->glBlobs;
380  typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
381  for (; glbit != glb.end() ; glbit++)
382  topblob->AddGreyLevelBlob(*glbit);
383 
384  topblob->SetBottomBifurcation(ssblob2->BottomBifurcation());
385  topblob->SetScaleMin(ssblob2->ScaleMin());
386 
387  int label2 = ssblob2->Label();
388  for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
389  if ((*itBlob)->Label() == label2) {
390  blobList.erase(itBlob);
391  break;
392  }
393  }
394 
395  }
396  else if ((*bifit)->Type() == SPLIT) {
397 
398  //cout << "split" << (*bifit)->TopBlobs().size()<< ";" << (*bifit)->BottomBlobs().size() << endl;
399  std::list<ScaleSpaceBlob<Site>*> blist = (*bifit)->TopBlobs();
400  typename std::list<ScaleSpaceBlob<Site>*>::iterator topit=blist.begin();
401  ScaleSpaceBlob<Site> *ssblob1, *ssblob2;
402  ssblob1 = *topit;
403  topit++;
404  ssblob2 = *topit;
405  /* float area1 = (*ssblob1).GetMeasurements().area;
406  float area2 = (*ssblob2).GetMeasurements().area; */
407  // cout << /*area1 << ";" << area2 << ";" <<*/ max(area1,area2)/min(area1,area2) << endl;
408  if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1) {
409  // LE BLOB2 EST PETIT : ON FUSIONNE BLOB1 AVEC BOTTOMBLOB ET BOTTOMBIF DE BLOB2 DEVIENT DISPARITION
410 
411  ssblob2->BottomBifurcation()->setType(DISAPPEAR);
412  ssblob2->BottomBifurcation()->BottomBlobs().clear();
413  ssblob2->BottomBifurcation()->TopBlobs().clear();
414  ssblob2->BottomBifurcation()->AddTopBlob(ssblob2);
415 
416  blist = (*bifit)->BottomBlobs();
417  ScaleSpaceBlob<Site> *topblob = *(blist.begin());
418  // RAJOUT DES GLB DE SSBLOB1
419  std::list<GreyLevelBlob<Site>*> glb = ssblob1->glBlobs;
420  typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
421  for (; glbit != glb.end() ; glbit++)
422  topblob->AddGreyLevelBlob(*glbit);
423 
424  topblob->SetTopBifurcation(ssblob1->TopBifurcation());
425  topblob->SetScaleMin(ssblob1->ScaleMin());
426 
427  int label1 = ssblob1->Label();
428  for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob)
429  if ((*itBlob)->Label() == label1) {
430  blobList.erase(itBlob);
431  break;
432  }
433  }
434  else if ((*ssblob1).GlBlobRep()->GetListePoints().size() == 1){
435  // LE BLOB1 EST PETIT : ON FUSIONNE BLOB2 AVEC BOTTOMBLOB ET BOTTOMBIF DE BLOB1 DEVIENT DISPARITION
436 
437  ssblob1->BottomBifurcation()->setType(DISAPPEAR);
438  ssblob1->BottomBifurcation()->BottomBlobs().clear();
439  ssblob1->BottomBifurcation()->TopBlobs().clear();
440  ssblob1->BottomBifurcation()->AddTopBlob(ssblob1);
441 
442  blist = (*bifit)->BottomBlobs();
443  ScaleSpaceBlob<Site> *topblob = * ( blist.begin() );
444  // RAJOUT DES GLB DE SSBLOB2
445  std::list<GreyLevelBlob<Site>*> glb = ssblob2->glBlobs;
446  typename std::list<GreyLevelBlob<Site>*>::iterator glbit = glb.begin();
447  for ( ; glbit != glb.end() ; glbit++ )
448  topblob->AddGreyLevelBlob( *glbit );
449 
450  topblob->SetTopBifurcation ( ssblob2->TopBifurcation() );
451  topblob->SetScaleMin ( ssblob2->ScaleMin() );
452 
453  int label2 = ssblob2->Label();
454  for ( itBlob = blobList.begin() ; itBlob != blobList.end() ; ++itBlob )
455  if ( (*itBlob)->Label() == label2 ) {
456  blobList.erase ( itBlob );
457  break;
458  }
459  }
460  }
461  }
462  merge_count = 0;
463  bifit = bifurcationList.begin();
464  for (bifit = bifurcationList.begin() ; bifit != bifurcationList.end() ; bifit++)
465  if ((*bifit)->Type() == MERGE)
466  merge_count ++;
467  //cout << "MERGECOUNT:" << merge_count << endl;
468 
469  // blob updates
470  std::cout << "ssblobs après : " << blobList.size() << std::endl;
471 
472  for (itBlob=blobList.begin(); itBlob != blobList.end(); ++itBlob) {
473  (*itBlob)->ComputeLifeTime();
474  (*itBlob)->ComputeScaleRep();
475  }
476 
477  ComputeBlobMeasurements(statFile);
478  }
479 
480  //---------------------------------------------------------------------
481 
482  template<typename Geom, typename Text>
483  void PrimalSketch<Geom,Text>::MatchScaleLevels(float t_up, float t_down, uint intersection_param)
484  {
485  float t, t_new;
486  int stat1=0, stat2=0, stat3=0;
487 
488  std::cout << "Matching " << t_up << " <-> " << t_down << std::endl;
489 
490  ScaleLevel<Geom,Text> *levelUp=_scaleSpace->Scale(t_up);
491  ScaleLevel<Geom,Text> *levelDown=_scaleSpace->Scale(t_down);
492  ScaleSpaceBlob<Site> *ssBlobUp, *ssBlobDown;
493  Bifurcation<Site> *bifurc;
494 
495  std::cout << "Blob detection at scale " << t_up << std::endl;
496  levelUp->DetectBlobs(_mask);
497  std::cout << "Blob detection at scale " << t_down << std::endl;
498  levelDown->DetectBlobs(_mask);
499 
500  std::cout << "\tt=" << t_up << " : " << levelUp->nbBlobs() << " blobs" << std::endl;
501  std::cout << "\tt=" << t_down << " : " << levelDown->nbBlobs() << " blobs" << std::endl;
502 
503  std::map<int, GreyLevelBlob<Site> *> listeGLB_up=levelUp->BlobList();
504  std::map<int, GreyLevelBlob<Site> *> listeGLB_down=levelDown->BlobList();
505 
506  // On essaye la maniere (simple mais) brutale
507 
508  std::cout << "\tComputing potential links" << std::endl;
509 
510  std::map<int, std::set<int> > matchUp;
511  std::map<int, std::set<int> > matchDown;
512  typename std::map<int, GreyLevelBlob<Site> *>::iterator
513  iteratorUp;
514  typename std::map<int, GreyLevelBlob<Site> *>::iterator
515  iteratorDown;
516  GreyLevelBlob<Site> *blobUp, *blobDown;
517  int labelUp, labelDown;
518 
519  iteratorUp=listeGLB_up.begin();
520  for (; iteratorUp!=listeGLB_up.end(); ++iteratorUp)
521  {
522  blobUp=(*iteratorUp).second; labelUp=(*iteratorUp).first;
523  matchUp[labelUp]=std::set<int>();
524  std::set<Site,ltstr_p3d<Site> > pointsUp=blobUp->GetListePoints();
525  iteratorDown=listeGLB_down.begin();
526  for (; iteratorDown!=listeGLB_down.end(); ++iteratorDown)
527  {
528  blobDown=(*iteratorDown).second; labelDown=(*iteratorDown).first;
529  if (matchDown.find(labelDown) == matchDown.end())
530  matchDown[labelDown]=std::set<int>();
531  // blobs intersect eachother
532  std::set<Site,ltstr_p3d<Site> > pointsInter;
533  std::set<Site,ltstr_p3d<Site> > pointsDown=blobDown->GetListePoints();
534  std::set_intersection(pointsDown.begin(), pointsDown.end(),
535  pointsUp.begin(), pointsUp.end(),
536  std::inserter(pointsInter, pointsInter.begin()),
537  ltstr_p3d<Site>() );
538 // if (2*pointsInter.size()/(pointsUp.size() + pointsDown.size()) >0.5)
539  if (pointsInter.size()>intersection_param)
540  {
541  matchUp[labelUp].insert(labelDown);
542  matchDown[labelDown].insert(labelUp);
543  }
544  }
545  }
546 
547  //-------------------------------------------------------------------
548  std::cout << "\tChecking potential refinements" << std::endl;
549  std::map<int, std::set<int> >::iterator itMatchUp=matchUp.begin();
550  std::map<int, std::set<int> >::iterator itMatchDown=matchDown.begin();
551  int refinement=0;
552  for (; itMatchUp!=matchUp.end(); ++itMatchUp)
553  {
554  int nUp=(*itMatchUp).second.size();
555  if (nUp > 2)
556  {
557  refinement=2;
558  } //---------------------------------------------------------------------
559  else if (nUp == 2)
560  {
561  std::set<int> setUp=(*itMatchUp).second;
562  std::set<int>::iterator itUp=setUp.begin();
563  for (; itUp!=setUp.end(); ++itUp)
564  if (matchDown[*itUp].size()>1)
565  refinement=1;
566  }
567  if (refinement>0) break;
568  }
569  if (refinement==0)
570  {
571  for (; itMatchDown!=matchDown.end(); ++itMatchDown)
572  {
573  if ((*itMatchDown).second.size() > 2)
574  refinement=3;
575  if (refinement>0) break;
576  }
577  }
578  t=exp((log(t_up) + log(t_down))/2.0); t_new=0;
579  int i=0;
580  while ((i*0.1) <= t) i++;
581  t_new=(i-1)*0.1;
582  if ((refinement>0) && ((t_up-t_new)>min_delta_t) && ((t_new-t_down)>min_delta_t))
583  {
584  std::cout << "\tRefinement needed (type " ;
585  switch(refinement)
586  {
587  case 1 :
588  std::cout << " 2 <-> 2)" << std::endl;
589  stat1++;
590  break;
591  case 2 :
592  std::cout << "split 1->3)" << std::endl;
593  stat2++;
594  break;
595  case 3 :
596  std::cout << "merge 3->1)" << std::endl;
597  stat3++;
598  break;
599  default :
600  std::cout << "unknown)" << std::endl;
601  break;
602  }
603  std::cout << "\tNew scale : " << t_new << std::endl;
604  MatchScaleLevels(t_up, t_new, intersection_param);
605  MatchScaleLevels(t_new, t_down, intersection_param);
606  std::cout << "\trefinement successful" << std::endl;
607  }
608  else
609  {
610  if (refinement >0) std::cout << "\tAccepting matching but refinement needed (type " << refinement << ")" << std::endl;
611  else std::cout << "\tMatching OK" << std::endl;
612  itMatchUp=matchUp.begin();
613  itMatchDown=matchDown.begin();
614 
615  //----- new (sort of) optimised version
616 
617  std::cout << "\tResolving matching - version 2" << std::endl;
618  std::vector<std::pair<std::set<int>, std::set<int> > > bifurcTable;
619  std::set<int> set1, set2;
620  for (; itMatchUp!=matchUp.end(); ++itMatchUp)
621  {
622  labelUp=(*itMatchUp).first;
623  std::set<int> setUp=(*itMatchUp).second;
624  if (setUp.size()==0)
625  {
626  //cr�ation d'une paire (apparition)
627  set1=std::set<int>();set1.insert(labelUp);
628  set2=setUp;
629  bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
630  }
631  else
632  {
633  int checkRepeat=0;
634  std::set<int>::iterator itSetUp;
635  for (itSetUp=setUp.begin(); (itSetUp!=setUp.end()&&(checkRepeat==0)); ++itSetUp)
636  {
637  int blobLab=(*itSetUp);
638  if (matchDown[blobLab].size() > 1)
639  checkRepeat=1;
640  }
641  if (checkRepeat==0)
642  {
643  // cr�ation d'une paire
644  set1=std::set<int>();set1.insert(labelUp);
645  set2=setUp;
646  bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
647  }
648  else
649  {
650  // doit on cr�er une nouvelle paire
651  // ou ins�rer dans une existante
652  std::vector<std::pair<std::set<int>, std::set<int> > >::iterator
653  itBifurc=bifurcTable.begin();
654  int flagMerge=0;
655  for (; (itBifurc!=bifurcTable.end())&&(flagMerge==0); ++itBifurc)
656  {
657  std::set<int> setTmp=(*itBifurc).second;
658  std::set<int> setInter;
659  std::set_intersection(setUp.begin(), setUp.end(),
660  setTmp.begin(), setTmp.end(),
661  std::inserter(setInter, setInter.begin()));
662  if (setInter.size()>0)
663  {
664  //on merge ici
665  (*itBifurc).first.insert(labelUp);
666  itSetUp=setUp.begin();
667  for (; itSetUp!=setUp.end(); ++itSetUp)
668  (*itBifurc).second.insert(*itSetUp);
669  flagMerge=1;
670  }
671  }
672  if (flagMerge==0)
673  {
674  //creation d'une paire
675  set1=std::set<int>();set1.insert(labelUp);
676  set2=setUp;
677  bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
678  }
679  }
680  }
681  }
682  for (; itMatchDown!=matchDown.end(); ++itMatchDown)
683  {
684  // cas des disparitions non trait� par la boucle pr�c�dente.
685  labelDown=(*itMatchDown).first;
686  std::set<int> setDown=(*itMatchDown).second;
687  if (setDown.size()==0)
688  {
689  set2=std::set<int>();set2.insert(labelDown);
690  set1=setDown;
691  bifurcTable.push_back(std::pair<std::set<int>, std::set<int> >(set1, set2));
692  }
693  }
694 
695  std::vector< std::pair<std::set<int>, std::set<int> > >::iterator itTable=bifurcTable.begin();
696  std::pair<std::set<int>, std::set<int> > pairSet;
697  //------ end of new version
698 
699 
700  std::cout << "\t Table -> Blobs and bifurcations" << std::endl;
701  // transformation table -> bifurcations et blobs
702 
703  for (; itTable!=bifurcTable.end(); ++itTable) {
704  pairSet=(*itTable);
705  set1=pairSet.first;
706  set2=pairSet.second;
707  if (set2.size()==0) // creation
708  {
709 // cout << "creation :" << flush;
710  labelUp = *(set1.begin());
711  blobUp = levelUp->Blob(labelUp);
712  ssBlobUp = GetSSBlobFromGLBlob(blobUp);
713  if ( ssBlobUp != NULL) {
714  bifurc = new Bifurcation<Site>(APPEAR, t_up, t_down);
715  ssBlobUp->SetScaleMin(t_down);
716  bifurc->AddTopBlob(ssBlobUp);
717 
718  ssBlobUp->SetBottomBifurcation(bifurc);
719  AddBifurcation(bifurc);
720  // Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
721  // bif = bifurc;
722  // cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
723  }
724  }
725  else if (set1.size()==0) // annihilation
726  {
727 // cout << "annihilation :" << flush;
728  labelDown=*(set2.begin());
729  blobDown=levelDown->Blob(labelDown);
730  ssBlobDown=new ScaleSpaceBlob<Site>( _subject, labelMax );
731  labelMax++;
732  ssBlobDown->SetScaleMax(t_down);
733  ssBlobDown->SetScaleMin(t_min);
734  ssBlobDown->AddGreyLevelBlob(blobDown);
735  bifurc = new Bifurcation<Site>(DISAPPEAR, t_up, t_down);
736  bifurc->AddBottomBlob(ssBlobDown);
737  ssBlobDown->SetTopBifurcation(bifurc);
738  AddBifurcation(bifurc);
739  AddBlob(ssBlobDown);
740 // Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
741 // bif = bifurc;
742 // cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
743  }
744  else
745  {
746  if ((set1.size()==1) && (set2.size()==1)) // blob growth
747  {
748  labelUp = *(set1.begin());
749  labelDown = *(set2.begin());
750  blobUp = levelUp->Blob ( labelUp );
751  ssBlobUp = GetSSBlobFromGLBlob ( blobUp );
752  if ( ssBlobUp != NULL) {
753  blobDown = levelDown->Blob ( labelDown );
754  ssBlobUp->AddGreyLevelBlob ( blobDown );
755  }
756  }
757  else if ((set1.size()==1) && (set2.size()==2)) // merge
758  {
759 // cout << "merge " << flush ;
760  bifurc = new Bifurcation<Site>(MERGE, t_up, t_down);
761  labelUp = *(set1.begin());
762  blobUp = levelUp->Blob(labelUp);
763  ssBlobUp = GetSSBlobFromGLBlob(blobUp);
764  if ( ssBlobUp != NULL) {
765  ssBlobUp->SetScaleMin(t_down);
766  bifurc->AddTopBlob(ssBlobUp);
767  ssBlobUp->SetBottomBifurcation(bifurc);
768  std::set<int>::iterator itSet2=set2.begin();
769  for (; itSet2!=set2.end(); ++itSet2)
770  {
771  labelDown=*itSet2;
772  blobDown=levelDown->Blob(labelDown);
773  ssBlobDown=new ScaleSpaceBlob<Site>(_subject, labelMax);
774  labelMax++;
775  ssBlobDown->SetScaleMax(t_down);
776  ssBlobDown->SetScaleMin(t_min);
777  ssBlobDown->AddGreyLevelBlob(blobDown);
778  bifurc->AddBottomBlob(ssBlobDown);
779 
780  AddBlob(ssBlobDown);
781  ssBlobDown->SetTopBifurcation(bifurc);
782  }
783  AddBifurcation(bifurc);
784  // Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
785  // bif = bifurc;
786  // cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
787  }
788  }
789  else if ( (set1.size()==2) && (set2.size()==1) ) // split
790  {
791 // cout << "split " << flush;
792  bifurc = new Bifurcation<Site>(SPLIT, t_up, t_down);
793  labelDown = *(set2.begin());
794  blobDown = levelDown->Blob(labelDown);
795  ssBlobDown = new ScaleSpaceBlob<Site>(_subject, labelMax);
796  labelMax++;
797  ssBlobDown->SetScaleMax(t_down);
798  ssBlobDown->SetScaleMin(t_min);
799  ssBlobDown->AddGreyLevelBlob(blobDown);
800  AddBlob(ssBlobDown);
801  ssBlobDown->SetTopBifurcation(bifurc);
802  bifurc->AddBottomBlob(ssBlobDown);
803  std::set<int>::iterator itSet1=set1.begin();
804  for (; itSet1!=set1.end(); ++itSet1)
805  {
806  labelUp = *itSet1;
807  blobUp = levelUp->Blob(labelUp);
808  ssBlobUp = GetSSBlobFromGLBlob(blobUp);
809  if ( ssBlobUp != NULL) {
810  ssBlobUp->SetScaleMin(t_down);
811  bifurc->AddTopBlob(ssBlobUp);
812  ssBlobUp->SetBottomBifurcation(bifurc);
813  }
814  }
815  AddBifurcation(bifurc);
816 // Bifurcation<SiteType<AimsSurface<3, Void> >::type> *bif;
817 // bif = bifurc;
818 // cout << bif->TopBlobs().size() << " " << bif->BottomBlobs().size() << " - " << flush;
819 
820  }
821  else // complex
822  {
823  bifurc=new Bifurcation<Site>(COMPLEX, t_up, t_down);
824  std::set<int>::iterator itSet2=set2.begin();
825  for (; itSet2!=set2.end(); ++itSet2)
826  {
827  labelDown = *itSet2;
828  blobDown = levelDown->Blob(labelDown);
829  ssBlobDown = new ScaleSpaceBlob<Site>(_subject, labelMax);
830  labelMax++;
831  ssBlobDown->SetScaleMax(t_down);
832  ssBlobDown->SetScaleMin(t_min);
833  ssBlobDown->AddGreyLevelBlob(blobDown);
834  bifurc->AddBottomBlob(ssBlobDown);
835  AddBlob(ssBlobDown);
836  ssBlobDown->SetTopBifurcation(bifurc);
837  }
838  std::set<int>::iterator itSet1=set1.begin();
839  for (; itSet1!=set1.end(); ++itSet1)
840  {
841  labelUp = *itSet1;
842  blobUp = levelUp->Blob(labelUp);
843  ssBlobUp = GetSSBlobFromGLBlob(blobUp);
844  if ( ssBlobUp != NULL) {
845  ssBlobUp->SetScaleMin(t_down);
846  bifurc->AddTopBlob(ssBlobUp);
847  ssBlobUp->SetBottomBifurcation(bifurc);
848  }
849  }
850  AddBifurcation(bifurc);
851  }
852  }
853  }
854  std::cout << "\t OK" << std::endl;
855 
856  }
857  }
858 
859 
860  //---------------------------------------------------------------------
861 
862  template<typename Geom, typename Text>
864  const std::string & statName )
865  {
866  ScaleSpaceBlob<Site> *ssBlob;
867  std::map<float, BlobMeasurements > statsM, statsSD;
868  FILE *fileStats;
869  float t;
870  float max_int, mean_int, max_cont, mean_cont, area, areamoy, areavar, tvalue, tv2;
871  BlobMeasurements *measurements;
872  GreyLevelBlob<Site> *glBlob, *glBlob1, *glBlob2;
873  int res;
874 
875  // Computing SSBlob measurements
876  // statM contains all the stats for normalizing GLBlobs measurements
877  // before integration across scales.
878 
879  std::cout << "Computing multiscale measurements for SSBlobs"
880  << std::endl;
881 
882  if (statName.length()>0)
883  {
884  if ((fileStats=fopen(statName.c_str(), "r"))==NULL)
885  {
886  std::cerr << "Problem : cannot open stat file " << statName
887  << std::endl;
888  exit(EXIT_FAILURE);
889  }
890  std::cout << "Found stat file, reading it" << std::endl;
891  while (!feof(fileStats))
892  {
893  res = fscanf(fileStats, "%f\n", &t); // !!! this defines the stat file format (simple)
894  res = fscanf(fileStats, "%f %f %f %f %f", &max_int, &mean_int, &max_cont, &mean_cont, &area);
895  measurements=new BlobMeasurements(max_int, mean_int, max_cont, mean_cont, area);
896  statsM.insert(std::pair<float, BlobMeasurements>(t, *measurements));
897  res = fscanf(fileStats, "%f %f %f %f %f", &max_int, &mean_int, &max_cont, &mean_cont, &area);
898  measurements=new BlobMeasurements(max_int, mean_int, max_cont, mean_cont, area);
899  statsSD.insert(std::pair<float, BlobMeasurements>(t, *measurements));
900  }
901  fclose(fileStats);
902  typename std::list<ScaleSpaceBlob<Site>*>::iterator itSSBlobs=blobList.begin();
903  std::cout << "processing GLblobs for normalisation" << std::endl;
904  for (; itSSBlobs!=blobList.end(); ++itSSBlobs)
905  {
906  ssBlob=*itSSBlobs;
907  typename std::list<GreyLevelBlob<Site>*>::iterator itglb;
908 
909  for (itglb=ssBlob->glBlobs.begin();itglb!=ssBlob->glBlobs.end();++itglb)
910  {
911  glBlob=*itglb;
912  t=glBlob->GetScale();
913  if (statsM.find(t)!=statsM.end())
914  {
915  glBlob->measurements.maxIntensity=(glBlob->measurements.maxIntensity-statsM[t].maxIntensity)/statsSD[t].maxIntensity;
916  glBlob->measurements.meanIntensity=(glBlob->measurements.meanIntensity-statsM[t].meanIntensity)/statsSD[t].meanIntensity;
917  glBlob->measurements.maxContrast=(glBlob->measurements.maxContrast-statsM[t].maxContrast)/statsSD[t].maxContrast;
918  glBlob->measurements.meanContrast=(glBlob->measurements.meanContrast-statsM[t].meanContrast)/statsSD[t].meanContrast;
919  glBlob->measurements.area=(glBlob->measurements.area-statsM[t].area)/statsSD[t].area;
920  if (glBlob->measurements.maxIntensity > 0)
922  else
923  glBlob->measurements.maxIntensity=exp(glBlob->measurements.maxIntensity);
924  if (glBlob->measurements.meanIntensity > 0)
926  else
928  if (glBlob->measurements.maxContrast > 0)
930  else
931  glBlob->measurements.maxContrast=exp(glBlob->measurements.maxContrast);
932  if (glBlob->measurements.meanContrast > 0)
934  else
935  glBlob->measurements.meanContrast=exp(glBlob->measurements.meanContrast);
936  if (glBlob->measurements.area > 0)
937  glBlob->measurements.area=1+glBlob->measurements.area;
938  else
939  glBlob->measurements.area=exp(glBlob->measurements.area);
940  }
941  else
942  {
943  std::cerr << "Problem during measurements normalisation !" << std::endl;
944  std::cerr << "Stat file " << statName << " incomplete, can't find scale " << t << std::endl;
945  exit(EXIT_FAILURE);
946  }
947  }
948  }
949  }
950  typename std::list< ScaleSpaceBlob<Site>* >::iterator itSSBlobs=blobList.begin();
951  std::cout << "Computing multi-scale measurements" << std::endl;
952 
953  for (; itSSBlobs!=blobList.end(); ++itSSBlobs)
954  {
955  ssBlob=*itSSBlobs;
956  typename std::list<GreyLevelBlob<Site>*>::iterator itGLBlobs=ssBlob->glBlobs.begin();
957  BlobMeasurements measure;
958  float t2,t1, dt;
959  max_int=0; mean_int=0; max_cont=0; mean_cont=0; area=0; areamoy=0; areavar=0; tvalue=0, tv2 = 0;
960  glBlob1=*itGLBlobs;
961  t1=glBlob1->GetScale();
962  dt=log(ssBlob->TopBifurcation()->tMid())-log(t1);
963  // DEBUG
964  if (dt< 0)
965  std::cout << "PROBLEME DE DT<0" << std::endl;
966  // FIN DEBUG
967  max_int += dt * (glBlob1->measurements.maxIntensity);
968  mean_int += dt * (glBlob1->measurements.meanIntensity);
969  max_cont += dt * (glBlob1->measurements.maxContrast);
970  mean_cont += dt * (glBlob1->measurements.meanContrast);
971  area += dt * (glBlob1->measurements.area);
972  tvalue += dt * (glBlob1->measurements.t);
973 
974 
975  ++itGLBlobs;
976  for (; itGLBlobs!=ssBlob->glBlobs.end(); ++itGLBlobs)
977  {
978  glBlob2=*itGLBlobs;
979  t2=glBlob2->GetScale();
980  dt=log(t1)-log(t2);
981  max_int += dt * (glBlob2->measurements.maxIntensity + glBlob1->measurements.maxIntensity)/2.0;
982  mean_int += dt * (glBlob2->measurements.meanIntensity + glBlob1->measurements.meanIntensity)/2.0;
983  max_cont += dt * (glBlob2->measurements.maxContrast + glBlob1->measurements.maxContrast)/2.0;
984  mean_cont += dt * (glBlob2->measurements.meanContrast + glBlob1->measurements.meanContrast)/2.0;
985  area += dt * (glBlob2->measurements.area + glBlob1->measurements.area)/2.0;
986  tvalue += dt * (glBlob2->measurements.t + glBlob1->measurements.t)/2.0;
987 
988 
989  glBlob1=glBlob2;
990  t1=t2;
991  }
992  t1=glBlob1->GetScale();
993  dt=log(t1) - log(ssBlob->BottomBifurcation()->tMid());
994  // DEBUG
995  if (dt< 0)
996  std::cout << "PROBLEME DE DT<0" << std::endl;
997  // FIN DEBUG
998  max_int += dt * (glBlob1->measurements.maxIntensity);
999  mean_int += dt * (glBlob1->measurements.meanIntensity);
1000  max_cont += dt * (glBlob1->measurements.maxContrast);
1001  mean_cont += dt * (glBlob1->measurements.meanContrast);
1002  area += dt * (glBlob1->measurements.area);
1003  tvalue += dt * (glBlob1->measurements.t);
1004 
1005 
1006  // Getting T-Value From Original Map...
1007 
1008  TexturedData<Geom, Text> ima=scaleSpace()->Scale(0.0)->Data();
1009 
1010  glBlob1=ssBlob->GlBlobRep();
1011  std::set<Site,ltstr_p3d<Site> > pixels;
1012  pixels=glBlob1->GetListePoints();
1013 
1014 
1015  typename std::set<Site, ltstr_p3d<Site> >::iterator itPix;
1016  float tvmax=-100.0;
1017  for (itPix=pixels.begin();itPix!=pixels.end();itPix++){
1018  if (ima.intensity(*itPix) > tvmax )
1019  tvmax = ima.intensity(*itPix);
1020  }
1021 
1022 
1023  measure=BlobMeasurements(max_int, mean_int, max_cont,mean_cont, area, tvalue, tvmax);
1024 
1025  areamoy +=area;
1026  areavar +=area*area;
1027  ssBlob->SetMeasurements(measure);
1028  }
1029  }
1030 
1031 
1032 }
1033 
1034 #endif
std::map< int, GreyLevelBlob< Site > * > BlobList()
Definition: scaleLevel.h:69
void ComputePrimalSketch(float tmin, float tMax, const std::string &statFile="", uint intersection_param=10)
Definition: primalSketch.h:180
BucketMap< Void > * mask(const BucketMap< Void > &src, const BucketMap< Void > &m, bool intersect=true)
std::list< ScaleSpaceBlob< Site > * > BlobSet()
Definition: primalSketch.h:112
PrimalSketch(const std::string &subject, ScaleSpace< Geom, Text > *scaleSpace, TexturedData< Geom, Text > *mask, int type)
Definition: primalSketch.h:91
std::list< ScaleSpaceBlob< T > * > TopBlobs()
std::set< T, ltstr_p3d< T > > & GetListePoints()
std::list< GreyLevelBlob< T > * > glBlobs
GreyLevelBlob< Site > * Blob(int label)
Definition: scaleLevel.h:73
std::list< Bifurcation< Site > * > BifurcationList()
Definition: primalSketch.h:110
ScaleSpace< Geom, Text > * _scaleSpace
Definition: primalSketch.h:67
GreyLevelBlob< T > * GlBlobRep()
BlobMeasurements measurements
SiteType< Geom >::type Site
Definition: primalSketch.h:62
std::string _subject
Definition: primalSketch.h:65
Bifurcation< T > * BottomBifurcation()
ScaleSpaceBlob< Site > * Blob(int label)
Definition: primalSketch.h:143
TexturedData< Geom, Text > * _mask
Definition: primalSketch.h:75
ScaleSpace< Geom, Text > * scaleSpace()
Definition: primalSketch.h:113
ScaleSpaceBlob< Site > * GetSSBlobFromGLBlob(GreyLevelBlob< Site > *glBlob)
Definition: primalSketch.h:156
void MatchScaleLevels(float t_up, float t_down, uint intersection_param)
Definition: primalSketch.h:483
std::list< ScaleSpaceBlob< T > * > BottomBlobs()
void setType(int t)
Definition: primalSketch.h:99
std::string Subject()
Definition: primalSketch.h:97
void setType(char type)
unsigned int uint
void AddGreyLevelBlob(GreyLevelBlob< T > *blob)
void SetTopBifurcation(Bifurcation< T > *bif)
void AddBifurcation(Bifurcation< Site > *bifurcation)
Definition: primalSketch.h:117
void SetMeasurements(BlobMeasurements meas)
void SetScaleMin(float t)
void ComputeBlobMeasurements(const std::string &statName="")
Definition: primalSketch.h:863
std::list< Bifurcation< Site > * > bifurcationList
Definition: primalSketch.h:68
void SetBottomBifurcation(Bifurcation< T > *bif)
std::list< ScaleSpaceBlob< Site > * > blobList
Definition: primalSketch.h:69
void SetScaleSpace(ScaleSpace< Geom, Text > *scaleSpace)
Definition: primalSketch.h:108
void AddTopBlob(ScaleSpaceBlob< T > *blob)
TexType< Text >::type Val
Definition: primalSketch.h:63
Bifurcation< T > * TopBifurcation()
void DetectBlobs(TexturedData< Geom, Text > *mask, char *stats=0)
PrimalSketch(const std::string &subject, ScaleSpace< Geom, Text > *scaleSpace, int type)
Definition: primalSketch.h:87
void AddBlob(ScaleSpaceBlob< Site > *blob)
Definition: primalSketch.h:115
void AddBottomBlob(ScaleSpaceBlob< T > *blob)
PrimalSketch(const std::string &subject, int type)
Definition: primalSketch.h:85
void SetScaleMax(float t)