]> Creatis software - clitk.git/blobdiff - segmentation/clitkExtractMediastinumFilter.txx
MS compiler compatibility issues
[clitk.git] / segmentation / clitkExtractMediastinumFilter.txx
index 854a8042b4c33c00d9258fe031a81aafe70ea683..b76e617946593c7faaca4924beb61f6d7a43e796 100644 (file)
 #include "clitkExtractMediastinumFilter.h"
 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
 #include "clitkSegmentationUtils.h"
+#include "clitkExtractAirwaysTreeInfoFilter.h"
+#include "clitkCropLikeImageFilter.h"
 
-// itk
+// std
 #include <deque>
+
+// itk
 #include "itkStatisticsLabelMapFilter.h"
 #include "itkLabelImageToStatisticsLabelMapFilter.h"
 #include "itkRegionOfInterestImageFilter.h"
 #include "RelativePositionPropImageFilter.h"
 
 //--------------------------------------------------------------------
-template <class TImageType>
-clitk::ExtractMediastinumFilter<TImageType>::
+template <class ImageType>
+clitk::ExtractMediastinumFilter<ImageType>::
 ExtractMediastinumFilter():
   clitk::FilterBase(),
-  itk::ImageToImageFilter<TImageType, TImageType>()
+  clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
+  itk::ImageToImageFilter<ImageType, ImageType>()
 {
-  this->SetNumberOfRequiredInputs(3);
+  this->SetNumberOfRequiredInputs(4);
 
   SetBackgroundValuePatient(0);
   SetBackgroundValueLung(0);
@@ -53,30 +58,33 @@ ExtractMediastinumFilter():
   SetForegroundValue(1);
 
   SetIntermediateSpacing(6);
-  SetFuzzyThreshold1(0.6);
-  SetFuzzyThreshold2(0.7);
+  SetFuzzyThreshold1(0.4);
+  SetFuzzyThreshold2(0.6);
+  SetFuzzyThreshold3(0.2);
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
-SetInputPatientLabelImage(const TImageType * image, ImagePixelType bg) {
-  this->SetNthInput(0, const_cast<TImageType *>(image));
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputPatientLabelImage(const ImageType * image, ImagePixelType bg) 
+{
+  this->SetNthInput(0, const_cast<ImageType *>(image));
   m_BackgroundValuePatient = bg;
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
-SetInputLungLabelImage(const TImageType * image, ImagePixelType bg, 
-                       ImagePixelType fgLeft, ImagePixelType fgRight) {
-  this->SetNthInput(1, const_cast<TImageType *>(image));
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputLungLabelImage(const ImageType * image, ImagePixelType bg, 
+                       ImagePixelType fgLeft, ImagePixelType fgRight) 
+{
+  this->SetNthInput(1, const_cast<ImageType *>(image));
   m_BackgroundValueLung = bg;
   m_ForegroundValueLeftLung = fgLeft;
   m_ForegroundValueRightLung = fgRight;
@@ -85,21 +93,34 @@ SetInputLungLabelImage(const TImageType * image, ImagePixelType bg,
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
-SetInputBonesLabelImage(const TImageType * image, ImagePixelType bg) {
-  this->SetNthInput(2, const_cast<TImageType *>(image));
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputBonesLabelImage(const ImageType * image, ImagePixelType bg) 
+{
+  this->SetNthInput(2, const_cast<ImageType *>(image));
   m_BackgroundValueBones = bg;
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
+void 
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputTracheaLabelImage(const ImageType * image, ImagePixelType bg) 
+{
+  this->SetNthInput(3, const_cast<ImageType *>(image));
+  m_BackgroundValueTrachea = bg;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
 template<class ArgsInfoType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
+clitk::ExtractMediastinumFilter<ImageType>::
 SetArgsInfo(ArgsInfoType mArgsInfo)
 {
   SetVerboseOption_GGO(mArgsInfo);
@@ -109,7 +130,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
 
   SetBackgroundValuePatient_GGO(mArgsInfo);
   SetBackgroundValueLung_GGO(mArgsInfo);
-  SetBackgroundValueBones_GGO(mArgsInfo);
+  SetBackgroundValueTrachea_GGO(mArgsInfo);
 
   SetForegroundValueLeftLung_GGO(mArgsInfo);
   SetForegroundValueRightLung_GGO(mArgsInfo);
@@ -117,16 +138,19 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
   SetIntermediateSpacing_GGO(mArgsInfo);
   SetFuzzyThreshold1_GGO(mArgsInfo);
   SetFuzzyThreshold2_GGO(mArgsInfo);
+  SetFuzzyThreshold3_GGO(mArgsInfo);
+
+  SetAFDBFilename_GGO(mArgsInfo);
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
+clitk::ExtractMediastinumFilter<ImageType>::
 GenerateOutputInformation() { 
-  ImagePointer input = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
+  ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
   ImagePointer outputImage = this->GetOutput(0);
   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
 }
@@ -134,43 +158,69 @@ GenerateOutputInformation() {
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
-GenerateInputRequestedRegion() {
+clitk::ExtractMediastinumFilter<ImageType>::
+GenerateInputRequestedRegion() 
+{
   // Call default
   Superclass::GenerateInputRequestedRegion();  
-  // Get input pointers
-  ImagePointer patient = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
-  ImagePointer lung    = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
-  ImagePointer bones   = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(2));
+
+  // Get input pointers  
+  LoadAFDB();
+  ImagePointer patient = GetAFDB()->template GetImage <ImageType>("patient");  
+  ImagePointer lung = GetAFDB()->template GetImage <ImageType>("lungs");  
+  ImagePointer bones = GetAFDB()->template GetImage <ImageType>("bones");  
+  ImagePointer trachea = GetAFDB()->template GetImage <ImageType>("trachea");  
     
   patient->SetRequestedRegion(patient->GetLargestPossibleRegion());
   lung->SetRequestedRegion(lung->GetLargestPossibleRegion());
   bones->SetRequestedRegion(bones->GetLargestPossibleRegion());
+  trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
 void 
-clitk::ExtractMediastinumFilter<TImageType>::
-GenerateData() {
+clitk::ExtractMediastinumFilter<ImageType>::
+GenerateData() 
+{
   // Get input pointers
-  ImageConstPointer patient = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(0));
-  ImageConstPointer lung   = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(1));
-  ImageConstPointer bones   = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(2));
+  ImagePointer patient = GetAFDB()->template GetImage <ImageType>("patient");  
+  ImagePointer lung = GetAFDB()->template GetImage <ImageType>("lungs");  
+  ImagePointer bones = GetAFDB()->template GetImage <ImageType>("bones");  
+  ImagePointer trachea = GetAFDB()->template GetImage <ImageType>("trachea");  
     
   // Get output pointer
   ImagePointer output;
 
+  // Step 0: Crop support (patient) to lung extend in RL
+  StartNewStep("Crop support like lungs along LR");
+  typedef clitk::CropLikeImageFilter<ImageType> CropFilterType;
+  typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+  cropFilter->SetInput(patient);
+  cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
+  cropFilter->Update();
+  output = cropFilter->GetOutput();
+  this->template StopCurrentStep<ImageType>(output);
+
+  // Step 0: Crop support (previous) to bones extend in AP
+  StartNewStep("Crop support like bones along AP");
+  cropFilter = CropFilterType::New();
+  cropFilter->SetInput(output);
+  cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
+  cropFilter->Update();
+  output = cropFilter->GetOutput();
+  this->template StopCurrentStep<ImageType>(output);
+
   // Step 1: patient minus lungs, minus bones
   StartNewStep("Patient contours minus lungs and minus bones");
   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
   boolFilter->InPlaceOn();
-  boolFilter->SetInput1(patient);
+  boolFilter->SetInput1(output);
   boolFilter->SetInput2(lung);    
   boolFilter->SetOperationType(BoolFilterType::AndNot);
   boolFilter->Update();    
@@ -180,79 +230,158 @@ GenerateData() {
   boolFilter->Update(); 
   output = boolFilter->GetOutput();
 
+  // Auto crop to gain support area
   output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue()); 
-  ////autoCropFilter->GetOutput();  typedef clitk::AutoCropFilter<ImageType> CropFilterType;
-  //typename CropFilterType::Pointer cropFilter = CropFilterType::New();
-  //cropFilter->SetInput(output);
-  //cropFilter->Update();   
-  //output = cropFilter->GetOutput();
-
-  this->template StopCurrentStep<TImageType>(output);
+  this->template StopCurrentStep<ImageType>(output);
 
   // Step 2: LR limits from lung (need separate lung ?)
   StartNewStep("Left/Right limits with lungs");
-  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType> RelPosFilterType;
+
+  /*
+  // WE DO NOT NEED THE FOLLOWING ?
+  // Get separate lung images to get only the right and left lung (because RelativePositionPropImageFilter only consider fg=1);
+  // (label must be '1' because right is greater than left).
+  ImagePointer right_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 2, 0);
+  ImagePointer left_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 1, 0);
+  writeImage<ImageType>(right_lung, "right.mhd");
+  writeImage<ImageType>(left_lung, "left.mhd");
+  */
+
+  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType> RelPosFilterType;
   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
   relPosFilter->VerboseStepOff();
   relPosFilter->WriteStepOff();
   relPosFilter->SetInput(output); 
+  //relPosFilter->SetInputObject(left_lung); 
   relPosFilter->SetInputObject(lung); 
-  relPosFilter->SetOrientationType(RelPosFilterType::LeftTo);
+  relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); // warning left lung is at right ;)
   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
   relPosFilter->Update();
   output = relPosFilter->GetOutput();
+  // writeImage<ImageType>(right_lung, "step4-left.mhd");
 
   relPosFilter->SetInput(output); 
   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
   relPosFilter->VerboseStepOff();
   relPosFilter->WriteStepOff();
+  //relPosFilter->SetInputObject(right_lung);
   relPosFilter->SetInputObject(lung); 
   relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
   relPosFilter->Update();   
   output = relPosFilter->GetOutput();
-  this->template StopCurrentStep<TImageType>(output);
+  this->template StopCurrentStep<ImageType>(output);
 
   // Step 3: AP limits from bones
   StartNewStep("Ant/Post limits with bones");
+  ImageConstPointer bones_ant;
+  ImageConstPointer bones_post;
+
+  // Find ant-post coordinate of trachea (assume the carena position is a
+  // good estimation of the ant-post position of the trachea)
+  ImagePointType carina;
+  LoadAFDB();
+  GetAFDB()->GetPoint3D("Carina", carina);
+  DD(carina);
+  ImageIndexType index_trachea;
+  bones->TransformPhysicalPointToIndex(carina, index_trachea);
+  DD(index_trachea);
+  
+  // Split bone image first into two parts (ant and post)
+  typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> ROIFilterType;
+  //  typedef itk::ExtractImageFilter<ImageType, ImageType> ROIFilterType;
+  typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
+  ImageRegionType region = bones->GetLargestPossibleRegion();
+  ImageSizeType size = region.GetSize();
+  DD(size);
+  size[1] =  index_trachea[1]; //size[1]/2.0;
+  DD(size);
+  region.SetSize(size);
+  roiFilter->SetInput(bones);
+  //  roiFilter->SetExtractionRegion(region);
+  roiFilter->SetRegionOfInterest(region);
+  roiFilter->ReleaseDataFlagOff();
+  roiFilter->Update();
+  bones_ant = roiFilter->GetOutput();
+  writeImage<ImageType>(bones_ant, "b_ant.mhd");
+  
+  //  roiFilter->ResetPipeline();// = ROIFilterType::New();  
+  roiFilter = ROIFilterType::New();  
+  ImageIndexType index = region.GetIndex();
+  index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
+  size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
+  DD(size);
+  region.SetIndex(index);
+  region.SetSize(size);
+  roiFilter->SetInput(bones);
+  //  roiFilter->SetExtractionRegion(region);
+  roiFilter->SetRegionOfInterest(region);
+  roiFilter->ReleaseDataFlagOff();
+  roiFilter->Update();
+  bones_post = roiFilter->GetOutput();
+  writeImage<ImageType>(bones_post, "b_post.mhd");
+
+  // Go ! 
   relPosFilter->SetCurrentStepNumber(0);
   relPosFilter->ResetPipeline();// = RelPosFilterType::New();
   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
   relPosFilter->VerboseStepOff();
   relPosFilter->WriteStepOff();
   relPosFilter->SetInput(output); 
-  relPosFilter->SetInputObject(bones); 
+  relPosFilter->SetInputObject(bones_post); 
   relPosFilter->SetOrientationType(RelPosFilterType::AntTo);
   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
   relPosFilter->Update();
+  output = relPosFilter->GetOutput();
+  writeImage<ImageType>(output, "post.mhd");
 
   relPosFilter->SetInput(relPosFilter->GetOutput()); 
   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
   relPosFilter->VerboseStepOff();
   relPosFilter->WriteStepOff();
-  relPosFilter->SetInputObject(bones); 
+  relPosFilter->SetInput(output); 
+  relPosFilter->SetInputObject(bones_ant); 
   relPosFilter->SetOrientationType(RelPosFilterType::PostTo);
   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
   relPosFilter->Update();   
   output = relPosFilter->GetOutput();
-  this->template StopCurrentStep<TImageType>(output);
+  this->template StopCurrentStep<ImageType>(output);
 
   // Get CCL
-  output = clitk::Labelize<TImageType>(output, GetBackgroundValue(), true, 100);
-  // output = RemoveLabels<TImageType>(output, BG, param->GetLabelsToRemove());
-  output = clitk::KeepLabels<TImageType>(output, GetBackgroundValue(), 
-                                  GetForegroundValue(), 1, 1, 0);
+  StartNewStep("Keep main connected component");
+  output = clitk::Labelize<ImageType>(output, GetBackgroundValue(), true, 500);
+  // output = RemoveLabels<ImageType>(output, BG, param->GetLabelsToRemove());
+  output = clitk::KeepLabels<ImageType>(output, GetBackgroundValue(), 
+                                        GetForegroundValue(), 1, 1, 0);
+  this->template StopCurrentStep<ImageType>(output);
+
+  // Step : Lower limits from lung (need separate lung ?)
+  StartNewStep("Lower limits with lungs");
+  relPosFilter = RelPosFilterType::New();
+  relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
+  relPosFilter->VerboseStepOff();
+  relPosFilter->WriteStepOff();
+  relPosFilter->SetInput(output); 
+  DD(output->GetLargestPossibleRegion().GetIndex());
+  //  relPosFilter->SetInputObject(left_lung); 
+  relPosFilter->SetInputObject(lung); 
+  relPosFilter->SetOrientationType(RelPosFilterType::SupTo);
+  relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
+  relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
+  relPosFilter->Update();
+  output = relPosFilter->GetOutput();
+  DD(output->GetLargestPossibleRegion());
 
   output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue()); 
-  //  cropFilter = CropFilterType::New();
-  //cropFilter->SetInput(output);
-  //cropFilter->Update();   
-  //output = cropFilter->GetOutput();
+  //  roiFilter = ROIFilterType::New();
+  //roiFilter->SetInput(output);
+  //roiFilter->Update();   
+  //output = roiFilter->GetOutput();
 
   // Final Step -> set output
   this->SetNthOutput(0, output);