segmentation/clitkExtractAirwaysTreeInfoFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ===========================================================================**/
  18
  19 #ifndef CLITKEXTRACTAIRWAYSTREEINFOSFILTER_TXX
  20 #define CLITKEXTRACTAIRWAYSTREEINFOSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkImageCommon.h"
  24 #include "clitkSetBackgroundImageFilter.h"
  25 #include "clitkSegmentationUtils.h"
  26 #include "clitkAutoCropFilter.h"
  27 #include "clitkExtractSliceFilter.h"
  28
  29 // itk
  30 #include "itkBinaryThresholdImageFilter.h"
  31 #include "itkConnectedComponentImageFilter.h"
  32 #include "itkRelabelComponentImageFilter.h"
  33 #include "itkOtsuThresholdImageFilter.h"
  34 #include "itkBinaryThinningImageFilter3D.h"
  35 #include "itkImageIteratorWithIndex.h"
  36
  37
  38 //--------------------------------------------------------------------
  39 template <class ImageType>
  40 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  41 ExtractAirwaysTreeInfoFilter():
  42   clitk::FilterBase(),
  43   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  44   itk::ImageToImageFilter<ImageType, ImageType>()
  45 {
  46   // Default global options
  47   this->SetNumberOfRequiredInputs(1);
  48   SetBackgroundValue(0);
  49   SetForegroundValue(1);
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 template <class ImageType>
  56 void
  57 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  58 SetInput(const ImageType * image)
  59 {
  60   this->SetNthInput(0, const_cast<ImageType *>(image));
  61 }
  62 //--------------------------------------------------------------------
  63
  64
  65 //--------------------------------------------------------------------
  66 template <class ImageType>
  67 void
  68 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  69 GenerateOutputInformation()
  70 {
  71   Superclass::GenerateOutputInformation();
  72   //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
  73
  74   // Get input pointer
  75   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
  76
  77   // Try to load the DB
  78   try {
  79     LoadAFDB();
  80   } catch (clitk::ExceptionObject e) {
  81     // Do nothing if not found, it will be used anyway to write the result
  82   }
  83
  84
  85 }
  86 //--------------------------------------------------------------------
  87
  88
  89 //--------------------------------------------------------------------
  90 template <class ImageType>
  91 void
  92 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  93 GenerateData()
  94 {
  95
  96   // Crop just abov lungs ?
  97
  98   // read skeleton if already exist
  99   bool foundSkeleton = true;
 100   try {
 101     skeleton = GetAFDB()->template GetImage <ImageType>("skeleton");
 102   }
 103   catch (clitk::ExceptionObject e) {
 104     DD("did not found skeleton");
 105     foundSkeleton = false;
 106   }
 107
 108   // Extract skeleton
 109   if (!foundSkeleton) {
 110     StartNewStep("Thinning filter (skeleton)");
 111     typedef itk::BinaryThinningImageFilter3D<ImageType, ImageType> ThinningFilterType;
 112     typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
 113     thinningFilter->SetInput(input); // input = trachea
 114     thinningFilter->Update();
 115     skeleton = thinningFilter->GetOutput();
 116     StopCurrentStep<ImageType>(skeleton);
 117     writeImage<ImageType>(skeleton, "skeleton.mhd");
 118   }
 119
 120   // Find first point for tracking
 121   StartNewStep("Find first point for tracking");
 122   typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 123   IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 124   it.GoToReverseBegin();
 125   while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 126     --it;
 127   }
 128   if (it.IsAtEnd()) {
 129     clitkExceptionMacro("first point in the trachea skeleton not found.");
 130   }
 131   typename ImageType::IndexType index = it.GetIndex();
 132   DD(index);
 133
 134   if (0) {
 135     // Initialize neighborhooditerator
 136     typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 137     typename NeighborhoodIteratorType::SizeType radius;
 138     radius.Fill(1);
 139     NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 140   }
 141
 142   // Find a label number different from BG and FG
 143   typename ImageType::PixelType label = GetForegroundValue()+1;
 144   while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 145   DD((int)label);
 146
 147   /*
 148     Tracking ?
 149     Goal : find position of C, RUL, RML, RLL, LUL, LLL bronchus
 150     Carina : ok "easy", track, slice by slice until 2 path into different label
 151     -> follow at Right
 152        -
 153     -> follow at Left
 154    */
 155
 156   // NEW Track from the first point
 157   StartNewStep("Start tracking");
 158   FullTreeNodeType n;
 159   n.index = index;
 160   skeleton->TransformIndexToPhysicalPoint(n.index, n.point); // à mettre dans TrackFromThisIndex2
 161   mFullSkeletonTree.set_head(n);
 162   StructuralTreeNodeType sn;
 163   sn.index = index;
 164   skeleton->TransformIndexToPhysicalPoint(sn.index, sn.point);
 165   mStructuralSkeletonTree.set_head(sn);
 166   TrackFromThisIndex2(index, skeleton, label, mFullSkeletonTree.begin(), mStructuralSkeletonTree.begin());
 167   StopCurrentStep();
 168
 169   // Reput FG instead of label in the skeleton image
 170   skeleton = clitk::SetBackground<ImageType, ImageType>(skeleton, skeleton, label, GetForegroundValue(), true);
 171
 172   // Debug
 173   typename StructuralTreeType::iterator sit = mStructuralSkeletonTree.begin();
 174   while (sit != mStructuralSkeletonTree.end()) {
 175     DD(sit->point);
 176     ++sit;
 177   }
 178
 179   // compute weight : n longueurs à chaque bifurfaction.
 180   // parcours de FullTreeNodeType, à partir de leaf, remonter de proche en proche la distance eucl.
 181   // par post order
 182
 183   // Init
 184   typename FullTreeType::iterator fit = mFullSkeletonTree.begin();
 185   while (fit != mFullSkeletonTree.end()) {
 186     fit->weight = 0.0;
 187     ++fit;
 188   }
 189
 190   DD("compute weight");
 191   typename FullTreeType::post_order_iterator pit = mFullSkeletonTree.begin_post();
 192   while (pit != mFullSkeletonTree.end_post()) {
 193     //DD(pit->point);
 194     /*
 195     if (pit != mFullSkeletonTree.begin()) {
 196       typename FullTreeType::iterator parent = mFullSkeletonTree.parent(pit);
 197       double d = pit->point.EuclideanDistanceTo(parent->point);
 198       // DD(parent->weight);
 199       //DD(d);
 200       parent->weight += d;
 201       if (pit.number_of_children() > 1) {
 202         DD(pit.number_of_children());
 203         DD(pit->point);
 204         DD(pit->weight);
 205         DD(parent->weight);
 206         DD(mFullSkeletonTree.child(pit, 0)->weight);
 207         DD(mFullSkeletonTree.child(pit, 1)->weight);
 208       }
 209     }
 210     */
 211     double previous = pit->weight;
 212     for(uint i=0; i<pit.number_of_children(); i++) {
 213       double d = pit->point.EuclideanDistanceTo(mFullSkeletonTree.child(pit, i)->point);
 214       //      pit->weight = pit->weight + mFullSkeletonTree.child(pit, i)->weight + d;
 215       if (i==0)
 216         pit->weight = pit->weight + mFullSkeletonTree.child(pit, i)->weight + d;
 217       if (i>0) {
 218         DD(pit.number_of_children());
 219         DD(pit->point);
 220         DD(pit->weight);
 221         DD(mFullSkeletonTree.child(pit, 0)->weight);
 222         DD(mFullSkeletonTree.child(pit, 1)->weight);
 223         pit->weight = std::max(pit->weight, previous+mFullSkeletonTree.child(pit, i)->weight + d);
 224       }
 225     }
 226     ++pit;
 227   }
 228
 229   DD("end");
 230   fit = mFullSkeletonTree.begin();
 231   while (fit != mFullSkeletonTree.end()) {
 232     std::cout << "p = " << fit->point
 233               << " " << fit->weight
 234               << " child= " << fit.number_of_children()
 235               << " -> ";
 236     for(uint i=0; i<fit.number_of_children(); i++) {
 237       std::cout << " " << mFullSkeletonTree.child(fit, i)->weight;
 238     }
 239     std::cout << std::endl;
 240     ++fit;
 241   }
 242
 243   /* Selection criteria ?
 244      - at least 2 child
 245      - from top to bottom
 246      - equilibr ?
 247    */
 248   DD("=========================================");
 249   fit = mFullSkeletonTree.begin();
 250   while (fit != mFullSkeletonTree.end()) {
 251     if (fit.number_of_children() > 1) {
 252       std::cout << "ppp = " << fit->point
 253                 << " " << fit->weight << std::endl;
 254       for(uint i=1; i<fit.number_of_children(); i++) {
 255         double w1 = mFullSkeletonTree.child(fit, i-1)->weight;
 256         double w2 = mFullSkeletonTree.child(fit, i)->weight;
 257         // if (w1 <= 0.1) break;
 258         // if (w2 <= 0.1) break;
 259         DD(w1);DD(w2);
 260         if (w1>w2) {
 261           double sw=w1;
 262           w1=w2; w2=sw;
 263         }
 264         DD(w1/w2);
 265         if (w1/w2 > 0.7) {
 266           DD("KEEP IT");
 267         }
 268       }
 269     }
 270     ++fit;
 271   }
 272
 273
 274   if (0) {
 275     // Track from the first point
 276     StartNewStep("Start tracking OLD");
 277     std::vector<BifurcationType> listOfBifurcations;
 278     m_SkeletonTree.set_head(index);
 279     TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
 280     DD("end track");
 281
 282   // Convert index into physical point coordinates
 283     for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 284       skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index,
 285                                               listOfBifurcations[i].point);
 286     }
 287
 288   // Search for the first slice that separate the bronchus
 289   // (carina). Labelize slice by slice, stop when the two points of
 290   // the skeleton ar not in the same connected component
 291   StartNewStep("Search for carina position");
 292   typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 293   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 294   extractSliceFilter->SetInput(input);
 295   extractSliceFilter->SetDirection(2);
 296   extractSliceFilter->Update();
 297   typedef typename ExtractSliceFilterType::SliceType SliceType;
 298   std::vector<typename SliceType::Pointer> mInputSlices;
 299   extractSliceFilter->GetOutputSlices(mInputSlices);
 300   DD(mInputSlices.size());
 301
 302
 303   DD("REDO !!!!!!!!!!!!");
 304   /**
 305      => chercher la bif qui a les plus important sous-arbres
 306      - iterate sur m_SkeletonTree depth first, remonter au parent, add value = distance
 307
 308      NONNNNN m_SkeletonTree n'est pas la structure mais l'ensemble avec ts le skeleton
 309
 310   TreeIterator itt = m_SkeletonTree.begin();
 311   while (itt != m_SkeletonTree.end()) {
 312     itt->weight = 0.0;
 313     ++itt;
 314   }
 315
 316   typename tree<NodeType>::post_order_iterator pit = m_SkeletonTree.begin_post();
 317   while (pit != m_SkeletonTree.end_post()) {
 318     typename tree<NodeType>::iterator parent = m_SkeletonTree.parent(pit);
 319     double d = pit->point.EuclideanDistanceTo(parent);
 320     DD(parent.weight);
 321     DD(d);
 322     parent.weight += d;
 323     ++it;
 324   }
 325
 326   itt = m_SkeletonTree.begin();
 327   while (itt != m_SkeletonTree.end()) {
 328     DD(itt->weight);
 329     ++itt;
 330   }
 331    **/
 332
 333
 334   // search for carina
 335   bool stop = false;
 336   int slice_index = listOfBifurcations[0].index[2]; // first slice from carina in skeleton
 337   int i=0;
 338   TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 0);
 339   TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 1);
 340   DD(firstIter.number_of_children());
 341   DD(secondIter.number_of_children());
 342   typename SliceType::IndexType in1;
 343   typename SliceType::IndexType in2;
 344   while (!stop) {
 345     //  Labelize the current slice
 346     typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
 347                                                            GetBackgroundValue(),
 348                                                            true,
 349                                                            0); // min component size=0
 350     DD(*firstIter);
 351     DD(*secondIter);
 352     // Check the value of the two skeleton points;
 353     in1[0] = (*firstIter)[0];
 354     in1[1] = (*firstIter)[1];
 355     typename SliceType::PixelType v1 = temp->GetPixel(in1);
 356     in2[0] = (*secondIter)[0];
 357     in2[1] = (*secondIter)[1];
 358     typename SliceType::PixelType v2 = temp->GetPixel(in2);
 359
 360     // Check the label value of the two points
 361     DD(slice_index);
 362     if (v1 != v2) {
 363       stop = true; // We found it !
 364     }
 365     else {
 366       // Check error
 367       if (slice_index == (int)(mInputSlices.size()-1)) {
 368         clitkExceptionMacro("Error while searching for carina, the two skeleton points are always in the same CC ... ???");
 369       }
 370       // Iterate
 371       i++;
 372       --slice_index;
 373       ++firstIter;
 374       ++secondIter;
 375     }
 376   }
 377   ImageIndexType carina_index; // middle position in X/Y
 378   carina_index[0] = lrint(in2[0] + in1[0])/2.0;
 379   carina_index[1] = lrint(in2[1] + in1[1])/2.0;
 380   carina_index[2] = slice_index;
 381   // Get physical coordinates
 382   ImagePointType carina_position;
 383   skeleton->TransformIndexToPhysicalPoint(carina_index,
 384                                           carina_position);
 385
 386   // Set and save Carina position
 387   if (GetVerboseStepFlag()) {
 388     std::cout << "\t Found carina at " << carina_position << " mm" << std::endl;
 389   }
 390   GetAFDB()->SetPoint3D("carina", carina_position);
 391
 392   // Write bifurcation (debug)
 393   for(uint i=0; i<listOfBifurcations.size(); i++) {
 394     GetAFDB()->SetPoint3D("Bif"+toString(i), listOfBifurcations[i].point);
 395   }
 396   }
 397
 398   // Set the output (skeleton);
 399   this->GraftOutput(skeleton); // not SetNthOutput
 400 }
 401 //--------------------------------------------------------------------
 402
 403
 404 //--------------------------------------------------------------------
 405 /**
 406    Start from the pixel "index" in the image "skeleton" and track the
 407    path from neighbor pixels. Put the tracked path in the tree at the
 408    currentNode position. Label is used to mark the FG pixel as already
 409    visited. Progress recursively when several neighbors are found.
 410  **/
 411 template <class ImageType>
 412 void
 413 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
 414 TrackFromThisIndex2(ImageIndexType index, ImagePointer skeleton,
 415                     ImagePixelType label,
 416                     typename FullTreeType::iterator currentNode,
 417                     typename StructuralTreeType::iterator currentSNode)
 418 {
 419   // Create NeighborhoodIterator
 420   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 421   typename NeighborhoodIteratorType::SizeType radius;
 422   radius.Fill(1);
 423   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 424
 425   // Track
 426   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 427   bool stop = false;
 428   while (!stop) {
 429     nit.SetLocation(index);
 430     nit.SetCenterPixel(label);
 431     listOfTrackedPoint.clear();
 432     for(unsigned int i=0; i<nit.Size(); i++) {
 433       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 434         if (nit.GetPixel(i) == GetForegroundValue()) {
 435           // if pixel value is FG, we keep this point
 436           listOfTrackedPoint.push_back(nit.GetIndex(i));
 437         }
 438       }
 439     }
 440     // If only neighbor is found, we keep the point and continue
 441     if (listOfTrackedPoint.size() == 1) {
 442       FullTreeNodeType n;
 443       index = n.index = listOfTrackedPoint[0];
 444       skeleton->TransformIndexToPhysicalPoint(n.index, n.point);
 445       currentNode = mFullSkeletonTree.append_child(currentNode, n);
 446       skeleton->SetPixel(n.index, label); // change label in skeleton image
 447     }
 448     else {
 449       if (listOfTrackedPoint.size() >= 2) {
 450         for(uint i=0; i<listOfTrackedPoint.size(); i++) {
 451           FullTreeNodeType n;
 452           n.index = listOfTrackedPoint[i];
 453           skeleton->TransformIndexToPhysicalPoint(n.index, n.point);
 454           typename FullTreeType::iterator node = mFullSkeletonTree.append_child(currentNode, n);
 455           StructuralTreeNodeType sn;
 456           sn.index = listOfTrackedPoint[i];
 457           skeleton->TransformIndexToPhysicalPoint(sn.index, sn.point);
 458           typename StructuralTreeType::iterator snode = mStructuralSkeletonTree.append_child(currentSNode, sn);
 459           TrackFromThisIndex2(listOfTrackedPoint[i], skeleton, label, node, snode);
 460         }
 461       }
 462       stop = true; // this it the end of the tracking
 463     } // end else
 464   } // end while stop
 465 }
 466 //--------------------------------------------------------------------
 467
 468 //--------------------------------------------------------------------
 469 template <class ImageType>
 470 void
 471 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
 472 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 473                    ImagePointer skeleton,
 474                    ImageIndexType index,
 475                    ImagePixelType label,
 476                    TreeIterator currentNode)
 477 {
 478   // Create NeighborhoodIterator
 479   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 480   typename NeighborhoodIteratorType::SizeType radius;
 481   radius.Fill(1);
 482   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 483
 484   // Track
 485   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 486   bool stop = false;
 487   while (!stop) {
 488     nit.SetLocation(index);
 489     nit.SetCenterPixel(label);
 490     listOfTrackedPoint.clear();
 491     for(unsigned int i=0; i<nit.Size(); i++) {
 492       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 493         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 494           listOfTrackedPoint.push_back(nit.GetIndex(i));
 495         }
 496       }
 497     }
 498     if (listOfTrackedPoint.size() == 1) {
 499       // Add this point to the current path
 500       currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 501       index = listOfTrackedPoint[0];
 502       skeleton->SetPixel(index, label); // change label in skeleton image
 503     }
 504     else {
 505       if (listOfTrackedPoint.size() == 2) {
 506         // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
 507         // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
 508         DD("BifurcationType");
 509         DD(listOfTrackedPoint[0]);
 510         DD(listOfTrackedPoint[1]);
 511         BifurcationType bif(index, label, label+1, label+2);
 512         bif.treeIter = currentNode;
 513         listOfBifurcations.push_back(bif);
 514         TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 515         TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
 516         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
 517         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
 518       }
 519       else {
 520         //DD(listOfTrackedPoint.size());
 521         if (listOfTrackedPoint.size() > 2) {
 522           //clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
 523           stop = true; // this it the end of the tracking
 524         }
 525         // Else this it the end of the tracking
 526       }
 527       stop = true;
 528     }
 529   }
 530 }
 531 //--------------------------------------------------------------------
 532
 533 /*TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 534   ImagePointer skeleton,
 535   ImageIndexType index,
 536   ImagePixelType label)
 537   {
 538   DD("TrackFromThisIndex");
 539   DD(index);
 540   DD((int)label);
 541   // Create NeighborhoodIterator
 542   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 543   typename NeighborhoodIteratorType::SizeType radius;
 544   radius.Fill(1);
 545   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 546
 547   // Track
 548   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 549   bool stop = false;
 550   while (!stop) {
 551   nit.SetLocation(index);
 552   // DD((int)nit.GetCenterPixel());
 553   nit.SetCenterPixel(label);
 554   listOfTrackedPoint.clear();
 555   for(unsigned int i=0; i<nit.Size(); i++) {
 556   if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 557   //          DD(nit.GetIndex(i));
 558   if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 559   // DD(nit.GetIndex(i));
 560   listOfTrackedPoint.push_back(nit.GetIndex(i));
 561   }
 562   }
 563   }
 564   // DD(listOfTrackedPoint.size());
 565   if (listOfTrackedPoint.size() == 1) {
 566   index = listOfTrackedPoint[0];
 567   }
 568   else {
 569   if (listOfTrackedPoint.size() == 2) {
 570   BifurcationType bif(index, label, label+1, label+2);
 571   listOfBifurcations.push_back(bif);
 572   TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 573   TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 574   }
 575   else {
 576   if (listOfTrackedPoint.size() > 2) {
 577   std::cerr << "too much bifurcation points ... ?" << std::endl;
 578   exit(0);
 579   }
 580   // Else this it the end of the tracking
 581   }
 582   stop = true;
 583   }
 584   }
 585   }
 586 */
 587 //--------------------------------------------------------------------
 588
 589
 590 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX