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

diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx
index 7dd27c7..596604b 100644
--- a/segmentation/clitkExtractLungFilter.txx
+++ b/segmentation/clitkExtractLungFilter.txx
@@ -47,6 +47,8 @@
 #include "itkOrientImageFilter.h"
 #include "itkSpatialOrientationAdapter.h"
 #include "itkImageDuplicator.h"
+#include "itkRelabelComponentImageFilter.h"
+
 #include <fcntl.h>
 
 //--------------------------------------------------------------------
@@ -67,6 +69,7 @@ ExtractLungFilter():
   SetForegroundValue(1);
   SetMinimalComponentSize(100);
   VerboseRegionGrowingFlagOff();
+  RemoveSmallLabelBeforeSeparationFlagOn();
 
   // Step 1 default values
   SetUpperThreshold(-300);
@@ -79,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);
@@ -127,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);
 }
@@ -254,6 +273,9 @@ GenerateOutputInformation()
   StartNewStep("Search for the trachea");
   SearchForTrachea();
   PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
+  if (m_Seeds.empty()) {
+    clitkExceptionMacro("No seeds for trachea... Aborting.");
+  }
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
@@ -428,7 +450,7 @@ GenerateOutputInformation()
     //--------------------------------------------------------------------
     //--------------------------------------------------------------------
     StartNewStep("Separate Left/Right lungs");
-      PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
+    PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
     // Initial label
     working_mask = Labelize<MaskImageType>(working_mask, 
                                           GetBackgroundValue(), 
@@ -443,9 +465,32 @@ GenerateOutputInformation()
     statisticsImageFilter->SetInput(working_mask);
     statisticsImageFilter->Update();
     unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
-    working_mask = statisticsImageFilter->GetOutput();	
-    
+    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();
+      }
+    }
   
     // Decompose the first label
     if (initialNumberOfLabels<2) {
@@ -542,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;
@@ -575,9 +620,6 @@ SearchForTracheaSeed(int skip)
 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
 {
   itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
-  std::cout << "orientation: " << std::hex << orientation << "; superior: " << std::hex << sup << std::endl;
-  std::cout << std::dec;
-
   bool primary = (orientation & 0x0000ff) == sup;
   bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
   bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
@@ -591,6 +633,9 @@ 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;
@@ -609,7 +654,7 @@ SearchForTracheaSeed2(int numberOfSlices)
     // threshold to isolate airawys and lungs
     typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
     threshold->SetLowerThreshold(-2000);
-    threshold->SetUpperThreshold(-400);
+    threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
     threshold->SetInput(working_input);
     threshold->Update();
     
@@ -623,6 +668,7 @@ SearchForTracheaSeed2(int numberOfSlices)
     opening->Update();
     
     typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
+    slicer->SetDirectionCollapseToIdentity();
     slicer->SetInput(opening->GetOutput());
     
     // label result
@@ -694,18 +740,17 @@ SearchForTracheaSeed2(int numberOfSlices)
         writer->SetFileName(file_name.str().c_str());
         writer->Update();
       }
-      
-      typename LabelImageType::LabelObjectContainerType shapes_map = label_map->GetLabelObjectContainer();
-      typename LabelImageType::LabelObjectContainerType::const_iterator s;
+
       typename ShapeLabelType::Pointer shape, max_e_shape;
-      double max_elongation = 0.5;
+      double max_elongation = GetMaxElongation();
       double max_size = size[0]*size[1]/128;
       double max_e = 0;
       int nshapes = 0;
-      for (s = shapes_map.begin(); s != shapes_map.end(); s++) {
-        shape = s->second;
-        if (shape->GetSize() > 150 && shape->GetSize() <= max_size) {
-          double e = 1 - 1/shape->GetBinaryElongation();
+      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++;
@@ -738,7 +783,7 @@ SearchForTracheaSeed2(int numberOfSlices)
         p2[2] = prev_e_centre[2];
         
         double mag = (p2 - p1).GetNorm();
-        if (GetVerboseStepFlag()) {
+        if (GetVerboseRegionGrowingFlag()) {
           cout.precision(3);
           cout << index[2] << ": ";
           cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
@@ -758,6 +803,11 @@ SearchForTracheaSeed2(int numberOfSlices)
         
         prev_e_centre= max_e_centre;
       }
+      else {
+        if (GetVerboseRegionGrowingFlag()) {
+          cout << "No shapes found at slice " << index[2] << std::endl;
+        }
+      }
     }
     
     size_t longest = 0;
@@ -770,9 +820,11 @@ SearchForTracheaSeed2(int numberOfSlices)
       }
     }
     
-    if (GetVerboseStepFlag()) 
-      std::cout << "seed at: " << trachea_centre << std::endl;
-    m_Seeds.push_back(trachea_centre);
+    if (longest > 0) {
+      if (GetVerboseRegionGrowingFlag()) 
+        std::cout << "seed at: " << trachea_centre << std::endl;
+      m_Seeds.push_back(trachea_centre);
+    }
   }
   
   return (m_Seeds.size() != 0);
@@ -794,7 +846,7 @@ TracheaRegionGrowing()
   f->SetUpper(GetUpperThresholdForTrachea());
   f->SetMinimumLowerThreshold(-2000);
   //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
-  f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
+  f->SetMaximumUpperThreshold(-300); // MAYBE TO CHANGE ???
   f->SetAdaptLowerBorder(false);
   f->SetAdaptUpperBorder(true);
   f->SetMinimumSize(5000); 
@@ -805,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");
@@ -859,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 = true;
-    if (SearchForTracheaSeed2(50)) {
+    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()) {
@@ -878,16 +945,17 @@ SearchForTrachea()
             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
           }
         }
-        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