X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkMotionMaskGenericFilter.txx;h=fd21b4629d4da99590b741e4fe7fa9ff634cfb50;hb=ae015f09e2fa0ebc736d24b37c9ed6c1ca0cb5b2;hp=414a2f69042a5494b2355472c5785f3d31604efd;hpb=08d1fd56ac1d08bd228d9e557f5472a395e9b708;p=clitk.git

diff --git a/segmentation/clitkMotionMaskGenericFilter.txx b/segmentation/clitkMotionMaskGenericFilter.txx
index 414a2f6..fd21b46 100644
--- a/segmentation/clitkMotionMaskGenericFilter.txx
+++ b/segmentation/clitkMotionMaskGenericFilter.txx
@@ -27,6 +27,8 @@
  *
  ===================================================*/
 
+#include "itkLabelImageToShapeLabelMapFilter.h"
+#include "itkShapeLabelObject.h"
 
 namespace clitk
 {
@@ -377,6 +379,8 @@ MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>
     connectFilter->SetBackgroundValue(0);
     connectFilter->SetFullyConnected(true);
     if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
+    connectFilter->Update();
+    if (m_Verbose) std::cout<<"found "<< connectFilter->GetObjectCount() << std::endl;
 
     //---------------------------------
     // Sort the labels according to size
@@ -384,7 +388,7 @@ MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>
     typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
     relabelFilter->SetInput(connectFilter->GetOutput());
     if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
-    // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/jef/tmp/labels.mhd");
+    // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/vdelmon/tmp/labels.mhd");
 
     //---------------------------------
     // Keep the label
@@ -478,7 +482,7 @@ MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimensi
 //-------------------------------------------------------------------
 template <unsigned int Dimension, class  PixelType>
 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
-MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
+MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low, typename itk::Image<InternalPixelType,Dimension>::Pointer lungs_low )
 {
 
   // ImageTypes
@@ -509,24 +513,85 @@ MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension
       std::cout<<"=========================================="<<std::endl;
     }
 
-    //---------------------------------
-    // Offset from center
-    //---------------------------------
-    typename itk::Vector<double,Dimension> offset;
-    if(m_ArgsInfo.offset_given== Dimension) {
-      for(unsigned int i=0; i<Dimension; i++)
-        offset[i]=m_ArgsInfo.offset_arg[i];
-    } else {
-      offset.Fill(0.);
-      offset[1]=-50;
+    itk::Vector<double,Dimension> centerEllips;
+    typename itk::Vector<double, Dimension> axes;
+    if (m_ArgsInfo.ellipseAutoDetect_flag) {
+      if(m_Verbose) {
+        std::cout<<"Auto-detecting intial ellipse..."<<std::endl;
+      }
+
+      typename InternalImageType::RegionType region, region_before = lungs_low->GetLargestPossibleRegion();
+      typename InternalImageType::SizeType size = region_before.GetSize();
+      size[2] /= 2;
+      region.SetSize(size);
+      lungs_low->SetRegions(region);
+      
+      // compute statistics of the (mirroed) lung region in order to guess the params of the ellipse
+      typedef itk::LabelImageToShapeLabelMapFilter<InternalImageType>  LabelImageToShapeLabelMapFilter;
+      typename LabelImageToShapeLabelMapFilter::Pointer label_filter = LabelImageToShapeLabelMapFilter::New();
+      label_filter->SetInput(lungs_low);
+      label_filter->Update();
+      
+      typename InternalImageType::SpacingType spacing = lungs_low->GetSpacing();
+      typedef typename LabelImageToShapeLabelMapFilter::OutputImageType::LabelObjectType LabelType;
+      const LabelType* label = dynamic_cast<LabelType*>(label_filter->GetOutput()->GetNthLabelObject(0));
+
+      // lungs' barycenter
+      typename LabelType::CentroidType lung_centroid = label->GetCentroid();
+
+      // try to guess ideal ellipse axes from the lungs' bounding box and centroid
+      // with some hard-coded "magic" constants...
+      typename LabelType::RegionType lung_bbox = label->GetBoundingBox();
+
+      axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2;
+      axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2;
+      axes[2] = 1.25*fabs(lung_centroid[2] - center[2]);
+
+      // ellipse's center in XY is the center of the lungs' bounding box
+      typename InternalImageType::PointType origin = lungs_low->GetOrigin();
+      centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0];
+      centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1];
+      centerEllips[2] = center[2];
+      
+      lungs_low->SetRegions(region_before);
+      
+      if(m_Verbose) {
+        std::cout << "Lungs'centroid at " << lung_centroid << std::endl;
+        std::cout << "Ellipse center at " << centerEllips << std::endl;
+        std::cout << "Ellipse axes as " << axes << std::endl;
+        std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl;
+      }
     }
-    itk::Vector<double,Dimension> centerEllips=center+offset;
-    if (m_Verbose) {
-      std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
-      std::cout<<centerEllips[0];
-      for (unsigned int i=1; i<Dimension; i++)
-        std::cout<<" "<<centerEllips[i];
-      std::cout<<std::endl;
+    else {
+      //---------------------------------
+      // Offset from center
+      //---------------------------------
+      typename itk::Vector<double,Dimension> offset;
+      if(m_ArgsInfo.offset_given== Dimension) {
+        for(unsigned int i=0; i<Dimension; i++)
+          offset[i]=m_ArgsInfo.offset_arg[i];
+      } else {
+        offset.Fill(0.);
+        offset[1]=-50;
+      }
+      centerEllips=center+offset;
+      if (m_Verbose) {
+        std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
+        std::cout<<centerEllips[0];
+        for (unsigned int i=1; i<Dimension; i++)
+          std::cout<<" "<<centerEllips[i];
+        std::cout<<std::endl;
+      }
+
+      // Axes
+      if (m_ArgsInfo.axes_given==Dimension)
+        for(unsigned int i=0; i<Dimension; i++)
+          axes[i]=m_ArgsInfo.axes_arg[i];
+      else {
+        axes[0]=100;
+        axes[1]=30;
+        axes[2]=150;
+      }
     }
 
     //---------------------------------
@@ -539,17 +604,6 @@ MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension
     ellips->Allocate();
     ellips->FillBuffer(0);
 
-    // Axes
-    typename itk::Vector<double, Dimension> axes;
-    if (m_ArgsInfo.axes_given==Dimension)
-      for(unsigned int i=0; i<Dimension; i++)
-        axes[i]=m_ArgsInfo.axes_arg[i];
-    else {
-      axes[0]=100;
-      axes[1]=30;
-      axes[2]=150;
-    }
-
     // Draw an ellips
     IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
     itEllips.GoToBegin();
@@ -725,7 +779,7 @@ MotionMaskGenericFilter::UpdateWithDimAndPixelType()
   //----------------------------------------------------------------------------------------------------
   // Ellipsoide initialization
   //----------------------------------------------------------------------------------------------------
-  typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
+  typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low,lungs_low);
 
   //----------------------------------------------------------------------------------------------------
   // Grow ellips
@@ -944,6 +998,10 @@ MotionMaskGenericFilter::UpdateWithDimAndPixelType()
     distanceMapImageFilter->SetInsideIsPositive(false);
     if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
     distanceMapImageFilter->Update();
+    if (m_ArgsInfo.writeDistMap_given) {
+      writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
+      
+    }
 
     //---------------------------------
     // Grow while monitoring detection point
@@ -1069,6 +1127,10 @@ MotionMaskGenericFilter::UpdateWithDimAndPixelType()
     distanceMapImageFilter->SetInsideIsPositive(false);
     if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
     distanceMapImageFilter->Update();
+    if (m_ArgsInfo.writeDistMap_given) {
+      writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
+      
+    }
 
     //---------------------------------
     // Grow while monitoring lung volume coverage