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