X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=d2bff5edd5962647512fd2d1070017b17a5dab34;hb=e2b37672d5ee8eafc7b1ac075f4e70596349f0c3;hp=1a945917be9149b76617dfca22c40de34f533db6;hpb=cf35c87a6dd45a7daf171d2adca9766846a0f127;p=clitk.git

diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx
index 1a94591..d2bff5e 100644
--- a/segmentation/clitkExtractLungFilter.txx
+++ b/segmentation/clitkExtractLungFilter.txx
@@ -37,6 +37,7 @@
 #include "itkImageIteratorWithIndex.h"
 #include "itkBinaryMorphologicalOpeningImageFilter.h"
 #include "itkBinaryMorphologicalClosingImageFilter.h"
+#include "itkConstantPadImageFilter.h"
 
 //--------------------------------------------------------------------
 template <class ImageType>
@@ -161,6 +162,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
 
@@ -314,6 +331,17 @@ GenerateOutputInformation()
   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
     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");
   }
@@ -325,6 +353,17 @@ GenerateOutputInformation()
   PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
   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);
 
   //--------------------------------------------------------------------
@@ -375,50 +414,51 @@ GenerateOutputInformation()
     StopCurrentStep<MaskImageType>(working_mask);
   }
 
-  //--------------------------------------------------------------------
-  //--------------------------------------------------------------------
-  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();	
+  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();	
+    
+    PrintMemory(GetVerboseMemoryFlag(), "After count label");
   
-  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;
@@ -591,7 +631,7 @@ ComputeTracheaVolume()
 
 //--------------------------------------------------------------------
 template <class ImageType>
-void 
+void
 clitk::ExtractLungFilter<ImageType>::
 SearchForTrachea()
 {
@@ -606,8 +646,8 @@ SearchForTrachea()
   double volume = 0.0;
   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
   while (!stop) {
-    stop = SearchForTracheaSeed(skip);
-    if (stop) {
+    stop = true;
+    if (SearchForTracheaSeed(skip)) {
       TracheaRegionGrowing();
       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
       if (GetWriteStepFlag()) {
@@ -619,7 +659,6 @@ SearchForTrachea()
           if (GetVerboseStepFlag()) {
             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
           }
-          stop = true; 
         }
         else {
           if (GetVerboseStepFlag()) {
@@ -632,6 +671,11 @@ SearchForTrachea()
           // 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 {
         stop = true;