X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=3de027cc3b8210d0dbbac0e5d2a24d1a4aa02cfa;hb=22681675b0ed087fd02fbd94968d5ea5251f5408;hp=2c7be4c3546fb80629a7b25af8667981582e27e5;hpb=e008d74b0ecdc4ca2eaae8c429901a78f9ef5c31;p=clitk.git

diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx
index 2c7be4c..3de027c 100644
--- a/segmentation/clitkExtractLungFilter.txx
+++ b/segmentation/clitkExtractLungFilter.txx
@@ -22,7 +22,7 @@
 // clitk
 #include "clitkImageCommon.h"
 #include "clitkSetBackgroundImageFilter.h"
-#include "clitkSegmentationFunctions.h"
+#include "clitkSegmentationUtils.h"
 #include "clitkAutoCropFilter.h"
 
 // itk
@@ -30,19 +30,27 @@
 #include "itkConnectedComponentImageFilter.h"
 #include "itkRelabelComponentImageFilter.h"
 #include "itkOtsuThresholdImageFilter.h"
+#include "itkBinaryThinningImageFilter3D.h"
+#include "itkImageIteratorWithIndex.h"
+#include "itkBinaryMorphologicalOpeningImageFilter.h"
+#include "itkBinaryMorphologicalClosingImageFilter.h"
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+template <class ImageType, class MaskImageType>
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
 ExtractLungFilter():
   clitk::FilterBase(),
-  itk::ImageToImageFilter<TInputImageType, TMaskImageType>()
+  itk::ImageToImageFilter<ImageType, MaskImageType>()
 {
+  SetNumberOfSteps(10);
+  m_MaxSeedNumber = 500;
+
   // Default global options
   this->SetNumberOfRequiredInputs(2);
   SetPatientMaskBackgroundValue(0);
   SetBackgroundValue(0); // Must be zero
   SetForegroundValue(1);
+  SetMinimalComponentSize(100);
 
   // Step 1 default values
   SetUpperThreshold(-300);
@@ -58,6 +66,7 @@ ExtractLungFilter():
   SetUpperThresholdForTrachea(-900);
   SetMultiplierForTrachea(5);
   SetThresholdStepSizeForTrachea(64);
+  SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
   
   // Step 3 default values
   SetNumberOfHistogramBins(500);
@@ -74,37 +83,41 @@ ExtractLungFilter():
   p3->SetLastKeep(2);
   p3->UseLastKeepOff();
   SetLabelizeParameters3(p3);
+  
+  // Step 5
+  FinalOpenCloseOff();
+  SetFinalOpenCloseRadius(1);
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
 void 
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
-SetInput(const TInputImageType * image) 
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SetInput(const ImageType * image) 
 {
-  this->SetNthInput(0, const_cast<TInputImageType *>(image));
+  this->SetNthInput(0, const_cast<ImageType *>(image));
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
 void 
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
-SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg ) 
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg ) 
 {
-  this->SetNthInput(1, const_cast<TMaskImageType *>(image));
+  this->SetNthInput(1, const_cast<MaskImageType *>(image));
   SetPatientMaskBackgroundValue(bg);
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
 void 
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
 AddSeed(InternalIndexType s) 
 { 
   m_Seeds.push_back(s);
@@ -113,10 +126,10 @@ AddSeed(InternalIndexType s)
 
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
 template<class ArgsInfoType>
 void 
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
 SetArgsInfo(ArgsInfoType mArgsInfo)
 {
   SetVerboseOption_GGO(mArgsInfo);
@@ -126,6 +139,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
 
   SetUpperThreshold_GGO(mArgsInfo);
   SetLowerThreshold_GGO(mArgsInfo);
+  SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
   SetLabelizeParameters1_GGO(mArgsInfo);
   if (!mArgsInfo.remove1_given) {
     GetLabelizeParameters1()->AddLabelToRemove(2); 
@@ -137,47 +151,55 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
   AddSeed_GGO(mArgsInfo);
 
+  SetMinimalComponentSize_GGO(mArgsInfo);
+  
   SetNumberOfHistogramBins_GGO(mArgsInfo);
   SetLabelizeParameters2_GGO(mArgsInfo);
 
   SetRadiusForTrachea_GGO(mArgsInfo);
   SetLabelizeParameters3_GGO(mArgsInfo);
+  
+  SetFinalOpenCloseRadius_GGO(mArgsInfo);
+  SetFinalOpenClose_GGO(mArgsInfo);
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
 void 
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
 GenerateOutputInformation() 
 { 
-  input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
   Superclass::GenerateOutputInformation();
-//   MaskImagePointer output = this->GetOutput(0);
 
   // Get input pointers
-  input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
-  patient = dynamic_cast<const TMaskImageType*>(itk::ProcessObject::GetInput(1));
+  patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
+  input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 
   // Check image
-  if (!HasSameSizeAndSpacing<TInputImageType, TMaskImageType>(input, patient)) {
-    this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
-    return;
+  if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
+    clitkExceptionMacro("the 'input' and 'patient' masks must have the same size & spacing.");
   }
   
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Set background to initial image");
-  working_input = SetBackground<TInputImageType, TMaskImageType>
+  working_input = SetBackground<ImageType, MaskImageType>
     (input, patient, GetPatientMaskBackgroundValue(), -1000);
-  StopCurrentStep<InputImageType>(working_input);
+  StopCurrentStep<ImageType>(working_input);
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Remove Air");
+  // Check threshold
+  if (m_UseLowerThreshold) {
+    if (m_LowerThreshold > m_UpperThreshold) {
+      clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
+    }
+  }
   // Threshold to get air
-  typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType; 
+  typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType; 
   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
   binarizeFilter->SetInput(working_input);
   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
@@ -189,8 +211,10 @@ GenerateOutputInformation()
   
   // Labelize and keep right labels
   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
+
   working_image = RemoveLabels<InternalImageType>
     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
+
   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
     (working_image, 
      GetBackgroundValue(), 
@@ -200,75 +224,22 @@ GenerateOutputInformation()
      GetLabelizeParameters1()->GetUseLastKeep());
  
   // Set Air to BG
-  working_input = SetBackground<TInputImageType, InternalImageType>
+  working_input = SetBackground<ImageType, InternalImageType>
     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 
   // End
-  StopCurrentStep<InputImageType>(working_input);
+  StopCurrentStep<ImageType>(working_input);
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
-  StartNewStep("Find the trachea");
-  if (m_Seeds.size() == 0) { // try to find seed
-    // Search seed (parameters = UpperThresholdForTrachea)
-    static const unsigned int Dim = InputImageType::ImageDimension;
-    typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
-    typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
-    typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
-    sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
-    sliceRegionSize[Dim-1]=5;
-    sliceRegion.SetSize(sliceRegionSize);
-    sliceRegion.SetIndex(sliceRegionIndex);
-    
-    typedef  itk::ImageRegionConstIterator<InputImageType> IteratorType;
-    IteratorType it(working_input, sliceRegion);
-    it.GoToBegin();
-    while (!it.IsAtEnd()) {
-      if(it.Get() < GetUpperThresholdForTrachea() ) {
-        AddSeed(it.GetIndex());
-      }
-      ++it;
-    }
-  }
-  
-  if (m_Seeds.size() != 0) {
-    // Explosion controlled region growing
-    typedef clitk::ExplosionControlledThresholdConnectedImageFilter<InputImageType, InternalImageType> ImageFilterType;
-    typename ImageFilterType::Pointer f= ImageFilterType::New();
-    f->SetInput(working_input);
-    f->SetVerbose(false);
-    f->SetLower(-2000);
-    f->SetUpper(GetUpperThresholdForTrachea());
-    f->SetMinimumLowerThreshold(-2000);
-    f->SetMaximumUpperThreshold(0);
-    f->SetAdaptLowerBorder(false);
-    f->SetAdaptUpperBorder(true);
-    f->SetMinimumSize(5000); 
-    f->SetReplaceValue(1);
-    f->SetMultiplier(GetMultiplierForTrachea());
-    f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
-    f->SetMinimumThresholdStepSize(1);
-    for(unsigned int i=0; i<m_Seeds.size();i++) {
-      //    std::cout<<"Adding seed " <<m_Seeds[i]<<"..."<<std::endl;
-      f->AddSeed(m_Seeds[i]);
-    }  
-    f->Update();
-    trachea_tmp = f->GetOutput();
-    // Set output
-    StopCurrentStep<InternalImageType>(trachea_tmp);
-
-  }
-  else { // Trachea not found
-    this->SetWarning("* WARNING * No seed found for trachea.");
-    // Set output
-    StopCurrentStep();
-  }
+  StartNewStep("Search for the trachea");
+  SearchForTrachea();
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Extract the lung with Otsu filter");
   // Automated Otsu thresholding and relabeling
-  typedef itk::OtsuThresholdImageFilter<InputImageType,InternalImageType> OtsuThresholdImageFilterType;
+  typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
   otsuFilter->SetInput(working_input);
   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
@@ -303,8 +274,8 @@ GenerateOutputInformation()
     working_image = SetBackground<InternalImageType, InternalImageType>
       (working_image, trachea_tmp, 1, -1);
   
-   // Dilate the trachea 
-    static const unsigned int Dim = InputImageType::ImageDimension;
+    // Dilate the trachea 
+    static const unsigned int Dim = ImageType::ImageDimension;
     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
     KernelType structuringElement;
     structuringElement.SetRadius(GetRadiusForTrachea());
@@ -339,7 +310,6 @@ GenerateOutputInformation()
     // Set output
     StopCurrentStep<InternalImageType>(working_image);
   }
-  
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
@@ -360,12 +330,43 @@ GenerateOutputInformation()
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Croping lung");
-  cropFilter = CropFilterType::New(); // Needed to reset pipeline
-  cropFilter->SetInput(working_image);
-  cropFilter->Update();   
-  working_image = cropFilter->GetOutput();
+  typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
+  cropFilter2->SetInput(working_image);
+  cropFilter2->Update();   
+  working_image = cropFilter2->GetOutput();
   StopCurrentStep<InternalImageType>(working_image);
 
+  //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // Final OpenClose
+  if (GetFinalOpenClose()) {
+    StartNewStep("Open/Close"); 
+
+    // Structuring element
+    typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
+    KernelType structuringElement;
+    structuringElement.SetRadius(GetFinalOpenCloseRadius());
+    structuringElement.CreateStructuringElement();
+	
+    // Open
+    typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType, KernelType> OpenFilterType;
+    typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
+    openFilter->SetInput(working_image);
+    openFilter->SetBackgroundValue(GetBackgroundValue());
+    openFilter->SetForegroundValue(GetForegroundValue());
+    openFilter->SetKernel(structuringElement);
+	
+    // Close
+    typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType, KernelType> CloseFilterType;
+    typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
+    closeFilter->SetInput(openFilter->GetOutput());
+    closeFilter->SetSafeBorder(true);
+    closeFilter->SetForegroundValue(GetForegroundValue());
+    closeFilter->SetKernel(structuringElement);
+    closeFilter->Update();
+    working_image = closeFilter->GetOutput();
+  }
+
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Separate Left/Right lungs");
@@ -384,10 +385,10 @@ GenerateOutputInformation()
   working_image = statisticsImageFilter->GetOutput();	
  
   // Decompose the first label
-  static const unsigned int Dim = InputImageType::ImageDimension;
+  static const unsigned int Dim = ImageType::ImageDimension;
   if (initialNumberOfLabels<2) {
     // Structuring element radius
-    typename InputImageType::SizeType radius;
+    typename ImageType::SizeType radius;
     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
@@ -405,7 +406,7 @@ GenerateOutputInformation()
     working_image = decomposeAndReconstructFilter->GetOutput();      
   }
 
-  // Retain labels (lungs)
+  // Retain labels ('1' is largset lung, so right. '2' is left)
   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
   thresholdFilter->SetInput(working_image);
@@ -414,26 +415,198 @@ GenerateOutputInformation()
   thresholdFilter->Update();
   working_image = thresholdFilter->GetOutput();
   StopCurrentStep<InternalImageType> (working_image);
+  
+  // Final Cast 
+  StartNewStep("Cast the lung mask"); 
+  typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
+  typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
+  caster->SetInput(working_image);
+  caster->Update();
+  output = caster->GetOutput();
+
+  // Update output info
+  this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
+
+
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
 void 
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
-GenerateData() {
-  // Final Cast 
-  typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
-  typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
-  caster->SetInput(working_image);
-  caster->Update();
-  // Set output
-  //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise
-  this->GraftOutput(caster->GetOutput());
-  return;
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+GenerateData() 
+{
+  // Set the output
+  this->GraftOutput(output); // not SetNthOutput
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+bool 
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SearchForTracheaSeed(int skip)
+{
+  if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)    
+    // Restart the search until a seed is found, skipping more and more slices
+    bool stop = false;
+    while (!stop) {
+      // Search seed (parameters = UpperThresholdForTrachea)
+      static const unsigned int Dim = ImageType::ImageDimension;
+      typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
+      typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
+      typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
+      sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
+      sliceRegionSize[Dim-1]=5;
+      sliceRegion.SetSize(sliceRegionSize);
+      sliceRegion.SetIndex(sliceRegionIndex);
+      typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
+      IteratorType it(working_input, sliceRegion);
+      it.GoToBegin();
+      while (!it.IsAtEnd()) {
+        if(it.Get() < GetUpperThresholdForTrachea() ) {
+          AddSeed(it.GetIndex());
+          // DD(it.GetIndex());
+        }
+        ++it;
+      }
+      
+      // if we do not found : restart
+      stop = (m_Seeds.size() != 0);
+      if (!stop) {
+	if (GetVerboseStep()) {
+	  std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
+	}
+        if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
+          // we want to skip more than a half of the image, it is probably a bug
+          std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
+          stop = true;
+        }
+        skip += 5;
+      }
+    }
+  }
+  return (m_Seeds.size() != 0);
 }
 //--------------------------------------------------------------------
 
   
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void 
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+TracheaRegionGrowing()
+{
+  // Explosion controlled region growing
+  typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
+  typename ImageFilterType::Pointer f= ImageFilterType::New();
+  f->SetInput(working_input);
+  f->SetVerbose(false);
+  f->SetLower(-2000);
+  f->SetUpper(GetUpperThresholdForTrachea());
+  f->SetMinimumLowerThreshold(-2000);
+  f->SetMaximumUpperThreshold(0);
+  f->SetAdaptLowerBorder(false);
+  f->SetAdaptUpperBorder(true);
+  f->SetMinimumSize(5000); 
+  f->SetReplaceValue(1);
+  f->SetMultiplier(GetMultiplierForTrachea());
+  f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
+  f->SetMinimumThresholdStepSize(1);
+  for(unsigned int i=0; i<m_Seeds.size();i++) {
+    f->AddSeed(m_Seeds[i]);
+    // DD(m_Seeds[i]);
+  }  
+  f->Update();
+
+  // take first (main) connected component
+  trachea_tmp = Labelize<InternalImageType>(f->GetOutput(), 
+					    GetBackgroundValue(), 
+					    true, 
+					    GetMinimalComponentSize());
+  trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp, 
+					      GetBackgroundValue(), 
+					      GetForegroundValue(), 
+					      1, 1, false);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+double 
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+ComputeTracheaVolume()
+{
+  typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
+  IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion());
+  iter.GoToBegin();
+  double volume = 0.0;
+  while (!iter.IsAtEnd()) {
+    if (iter.Get() == this->GetForegroundValue()) volume++;
+    ++iter;
+  }
+  
+  double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2];
+  return volume*voxelsize;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void 
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SearchForTrachea()
+{
+  // Search for seed among n slices, skip some slices before starting
+  // if not found -> skip more and restart 
+  
+  // when seed found : perform region growing
+  // compute trachea volume
+  // if volume not plausible  -> skip more slices and restart 
+
+  bool stop = false;
+  double volume = 0.0;
+  int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
+  while (!stop) {
+    stop = SearchForTracheaSeed(skip);
+    if (stop) {
+      TracheaRegionGrowing();
+      volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
+      if ((volume > 10) && (volume < 55 )) { // it is ok
+        // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
+	if (GetVerboseStep()) {
+	  std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
+	}
+        stop = true; 
+      }
+      else {
+	if (GetVerboseStep()) {
+	  std::cout << "\t The volume of the trachea (" << volume 
+		    << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds." 
+		    << std::endl;
+	}
+        skip += 5;
+        stop = false;
+	// empty the list of seed
+	m_Seeds.clear();
+      }
+    }
+  }
+
+  if (volume != 0.0) {
+    // Set output
+    StopCurrentStep<InternalImageType>(trachea_tmp);
+  }
+  else { // Trachea not found
+    this->SetWarning("* WARNING * No seed found for trachea.");
+    StopCurrentStep();
+  }
+}
+//--------------------------------------------------------------------
+
 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX