segmentation/clitkExtractLungFilter.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 CLITKEXTRACTLUNGSFILTER_TXX
  20 #define CLITKEXTRACTLUNGSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkImageCommon.h"
  24 #include "clitkSetBackgroundImageFilter.h"
  25 #include "clitkSegmentationUtils.h"
  26 #include "clitkAutoCropFilter.h"
  27
  28 // itk
  29 #include "itkBinaryThresholdImageFilter.h"
  30 #include "itkConnectedComponentImageFilter.h"
  31 #include "itkRelabelComponentImageFilter.h"
  32 #include "itkOtsuThresholdImageFilter.h"
  33 #include "itkBinaryThinningImageFilter3D.h"
  34 #include "itkImageIteratorWithIndex.h"
  35
  36
  37 //--------------------------------------------------------------------
  38 template <class ImageType, class MaskImageType>
  39 clitk::ExtractLungFilter<ImageType, MaskImageType>::
  40 ExtractLungFilter():
  41   clitk::FilterBase(),
  42   itk::ImageToImageFilter<ImageType, MaskImageType>()
  43 {
  44   SetNumberOfSteps(10);
  45   m_MaxSeedNumber = 500;
  46
  47   // Default global options
  48   this->SetNumberOfRequiredInputs(2);
  49   SetPatientMaskBackgroundValue(0);
  50   SetBackgroundValue(0); // Must be zero
  51   SetForegroundValue(1);
  52   SetMinimalComponentSize(100);
  53
  54   // Step 1 default values
  55   SetUpperThreshold(-300);
  56   SetLowerThreshold(-1000);
  57   UseLowerThresholdOff();
  58   LabelParamType * p1 = new LabelParamType;
  59   p1->SetFirstKeep(1);
  60   p1->UseLastKeepOff();
  61   p1->AddLabelToRemove(2);
  62   SetLabelizeParameters1(p1);
  63
  64   // Step 2 default values
  65   SetUpperThresholdForTrachea(-900);
  66   SetMultiplierForTrachea(5);
  67   SetThresholdStepSizeForTrachea(64);
  68
  69   // Step 3 default values
  70   SetNumberOfHistogramBins(500);
  71   LabelParamType * p2 = new LabelParamType;
  72   p2->SetFirstKeep(1);
  73   p2->UseLastKeepOff();
  74   // p->AddLabelToRemove(2); // No label to remove by default
  75   SetLabelizeParameters2(p2);
  76
  77   // Step 4 default values
  78   SetRadiusForTrachea(1);
  79   LabelParamType * p3 = new LabelParamType;
  80   p3->SetFirstKeep(1);
  81   p3->SetLastKeep(2);
  82   p3->UseLastKeepOff();
  83   SetLabelizeParameters3(p3);
  84
  85   // Step 5 : find bronchial bifurcations
  86   FindBronchialBifurcationsOn();
  87 }
  88 //--------------------------------------------------------------------
  89
  90
  91 //--------------------------------------------------------------------
  92 template <class ImageType, class MaskImageType>
  93 void
  94 clitk::ExtractLungFilter<ImageType, MaskImageType>::
  95 SetInput(const ImageType * image)
  96 {
  97   this->SetNthInput(0, const_cast<ImageType *>(image));
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 template <class ImageType, class MaskImageType>
 104 void
 105 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 106 SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg )
 107 {
 108   this->SetNthInput(1, const_cast<MaskImageType *>(image));
 109   SetPatientMaskBackgroundValue(bg);
 110 }
 111 //--------------------------------------------------------------------
 112
 113
 114 //--------------------------------------------------------------------
 115 template <class ImageType, class MaskImageType>
 116 void
 117 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 118 AddSeed(InternalIndexType s)
 119 {
 120   m_Seeds.push_back(s);
 121 }
 122 //--------------------------------------------------------------------
 123
 124
 125 //--------------------------------------------------------------------
 126 template <class ImageType, class MaskImageType>
 127 template<class ArgsInfoType>
 128 void
 129 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 130 SetArgsInfo(ArgsInfoType mArgsInfo)
 131 {
 132   SetVerboseOption_GGO(mArgsInfo);
 133   SetVerboseStep_GGO(mArgsInfo);
 134   SetWriteStep_GGO(mArgsInfo);
 135   SetVerboseWarningOff_GGO(mArgsInfo);
 136
 137   SetUpperThreshold_GGO(mArgsInfo);
 138   SetLowerThreshold_GGO(mArgsInfo);
 139   SetLabelizeParameters1_GGO(mArgsInfo);
 140   if (!mArgsInfo.remove1_given) {
 141     GetLabelizeParameters1()->AddLabelToRemove(2);
 142     if (GetVerboseOption()) VerboseOption("remove1", 2);
 143   }
 144
 145   SetUpperThresholdForTrachea_GGO(mArgsInfo);
 146   SetMultiplierForTrachea_GGO(mArgsInfo);
 147   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
 148   AddSeed_GGO(mArgsInfo);
 149
 150   SetMinimalComponentSize_GGO(mArgsInfo);
 151
 152   SetNumberOfHistogramBins_GGO(mArgsInfo);
 153   SetLabelizeParameters2_GGO(mArgsInfo);
 154
 155   SetRadiusForTrachea_GGO(mArgsInfo);
 156   SetLabelizeParameters3_GGO(mArgsInfo);
 157 }
 158 //--------------------------------------------------------------------
 159
 160
 161 //--------------------------------------------------------------------
 162 template <class ImageType, class MaskImageType>
 163 void
 164 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 165 GenerateOutputInformation()
 166 {
 167   Superclass::GenerateOutputInformation(); // Needed  ??
 168   this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
 169
 170   // Get input pointers
 171   patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
 172   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 173
 174   // Check image
 175   if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
 176     this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
 177     return;
 178   }
 179
 180   //--------------------------------------------------------------------
 181   //--------------------------------------------------------------------
 182   StartNewStepOrStop("Set background to initial image");
 183   working_input = SetBackground<ImageType, MaskImageType>
 184     (input, patient, GetPatientMaskBackgroundValue(), -1000);
 185   StopCurrentStep<ImageType>(working_input);
 186
 187   //--------------------------------------------------------------------
 188   //--------------------------------------------------------------------
 189   StartNewStepOrStop("Remove Air");
 190   // Threshold to get air
 191   typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
 192   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 193   binarizeFilter->SetInput(working_input);
 194   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 195   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 196   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 197   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 198   binarizeFilter->Update();
 199   working_image = binarizeFilter->GetOutput();
 200
 201   // Labelize and keep right labels
 202   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 203
 204   working_image = RemoveLabels<InternalImageType>
 205     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 206
 207   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 208     (working_image,
 209      GetBackgroundValue(),
 210      GetForegroundValue(),
 211      GetLabelizeParameters1()->GetFirstKeep(),
 212      GetLabelizeParameters1()->GetLastKeep(),
 213      GetLabelizeParameters1()->GetUseLastKeep());
 214
 215   // Set Air to BG
 216   working_input = SetBackground<ImageType, InternalImageType>
 217     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 218
 219   // End
 220   StopCurrentStep<ImageType>(working_input);
 221
 222   //--------------------------------------------------------------------
 223   //--------------------------------------------------------------------
 224   StartNewStepOrStop("Find the trachea");
 225   if (m_Seeds.size() == 0) { // try to find seed
 226     // Search seed (parameters = UpperThresholdForTrachea)
 227     static const unsigned int Dim = ImageType::ImageDimension;
 228     typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 229     typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 230     typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 231     sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
 232     sliceRegionSize[Dim-1]=5;
 233     sliceRegion.SetSize(sliceRegionSize);
 234     sliceRegion.SetIndex(sliceRegionIndex);
 235     typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 236     IteratorType it(working_input, sliceRegion);
 237     it.GoToBegin();
 238     while (!it.IsAtEnd()) {
 239       if(it.Get() < GetUpperThresholdForTrachea() ) {
 240         AddSeed(it.GetIndex());
 241       }
 242       ++it;
 243     }
 244   }
 245
 246   if (m_Seeds.size() != 0) {
 247     // Explosion controlled region growing
 248     typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
 249     typename ImageFilterType::Pointer f= ImageFilterType::New();
 250     f->SetInput(working_input);
 251     f->SetVerbose(false);
 252     f->SetLower(-2000);
 253     f->SetUpper(GetUpperThresholdForTrachea());
 254     f->SetMinimumLowerThreshold(-2000);
 255     f->SetMaximumUpperThreshold(0);
 256     f->SetAdaptLowerBorder(false);
 257     f->SetAdaptUpperBorder(true);
 258     f->SetMinimumSize(5000);
 259     f->SetReplaceValue(1);
 260     f->SetMultiplier(GetMultiplierForTrachea());
 261     f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 262     f->SetMinimumThresholdStepSize(1);
 263     for(unsigned int i=0; i<m_Seeds.size();i++) {
 264       f->AddSeed(m_Seeds[i]);
 265     }
 266     f->Update();
 267     trachea_tmp = f->GetOutput();
 268     // Set output
 269     StopCurrentStep<InternalImageType>(trachea_tmp);
 270   }
 271   else { // Trachea not found
 272     this->SetWarning("* WARNING * No seed found for trachea.");
 273     StopCurrentStep();
 274   }
 275
 276   //--------------------------------------------------------------------
 277   //--------------------------------------------------------------------
 278   StartNewStepOrStop("Extract the lung with Otsu filter");
 279   // Automated Otsu thresholding and relabeling
 280   typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
 281   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 282   otsuFilter->SetInput(working_input);
 283   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 284   otsuFilter->SetInsideValue(this->GetForegroundValue());
 285   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 286   otsuFilter->Update();
 287   working_image =  otsuFilter->GetOutput();
 288
 289   // Set output
 290   StopCurrentStep<InternalImageType>(working_image);
 291
 292   //--------------------------------------------------------------------
 293   //--------------------------------------------------------------------
 294   StartNewStepOrStop("Select labels");
 295   // Keep right labels
 296   working_image = LabelizeAndSelectLabels<InternalImageType>
 297     (working_image,
 298      GetBackgroundValue(),
 299      GetForegroundValue(),
 300      false,
 301      GetMinimalComponentSize(),
 302      GetLabelizeParameters2());
 303
 304   // Set output
 305   StopCurrentStep<InternalImageType>(working_image);
 306
 307   //--------------------------------------------------------------------
 308   //--------------------------------------------------------------------
 309   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 310     StartNewStepOrStop("Remove the trachea");
 311     // Set the trachea
 312     working_image = SetBackground<InternalImageType, InternalImageType>
 313       (working_image, trachea_tmp, 1, -1);
 314
 315    // Dilate the trachea
 316     static const unsigned int Dim = ImageType::ImageDimension;
 317     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 318     KernelType structuringElement;
 319     structuringElement.SetRadius(GetRadiusForTrachea());
 320     structuringElement.CreateStructuringElement();
 321     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 322     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 323     dilateFilter->SetBoundaryToForeground(false);
 324     dilateFilter->SetKernel(structuringElement);
 325     dilateFilter->SetBackgroundValue (1);
 326     dilateFilter->SetForegroundValue (-1);
 327     dilateFilter->SetInput (working_image);
 328     dilateFilter->Update();
 329     working_image = dilateFilter->GetOutput();
 330
 331     // Set trachea with dilatation
 332     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 333       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 334
 335     // Remove the trachea
 336     working_image = SetBackground<InternalImageType, InternalImageType>
 337       (working_image, working_image, -1, this->GetBackgroundValue());
 338
 339     // Label
 340     working_image = LabelizeAndSelectLabels<InternalImageType>
 341       (working_image,
 342        GetBackgroundValue(),
 343        GetForegroundValue(),
 344        false,
 345        GetMinimalComponentSize(),
 346        GetLabelizeParameters3());
 347
 348     // Set output
 349     StopCurrentStep<InternalImageType>(working_image);
 350   }
 351
 352   //--------------------------------------------------------------------
 353   //--------------------------------------------------------------------
 354   typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 355   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 356   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 357     StartNewStepOrStop("Croping trachea");
 358     cropFilter->SetInput(trachea_tmp);
 359     cropFilter->Update(); // Needed
 360     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 361     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 362     caster->SetInput(cropFilter->GetOutput());
 363     caster->Update();
 364     trachea = caster->GetOutput();
 365     StopCurrentStep<MaskImageType>(trachea);
 366   }
 367
 368   //--------------------------------------------------------------------
 369   //--------------------------------------------------------------------
 370   StartNewStepOrStop("Croping lung");
 371   typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
 372   cropFilter2->SetInput(working_image);
 373   cropFilter2->Update();
 374   working_image = cropFilter2->GetOutput();
 375   StopCurrentStep<InternalImageType>(working_image);
 376
 377   //--------------------------------------------------------------------
 378   //--------------------------------------------------------------------
 379   StartNewStepOrStop("Separate Left/Right lungs");
 380   // Initial label
 381   working_image = Labelize<InternalImageType>(working_image,
 382                                               GetBackgroundValue(),
 383                                               false,
 384                                               GetMinimalComponentSize());
 385
 386   // Count the labels
 387   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 388   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 389   statisticsImageFilter->SetInput(working_image);
 390   statisticsImageFilter->Update();
 391   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 392   working_image = statisticsImageFilter->GetOutput();
 393
 394   // Decompose the first label
 395   static const unsigned int Dim = ImageType::ImageDimension;
 396   if (initialNumberOfLabels<2) {
 397     // Structuring element radius
 398     typename ImageType::SizeType radius;
 399     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 400     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 401     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 402     decomposeAndReconstructFilter->SetInput(working_image);
 403     decomposeAndReconstructFilter->SetVerbose(false);
 404     decomposeAndReconstructFilter->SetRadius(radius);
 405     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 406     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 407     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 408     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 409     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 410     decomposeAndReconstructFilter->SetFullyConnected(true);
 411     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 412     decomposeAndReconstructFilter->Update();
 413     working_image = decomposeAndReconstructFilter->GetOutput();
 414   }
 415
 416   // Retain labels (lungs)
 417   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 418   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 419   thresholdFilter->SetInput(working_image);
 420   thresholdFilter->ThresholdAbove(2);
 421   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 422   thresholdFilter->Update();
 423   working_image = thresholdFilter->GetOutput();
 424   StopCurrentStep<InternalImageType> (working_image);
 425
 426   // Final Cast
 427   StartNewStepOrStop("Final cast");
 428   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 429   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 430   caster->SetInput(working_image);
 431   caster->Update();
 432   output = caster->GetOutput();
 433
 434   // Update output info
 435   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 436
 437   // Try to extract bifurcation in the trachea (bronchi)
 438   // STILL EXPERIMENTAL
 439   if (GetFindBronchialBifurcations()) {
 440     StartNewStepOrStop("Find bronchial bifurcations");
 441     // Step 1 : extract skeleton
 442     // Define the thinning filter
 443     typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
 444     typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
 445     thinningFilter->SetInput(trachea);
 446     thinningFilter->Update();
 447     typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
 448     writeImage<MaskImageType>(skeleton, "skeleton.mhd");
 449
 450     // Step 2 : tracking
 451     DD("tracking");
 452
 453     // Step 2.1 : find first point for tracking
 454     typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
 455     IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 456     it.GoToReverseBegin();
 457     while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 458       --it;
 459     }
 460     if (it.IsAtEnd()) {
 461       this->SetLastError("ERROR: first point in the skeleton not found ! Abord");
 462       return;
 463     }
 464     DD(skeleton->GetLargestPossibleRegion().GetIndex());
 465     typename MaskImageType::IndexType index = it.GetIndex();
 466     DD(index);
 467
 468     // Step 2.2 : initialize neighborhooditerator
 469     typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 470     typename NeighborhoodIteratorType::SizeType radius;
 471     radius.Fill(1);
 472     NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 473     DD(nit.GetSize());
 474     DD(nit.Size());
 475
 476     // Find first label number (must be different from BG and FG)
 477     typename MaskImageType::PixelType label = GetForegroundValue()+1;
 478     while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 479     DD(label);
 480
 481     // Track from the first point
 482     std::vector<BifurcationType> listOfBifurcations;
 483     TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
 484     DD("end track");
 485     DD(listOfBifurcations.size());
 486     writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
 487
 488     for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 489       typename MaskImageType::PointType p;
 490       skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
 491       DD(p);
 492     }
 493
 494   }
 495 }
 496 //--------------------------------------------------------------------
 497
 498
 499 //--------------------------------------------------------------------
 500 template <class ImageType, class MaskImageType>
 501 void
 502 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 503 GenerateData() {
 504   // Do not put some "startnewstep" here, because the object if
 505   // modified and the filter's pipeline it do two times. But it is
 506   // required to quit if MustStop was set before.
 507   if (GetMustStop()) return;
 508
 509   // If everything goes well, set the output
 510   this->GraftOutput(output); // not SetNthOutput
 511 }
 512 //--------------------------------------------------------------------
 513
 514
 515 //--------------------------------------------------------------------
 516 template <class ImageType, class MaskImageType>
 517 void
 518 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 519 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 520                    MaskImagePointer skeleton,
 521                    MaskImageIndexType index,
 522                    MaskImagePixelType label) {
 523   DD("TrackFromThisIndex");
 524   DD(index);
 525   DD((int)label);
 526   // Create NeighborhoodIterator
 527   typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 528   typename NeighborhoodIteratorType::SizeType radius;
 529   radius.Fill(1);
 530   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 531
 532   // Track
 533   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 534   bool stop = false;
 535   while (!stop) {
 536     nit.SetLocation(index);
 537     // DD((int)nit.GetCenterPixel());
 538     nit.SetCenterPixel(label);
 539     listOfTrackedPoint.clear();
 540     for(unsigned int i=0; i<nit.Size(); i++) {
 541       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 542         //          DD(nit.GetIndex(i));
 543         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 544           // DD(nit.GetIndex(i));
 545           listOfTrackedPoint.push_back(nit.GetIndex(i));
 546         }
 547       }
 548     }
 549     // DD(listOfTrackedPoint.size());
 550     if (listOfTrackedPoint.size() == 1) {
 551       index = listOfTrackedPoint[0];
 552     }
 553     else {
 554       if (listOfTrackedPoint.size() == 2) {
 555         BifurcationType bif(index, label, label+1, label+2);
 556         listOfBifurcations.push_back(bif);
 557         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 558         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 559       }
 560       else {
 561         if (listOfTrackedPoint.size() > 2) {
 562           std::cerr << "too much bifurcation points ... ?" << std::endl;
 563           exit(0);
 564         }
 565         // Else this it the end of the tracking
 566       }
 567       stop = true;
 568     }
 569   }
 570 }
 571 //--------------------------------------------------------------------
 572
 573
 574 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX