X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=596604b61f44b1f6e1c8c39115f3888a4f069734;hb=HEAD;hp=67d3f5e641c08bb357a57c0b51fc17805287de45;hpb=573d80d0f7a17607d2ee883c21c940c0ba020282;p=clitk.git

diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx
index 67d3f5e..596604b 100644
--- a/segmentation/clitkExtractLungFilter.txx
+++ b/segmentation/clitkExtractLungFilter.txx
@@ -37,6 +37,19 @@
 #include "itkImageIteratorWithIndex.h"
 #include "itkBinaryMorphologicalOpeningImageFilter.h"
 #include "itkBinaryMorphologicalClosingImageFilter.h"
+#include "itkConstantPadImageFilter.h"
+#include <itkBinaryBallStructuringElement.h>
+#include "itkStatisticsLabelObject.h"
+#include "itkLabelMap.h"
+#include "itkLabelImageToShapeLabelMapFilter.h"
+#include "itkLabelMapToLabelImageFilter.h"
+#include "itkExtractImageFilter.h"
+#include "itkOrientImageFilter.h"
+#include "itkSpatialOrientationAdapter.h"
+#include "itkImageDuplicator.h"
+#include "itkRelabelComponentImageFilter.h"
+
+#include <fcntl.h>
 
 //--------------------------------------------------------------------
 template <class ImageType>
@@ -56,6 +69,7 @@ ExtractLungFilter():
   SetForegroundValue(1);
   SetMinimalComponentSize(100);
   VerboseRegionGrowingFlagOff();
+  RemoveSmallLabelBeforeSeparationFlagOn();
 
   // Step 1 default values
   SetUpperThreshold(-300);
@@ -68,11 +82,15 @@ ExtractLungFilter():
   SetLabelizeParameters1(p1);
 
   // Step 2 default values
-  SetUpperThresholdForTrachea(-900);
+  SetTracheaSeedAlgorithm(0);
+  SetUpperThresholdForTrachea(-700);
   SetMultiplierForTrachea(5);
   SetThresholdStepSizeForTrachea(64);
   SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
   TracheaVolumeMustBeCheckedFlagOn();
+  SetNumSlices(50);
+  SetMaxElongation(0.5);
+  SetSeedPreProcessingThreshold(-400); 
   
   // Step 3 default values
   SetNumberOfHistogramBins(500);
@@ -96,6 +114,7 @@ ExtractLungFilter():
   
   // Step 6
   FillHolesFlagOn();
+  AutoCropOn();
 }
 //--------------------------------------------------------------------
 
@@ -115,7 +134,19 @@ SetInput(const ImageType * image)
 template <class ImageType>
 void 
 clitk::ExtractLungFilter<ImageType>::
-AddSeed(InternalIndexType s) 
+//AddSeed(InternalIndexType s) 
+AddSeed(InputImagePointType s) 
+{ 
+  m_SeedsInMM.push_back(s);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::ExtractLungFilter<ImageType>::
+AddSeedInPixels(InternalIndexType s) 
 { 
   m_Seeds.push_back(s);
 }
@@ -160,6 +191,22 @@ GenerateOutputInformation()
   StartNewStep("Set background to initial image");
   working_input = SetBackground<ImageType, MaskImageType>
     (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
+
+  // Pad images with air to prevent patient touching the image border
+  static const unsigned int Dim = ImageType::ImageDimension;
+  typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
+  typename PadFilterType::Pointer padFilter = PadFilterType::New();
+  padFilter->SetInput(working_input);
+  padFilter->SetConstant(-1000);
+  typename ImageType::SizeType bounds;
+  for (unsigned i = 0; i < Dim - 1; ++i)
+    bounds[i] = 1;
+  bounds[Dim - 1] = 0;
+  padFilter->SetPadLowerBound(bounds);
+  padFilter->SetPadUpperBound(bounds);
+  padFilter->Update();
+  working_input = padFilter->GetOutput();
+
   StopCurrentStep<ImageType>(working_input);
   PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
 
@@ -226,6 +273,9 @@ GenerateOutputInformation()
   StartNewStep("Search for the trachea");
   SearchForTrachea();
   PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
+  if (m_Seeds.empty()) {
+    clitkExceptionMacro("No seeds for trachea... Aborting.");
+  }
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
@@ -311,7 +361,19 @@ GenerateOutputInformation()
   //--------------------------------------------------------------------
   PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
-    trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
+    if (GetAutoCrop())
+      trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
+    else
+    {
+      // Remove Padding region
+      typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
+      typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+      cropFilter->SetInput(trachea);
+      cropFilter->SetLowerBoundaryCropSize(bounds);
+      cropFilter->SetUpperBoundaryCropSize(bounds);
+      cropFilter->Update();
+      trachea = cropFilter->GetOutput();
+    }
     StopCurrentStep<MaskImageType>(trachea);  
     PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
   }
@@ -321,7 +383,19 @@ GenerateOutputInformation()
   //--------------------------------------------------------------------
   StartNewStep("Cropping lung");
   PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
-  working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
+  if (GetAutoCrop())
+    working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
+  else
+  {
+    // Remove Padding region
+    typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
+    typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+    cropFilter->SetInput(working_mask);
+    cropFilter->SetLowerBoundaryCropSize(bounds);
+    cropFilter->SetUpperBoundaryCropSize(bounds);
+    cropFilter->Update();
+    working_mask = cropFilter->GetOutput();
+  }
   StopCurrentStep<MaskImageType>(working_mask);
 
   //--------------------------------------------------------------------
@@ -372,50 +446,74 @@ GenerateOutputInformation()
     StopCurrentStep<MaskImageType>(working_mask);
   }
 
-  //--------------------------------------------------------------------
-  //--------------------------------------------------------------------
-  StartNewStep("Separate Left/Right lungs");
+  if (GetSeparateLungsFlag()) {
+    //--------------------------------------------------------------------
+    //--------------------------------------------------------------------
+    StartNewStep("Separate Left/Right lungs");
     PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
-  // Initial label
-  working_mask = Labelize<MaskImageType>(working_mask, 
-                                         GetBackgroundValue(), 
-                                         false, 
-                                         GetMinimalComponentSize());
-
-  PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
-
-  // Count the labels
-  typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
-  typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
-  statisticsImageFilter->SetInput(working_mask);
-  statisticsImageFilter->Update();
-  unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
-  working_mask = statisticsImageFilter->GetOutput();	
+    // Initial label
+    working_mask = Labelize<MaskImageType>(working_mask, 
+                                          GetBackgroundValue(), 
+                                          false, 
+                                          GetMinimalComponentSize());
+
+    PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
+
+    // Count the labels
+    typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
+    typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
+    statisticsImageFilter->SetInput(working_mask);
+    statisticsImageFilter->Update();
+    unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
+    working_mask = statisticsImageFilter->GetOutput();	    
+    PrintMemory(GetVerboseMemoryFlag(), "After count label");
+
+    // If already 2 labels, but a too big differences, remove the
+    // smalest one (sometimes appends with the stomach
+    if (initialNumberOfLabels >1) {
+      if (GetRemoveSmallLabelBeforeSeparationFlag()) {
+        typedef itk::RelabelComponentImageFilter<MaskImageType, MaskImageType> RelabelFilterType;
+        typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
+        relabelFilter->SetInput(working_mask);
+        relabelFilter->SetMinimumObjectSize(10);
+        relabelFilter->Update();
+        const std::vector<float> & a = relabelFilter->GetSizeOfObjectsInPhysicalUnits();
+        std::vector<MaskImagePixelType> remove_label;
+        for(unsigned int i=1; i<a.size(); i++) {
+          if (a[i] < 0.5*a[0]) { // more than 0.5 difference
+            remove_label.push_back(i+1); // label zero is BG
+          }
+        }
+        working_mask = 
+          clitk::RemoveLabels<MaskImageType>(working_mask, GetBackgroundValue(), remove_label);
+        statisticsImageFilter->SetInput(working_mask);
+        statisticsImageFilter->Update();
+        initialNumberOfLabels = statisticsImageFilter->GetMaximum();
+      }
+    }
   
-  PrintMemory(GetVerboseMemoryFlag(), "After count label");
- 
-  // Decompose the first label
-  static const unsigned int Dim = ImageType::ImageDimension;
-  if (initialNumberOfLabels<2) {
-    // Structuring element radius
-    typename ImageType::SizeType radius;
-    for (unsigned int i=0;i<Dim;i++) radius[i]=1;
-    typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
-    typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
-    decomposeAndReconstructFilter->SetInput(working_mask);
-    decomposeAndReconstructFilter->SetVerbose(false);
-    decomposeAndReconstructFilter->SetRadius(radius);
-    decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
-    decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
-    decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
-    decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
-    decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
-    decomposeAndReconstructFilter->SetFullyConnected(true);
-    decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
-    decomposeAndReconstructFilter->Update();
-    working_mask = decomposeAndReconstructFilter->GetOutput();      
+    // Decompose the first label
+    if (initialNumberOfLabels<2) {
+      // Structuring element radius
+      typename ImageType::SizeType radius;
+      for (unsigned int i=0;i<Dim;i++) radius[i]=1;
+      typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
+      typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
+      decomposeAndReconstructFilter->SetInput(working_mask);
+      decomposeAndReconstructFilter->SetVerbose(false);
+      decomposeAndReconstructFilter->SetRadius(radius);
+      decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
+      decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
+      decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
+      decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
+      decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
+      decomposeAndReconstructFilter->SetFullyConnected(true);
+      decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
+      decomposeAndReconstructFilter->Update();
+      working_mask = decomposeAndReconstructFilter->GetOutput();      
+    }
+    PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
   }
-  PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
 
   // Retain labels ('1' is largset lung, so right. '2' is left)
   typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
@@ -489,7 +587,7 @@ SearchForTracheaSeed(int skip)
       it.GoToBegin();
       while (!it.IsAtEnd()) {
         if(it.Get() < GetUpperThresholdForTrachea() ) {
-          AddSeed(it.GetIndex());
+          AddSeedInPixels(it.GetIndex());
           // DD(it.GetIndex());
         }
         ++it;
@@ -518,7 +616,221 @@ SearchForTracheaSeed(int skip)
 }
 //--------------------------------------------------------------------
 
+
+bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
+{
+  itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
+  bool primary = (orientation & 0x0000ff) == sup;
+  bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
+  bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
+  return primary || secondary || tertiary;
+}
+
+//--------------------------------------------------------------------
+template <class ImageType>
+bool 
+clitk::ExtractLungFilter<ImageType>::
+SearchForTracheaSeed2(int numberOfSlices)
+{
+  if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)    
+    if (GetVerboseRegionGrowingFlag())
+      std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
+    
+    typedef unsigned char MaskPixelType;
+    typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
+    typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
+    typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
+    typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
+    typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
+    typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
+    typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
+    typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
+    typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
+    typedef itk::LabelMap<ShapeLabelType> LabelImageType;
+    typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
+    typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
+    typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
+    
+    // threshold to isolate airawys and lungs
+    typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
+    threshold->SetLowerThreshold(-2000);
+    threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
+    threshold->SetInput(working_input);
+    threshold->Update();
+    
+    KernelType kernel_closing, kernel_opening;
+    
+    // remove small noise
+    typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
+    kernel_opening.SetRadius(1);
+    opening->SetKernel(kernel_opening);
+    opening->SetInput(threshold->GetOutput());
+    opening->Update();
+    
+    typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
+    slicer->SetDirectionCollapseToIdentity();
+    slicer->SetInput(opening->GetOutput());
+    
+    // label result
+    typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
+    label_filter->SetInput(slicer->GetOutput());
+    
+    // extract shape information from labels
+    typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
+    label_to_map_filter->SetInput(label_filter->GetOutput());
+    
+    typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
+    typename FileWriterType::Pointer writer = FileWriterType::New();
+    
+    typename ImageType::IndexType index;
+    typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
+    typename ImageType::SizeType size = region.GetSize();
+    typename ImageType::SpacingType spacing = working_input->GetSpacing();
+    typename ImageType::PointType origin = working_input->GetOrigin();
+
+    int nslices = min(numberOfSlices, size[2]);
+    int start = 0, increment = 1;
+    itk::SpatialOrientationAdapter orientation;
+    typename ImageType::DirectionType dir = working_input->GetDirection();
+    if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
+      start = size[2]-1;
+      increment = -1;
+    }
+    
+    typename MaskImageType::PointType image_centre;
+    image_centre[0] = size[0]/2;
+    image_centre[1] = size[1]/2;
+    image_centre[2] = 0;
+  
+    typedef InternalIndexType SeedType;
+    SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
+    typedef std::list<SeedType> PointListType;
+    typedef std::list<PointListType> SequenceListType;
+    PointListType* current_sequence = NULL;
+    SequenceListType sequence_list;
+
+    prev_e_centre.Fill(0);
+    std::ostringstream file_name;
+    index[0] = index[1] = 0;
+    size[0] = size[1] = 512;
+    size[2] = 0;
+    while (nslices--) {
+      index[2] = start;
+      start += increment;
+      
+      region.SetIndex(index);
+      region.SetSize(size);
+      slicer->SetExtractionRegion(region);
+      slicer->Update();
+      label_filter->SetInput(slicer->GetOutput());
+      label_filter->Update();
+
+      label_to_map_filter->SetInput(label_filter->GetOutput());
+      label_to_map_filter->Update();
+      typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
+
+      if (GetWriteStepFlag()) {
+        map_to_label_filter->SetInput(label_map);
+        writer->SetInput(map_to_label_filter->GetOutput());
+        file_name.str("");
+        file_name << "labels_";
+        file_name.width(3);
+        file_name.fill('0');
+        file_name << index[2] << ".mhd";
+        writer->SetFileName(file_name.str().c_str());
+        writer->Update();
+      }
+
+      typename ShapeLabelType::Pointer shape, max_e_shape;
+      double max_elongation = GetMaxElongation();
+      double max_size = size[0]*size[1]/128;
+      double max_e = 0;
+      int nshapes = 0;
+      unsigned int nlables = label_map->GetNumberOfLabelObjects();
+      for (unsigned int j = 0; j < nlables; j++) {
+        shape = label_map->GetNthLabelObject(j);
+        if (shape->Size() > 150 && shape->Size() <= max_size) {
+          double e = 1 - 1/shape->GetElongation();
+          //double area = 1 - r->Area() ;
+          if (e < max_elongation) {
+            nshapes++;
+            shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
+            shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
+            shape_centre[2] = index[2];
+            //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
+            double dx = shape_centre[0] - image_centre[0];
+            double d = 1 - dx*2/size[0];
+            e = e + d;
+            if (e > max_e)
+            {
+              max_e = e;
+              max_e_shape = shape;
+              max_e_centre = shape_centre;
+            }
+          }
+        }
+      }
+      
+      if (nshapes > 0)
+      {
+        itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
+        p1[0] = max_e_centre[0];
+        p1[1] = max_e_centre[1];
+        p1[2] = max_e_centre[2];
+        
+        p2[0] = prev_e_centre[0];
+        p2[1] = prev_e_centre[1];
+        p2[2] = prev_e_centre[2];
+        
+        double mag = (p2 - p1).GetNorm();
+        if (GetVerboseRegionGrowingFlag()) {
+          cout.precision(3);
+          cout << index[2] << ": ";
+          cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
+          cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
+          cout << "mag(" << mag << "); " << endl;
+        }
+        
+        if (mag > 5)
+        {
+          PointListType point_list;
+          point_list.push_back(max_e_centre);
+          sequence_list.push_back(point_list);
+          current_sequence = &sequence_list.back();
+        }
+        else if (current_sequence)
+          current_sequence->push_back(max_e_centre);
+        
+        prev_e_centre= max_e_centre;
+      }
+      else {
+        if (GetVerboseRegionGrowingFlag()) {
+          cout << "No shapes found at slice " << index[2] << std::endl;
+        }
+      }
+    }
+    
+    size_t longest = 0;
+    for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
+    {
+      if (s->size() > longest)
+      {
+        longest = s->size();
+        trachea_centre = s->front();
+      }
+    }
+    
+    if (longest > 0) {
+      if (GetVerboseRegionGrowingFlag()) 
+        std::cout << "seed at: " << trachea_centre << std::endl;
+      m_Seeds.push_back(trachea_centre);
+    }
+  }
   
+  return (m_Seeds.size() != 0);
+}
+//--------------------------------------------------------------------
+
 //--------------------------------------------------------------------
 template <class ImageType>
 void 
@@ -534,7 +846,7 @@ TracheaRegionGrowing()
   f->SetUpper(GetUpperThresholdForTrachea());
   f->SetMinimumLowerThreshold(-2000);
   //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
-  f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
+  f->SetMaximumUpperThreshold(-300); // MAYBE TO CHANGE ???
   f->SetAdaptLowerBorder(false);
   f->SetAdaptUpperBorder(true);
   f->SetMinimumSize(5000); 
@@ -545,7 +857,7 @@ TracheaRegionGrowing()
   f->SetVerbose(GetVerboseRegionGrowingFlag());
   for(unsigned int i=0; i<m_Seeds.size();i++) {
     f->AddSeed(m_Seeds[i]);
-    // DD(m_Seeds[i]);
+    //DD(m_Seeds[i]);
   }  
   f->Update();
   PrintMemory(GetVerboseMemoryFlag(), "After RG update");
@@ -588,7 +900,7 @@ ComputeTracheaVolume()
 
 //--------------------------------------------------------------------
 template <class ImageType>
-void 
+void
 clitk::ExtractLungFilter<ImageType>::
 SearchForTrachea()
 {
@@ -599,12 +911,27 @@ SearchForTrachea()
   // compute trachea volume
   // if volume not plausible  -> skip more slices and restart 
 
+  // If initial seed, convert from mm to pixels
+  if (m_SeedsInMM.size() > 0) {
+    for(unsigned int i=0; i<m_SeedsInMM.size(); i++) {
+      InputImageIndexType index;
+      working_input->TransformPhysicalPointToIndex(m_SeedsInMM[i], index);
+      m_Seeds.push_back(index);
+    }
+  }
+
+  bool has_seed;
   bool stop = false;
   double volume = 0.0;
   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
   while (!stop) {
-    stop = SearchForTracheaSeed(skip);
-    if (stop) {
+    stop = true;
+    if (GetTracheaSeedAlgorithm() == 0)
+      has_seed = SearchForTracheaSeed(skip);
+    else
+      has_seed = SearchForTracheaSeed2(GetNumSlices());
+    
+    if (has_seed) {
       TracheaRegionGrowing();
       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
       if (GetWriteStepFlag()) {
@@ -613,21 +940,27 @@ SearchForTrachea()
       if (GetTracheaVolumeMustBeCheckedFlag()) {
         if ((volume > 10) && (volume < 65 )) { // depend on image size ...
           // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
-          if (GetVerboseStepFlag()) {
+          if (GetVerboseStepFlag())
+          {
             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
           }
-          stop = true; 
         }
-        else {
-          if (GetVerboseStepFlag()) {
-            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();
+        else 
+          if (GetTracheaSeedAlgorithm() == 0) {
+            if (GetVerboseStepFlag()) {
+              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 (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 << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
+          stop = true;
         }
       }
       else {