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();
 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   // Check threshold
 191   if (m_UseLowerThreshold) {
 192     if (m_LowerThreshold > m_UpperThreshold) {
 193       this->SetLastError("ERROR: lower threshold cannot be greater than upper threshold.");
 194       return;
 195     }
 196   }
 197   // Threshold to get air
 198   typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
 199   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 200   binarizeFilter->SetInput(working_input);
 201   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 202   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 203   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 204   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 205   binarizeFilter->Update();
 206   working_image = binarizeFilter->GetOutput();
 207
 208   // Labelize and keep right labels
 209   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 210
 211   working_image = RemoveLabels<InternalImageType>
 212     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 213
 214   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 215     (working_image,
 216      GetBackgroundValue(),
 217      GetForegroundValue(),
 218      GetLabelizeParameters1()->GetFirstKeep(),
 219      GetLabelizeParameters1()->GetLastKeep(),
 220      GetLabelizeParameters1()->GetUseLastKeep());
 221
 222   // Set Air to BG
 223   working_input = SetBackground<ImageType, InternalImageType>
 224     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 225
 226   // End
 227   StopCurrentStep<ImageType>(working_input);
 228
 229   //--------------------------------------------------------------------
 230   //--------------------------------------------------------------------
 231   StartNewStepOrStop("Find the trachea");
 232   if (m_Seeds.size() == 0) { // try to find seed
 233     // Search seed (parameters = UpperThresholdForTrachea)
 234     static const unsigned int Dim = ImageType::ImageDimension;
 235     typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 236     typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 237     typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 238     sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
 239     sliceRegionSize[Dim-1]=5;
 240     sliceRegion.SetSize(sliceRegionSize);
 241     sliceRegion.SetIndex(sliceRegionIndex);
 242     typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 243     IteratorType it(working_input, sliceRegion);
 244     it.GoToBegin();
 245     while (!it.IsAtEnd()) {
 246       if(it.Get() < GetUpperThresholdForTrachea() ) {
 247         AddSeed(it.GetIndex());
 248       }
 249       ++it;
 250     }
 251   }
 252
 253   if (m_Seeds.size() != 0) {
 254     // Explosion controlled region growing
 255     typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
 256     typename ImageFilterType::Pointer f= ImageFilterType::New();
 257     f->SetInput(working_input);
 258     f->SetVerbose(false);
 259     f->SetLower(-2000);
 260     f->SetUpper(GetUpperThresholdForTrachea());
 261     f->SetMinimumLowerThreshold(-2000);
 262     f->SetMaximumUpperThreshold(0);
 263     f->SetAdaptLowerBorder(false);
 264     f->SetAdaptUpperBorder(true);
 265     f->SetMinimumSize(5000);
 266     f->SetReplaceValue(1);
 267     f->SetMultiplier(GetMultiplierForTrachea());
 268     f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 269     f->SetMinimumThresholdStepSize(1);
 270     for(unsigned int i=0; i<m_Seeds.size();i++) {
 271       f->AddSeed(m_Seeds[i]);
 272     }
 273     f->Update();
 274     trachea_tmp = f->GetOutput();
 275     // Set output
 276     StopCurrentStep<InternalImageType>(trachea_tmp);
 277   }
 278   else { // Trachea not found
 279     this->SetWarning("* WARNING * No seed found for trachea.");
 280     StopCurrentStep();
 281   }
 282
 283   //--------------------------------------------------------------------
 284   //--------------------------------------------------------------------
 285   StartNewStepOrStop("Extract the lung with Otsu filter");
 286   // Automated Otsu thresholding and relabeling
 287   typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
 288   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 289   otsuFilter->SetInput(working_input);
 290   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 291   otsuFilter->SetInsideValue(this->GetForegroundValue());
 292   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 293   otsuFilter->Update();
 294   working_image =  otsuFilter->GetOutput();
 295
 296   // Set output
 297   StopCurrentStep<InternalImageType>(working_image);
 298
 299   //--------------------------------------------------------------------
 300   //--------------------------------------------------------------------
 301   StartNewStepOrStop("Select labels");
 302   // Keep right labels
 303   working_image = LabelizeAndSelectLabels<InternalImageType>
 304     (working_image,
 305      GetBackgroundValue(),
 306      GetForegroundValue(),
 307      false,
 308      GetMinimalComponentSize(),
 309      GetLabelizeParameters2());
 310
 311   // Set output
 312   StopCurrentStep<InternalImageType>(working_image);
 313
 314   //--------------------------------------------------------------------
 315   //--------------------------------------------------------------------
 316   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 317     StartNewStepOrStop("Remove the trachea");
 318     // Set the trachea
 319     working_image = SetBackground<InternalImageType, InternalImageType>
 320       (working_image, trachea_tmp, 1, -1);
 321
 322    // Dilate the trachea
 323     static const unsigned int Dim = ImageType::ImageDimension;
 324     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 325     KernelType structuringElement;
 326     structuringElement.SetRadius(GetRadiusForTrachea());
 327     structuringElement.CreateStructuringElement();
 328     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 329     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 330     dilateFilter->SetBoundaryToForeground(false);
 331     dilateFilter->SetKernel(structuringElement);
 332     dilateFilter->SetBackgroundValue (1);
 333     dilateFilter->SetForegroundValue (-1);
 334     dilateFilter->SetInput (working_image);
 335     dilateFilter->Update();
 336     working_image = dilateFilter->GetOutput();
 337
 338     // Set trachea with dilatation
 339     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 340       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 341
 342     // Remove the trachea
 343     working_image = SetBackground<InternalImageType, InternalImageType>
 344       (working_image, working_image, -1, this->GetBackgroundValue());
 345
 346     // Label
 347     working_image = LabelizeAndSelectLabels<InternalImageType>
 348       (working_image,
 349        GetBackgroundValue(),
 350        GetForegroundValue(),
 351        false,
 352        GetMinimalComponentSize(),
 353        GetLabelizeParameters3());
 354
 355     // Set output
 356     StopCurrentStep<InternalImageType>(working_image);
 357   }
 358
 359   //--------------------------------------------------------------------
 360   //--------------------------------------------------------------------
 361   typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 362   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 363   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 364     StartNewStepOrStop("Croping trachea");
 365     cropFilter->SetInput(trachea_tmp);
 366     cropFilter->Update(); // Needed
 367     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 368     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 369     caster->SetInput(cropFilter->GetOutput());
 370     caster->Update();
 371     trachea = caster->GetOutput();
 372     StopCurrentStep<MaskImageType>(trachea);
 373   }
 374
 375   //--------------------------------------------------------------------
 376   //--------------------------------------------------------------------
 377   StartNewStepOrStop("Croping lung");
 378   typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
 379   cropFilter2->SetInput(working_image);
 380   cropFilter2->Update();
 381   working_image = cropFilter2->GetOutput();
 382   StopCurrentStep<InternalImageType>(working_image);
 383
 384   //--------------------------------------------------------------------
 385   //--------------------------------------------------------------------
 386   StartNewStepOrStop("Separate Left/Right lungs");
 387   // Initial label
 388   working_image = Labelize<InternalImageType>(working_image,
 389                                               GetBackgroundValue(),
 390                                               false,
 391                                               GetMinimalComponentSize());
 392
 393   // Count the labels
 394   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 395   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 396   statisticsImageFilter->SetInput(working_image);
 397   statisticsImageFilter->Update();
 398   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 399   working_image = statisticsImageFilter->GetOutput();
 400
 401   // Decompose the first label
 402   static const unsigned int Dim = ImageType::ImageDimension;
 403   if (initialNumberOfLabels<2) {
 404     // Structuring element radius
 405     typename ImageType::SizeType radius;
 406     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 407     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 408     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 409     decomposeAndReconstructFilter->SetInput(working_image);
 410     decomposeAndReconstructFilter->SetVerbose(false);
 411     decomposeAndReconstructFilter->SetRadius(radius);
 412     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 413     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 414     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 415     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 416     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 417     decomposeAndReconstructFilter->SetFullyConnected(true);
 418     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 419     decomposeAndReconstructFilter->Update();
 420     working_image = decomposeAndReconstructFilter->GetOutput();
 421   }
 422
 423   // Retain labels (lungs)
 424   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 425   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 426   thresholdFilter->SetInput(working_image);
 427   thresholdFilter->ThresholdAbove(2);
 428   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 429   thresholdFilter->Update();
 430   working_image = thresholdFilter->GetOutput();
 431   StopCurrentStep<InternalImageType> (working_image);
 432
 433   // Final Cast
 434   StartNewStepOrStop("Final cast");
 435   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 436   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 437   caster->SetInput(working_image);
 438   caster->Update();
 439   output = caster->GetOutput();
 440
 441   // Update output info
 442   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 443
 444   // Try to extract bifurcation in the trachea (bronchi)
 445   // STILL EXPERIMENTAL
 446   if (GetFindBronchialBifurcations()) {
 447     StartNewStepOrStop("Find bronchial bifurcations");
 448     // Step 1 : extract skeleton
 449     // Define the thinning filter
 450     typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
 451     typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
 452     thinningFilter->SetInput(trachea);
 453     thinningFilter->Update();
 454     typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
 455     writeImage<MaskImageType>(skeleton, "skeleton.mhd");
 456
 457     // Step 2 : tracking
 458     DD("tracking");
 459
 460     // Step 2.1 : find first point for tracking
 461     typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
 462     IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 463     it.GoToReverseBegin();
 464     while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 465       --it;
 466     }
 467     if (it.IsAtEnd()) {
 468       this->SetLastError("ERROR: first point in the skeleton not found ! Abort");
 469       return;
 470     }
 471     DD(skeleton->GetLargestPossibleRegion().GetIndex());
 472     typename MaskImageType::IndexType index = it.GetIndex();
 473     DD(index);
 474
 475     // Step 2.2 : initialize neighborhooditerator
 476     typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 477     typename NeighborhoodIteratorType::SizeType radius;
 478     radius.Fill(1);
 479     NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 480     DD(nit.GetSize());
 481     DD(nit.Size());
 482
 483     // Find first label number (must be different from BG and FG)
 484     typename MaskImageType::PixelType label = GetForegroundValue()+1;
 485     while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 486     DD(label);
 487
 488     // Track from the first point
 489     std::vector<BifurcationType> listOfBifurcations;
 490     TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
 491     DD("end track");
 492     DD(listOfBifurcations.size());
 493     writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
 494
 495     for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 496       typename MaskImageType::PointType p;
 497       skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
 498       DD(p);
 499     }
 500
 501   }
 502 }
 503 //--------------------------------------------------------------------
 504
 505
 506 //--------------------------------------------------------------------
 507 template <class ImageType, class MaskImageType>
 508 void
 509 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 510 GenerateData() {
 511   // Do not put some "startnewstep" here, because the object if
 512   // modified and the filter's pipeline it do two times. But it is
 513   // required to quit if MustStop was set before.
 514   if (GetMustStop()) return;
 515
 516   // If everything goes well, set the output
 517   this->GraftOutput(output); // not SetNthOutput
 518 }
 519 //--------------------------------------------------------------------
 520
 521
 522 //--------------------------------------------------------------------
 523 template <class ImageType, class MaskImageType>
 524 void
 525 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 526 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 527                    MaskImagePointer skeleton,
 528                    MaskImageIndexType index,
 529                    MaskImagePixelType label) {
 530   DD("TrackFromThisIndex");
 531   DD(index);
 532   DD((int)label);
 533   // Create NeighborhoodIterator
 534   typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 535   typename NeighborhoodIteratorType::SizeType radius;
 536   radius.Fill(1);
 537   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 538
 539   // Track
 540   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 541   bool stop = false;
 542   while (!stop) {
 543     nit.SetLocation(index);
 544     // DD((int)nit.GetCenterPixel());
 545     nit.SetCenterPixel(label);
 546     listOfTrackedPoint.clear();
 547     for(unsigned int i=0; i<nit.Size(); i++) {
 548       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 549         //          DD(nit.GetIndex(i));
 550         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 551           // DD(nit.GetIndex(i));
 552           listOfTrackedPoint.push_back(nit.GetIndex(i));
 553         }
 554       }
 555     }
 556     // DD(listOfTrackedPoint.size());
 557     if (listOfTrackedPoint.size() == 1) {
 558       index = listOfTrackedPoint[0];
 559     }
 560     else {
 561       if (listOfTrackedPoint.size() == 2) {
 562         BifurcationType bif(index, label, label+1, label+2);
 563         listOfBifurcations.push_back(bif);
 564         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 565         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 566       }
 567       else {
 568         if (listOfTrackedPoint.size() > 2) {
 569           std::cerr << "too much bifurcation points ... ?" << std::endl;
 570           exit(0);
 571         }
 572         // Else this it the end of the tracking
 573       }
 574       stop = true;
 575     }
 576   }
 577 }
 578 //--------------------------------------------------------------------
 579
 580
 581 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX