segmentation/clitkExtractMediastinumFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ======================================================================-====*/
  18
  19 #ifndef CLITKEXTRACTMEDIASTINUMFILTER_TXX
  20 #define CLITKEXTRACTMEDIASTINUMFILTER_TXX
  21
  22 // clitk
  23 #include "clitkCommon.h"
  24 #include "clitkExtractMediastinumFilter.h"
  25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
  26 #include "clitkSegmentationUtils.h"
  27 #include "clitkExtractAirwayTreeInfoFilter.h"
  28
  29 // itk
  30 #include <deque>
  31 #include "itkStatisticsLabelMapFilter.h"
  32 #include "itkLabelImageToStatisticsLabelMapFilter.h"
  33 #include "itkRegionOfInterestImageFilter.h"
  34 #include "itkBinaryThresholdImageFilter.h"
  35
  36 // itk ENST
  37 #include "RelativePositionPropImageFilter.h"
  38
  39 //--------------------------------------------------------------------
  40 template <class ImageType>
  41 clitk::ExtractMediastinumFilter<ImageType>::
  42 ExtractMediastinumFilter():
  43   clitk::FilterBase(),
  44   itk::ImageToImageFilter<ImageType, ImageType>()
  45 {
  46   this->SetNumberOfRequiredInputs(4);
  47
  48   SetBackgroundValuePatient(0);
  49   SetBackgroundValueLung(0);
  50   SetBackgroundValueBones(0);
  51   SetForegroundValueLeftLung(1);
  52   SetForegroundValueRightLung(2);
  53   SetBackgroundValue(0);
  54   SetForegroundValue(1);
  55
  56   SetIntermediateSpacing(6);
  57   SetFuzzyThreshold1(0.6);
  58   SetFuzzyThreshold2(0.7);
  59 }
  60 //--------------------------------------------------------------------
  61
  62
  63 //--------------------------------------------------------------------
  64 template <class ImageType>
  65 void
  66 clitk::ExtractMediastinumFilter<ImageType>::
  67 SetInputPatientLabelImage(const ImageType * image, ImagePixelType bg)
  68 {
  69   this->SetNthInput(0, const_cast<ImageType *>(image));
  70   m_BackgroundValuePatient = bg;
  71 }
  72 //--------------------------------------------------------------------
  73
  74
  75 //--------------------------------------------------------------------
  76 template <class ImageType>
  77 void
  78 clitk::ExtractMediastinumFilter<ImageType>::
  79 SetInputLungLabelImage(const ImageType * image, ImagePixelType bg,
  80                        ImagePixelType fgLeft, ImagePixelType fgRight)
  81 {
  82   this->SetNthInput(1, const_cast<ImageType *>(image));
  83   m_BackgroundValueLung = bg;
  84   m_ForegroundValueLeftLung = fgLeft;
  85   m_ForegroundValueRightLung = fgRight;
  86 }
  87 //--------------------------------------------------------------------
  88
  89
  90 //--------------------------------------------------------------------
  91 template <class ImageType>
  92 void
  93 clitk::ExtractMediastinumFilter<ImageType>::
  94 SetInputBonesLabelImage(const ImageType * image, ImagePixelType bg)
  95 {
  96   this->SetNthInput(2, const_cast<ImageType *>(image));
  97   m_BackgroundValueBones = bg;
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 template <class ImageType>
 104 void
 105 clitk::ExtractMediastinumFilter<ImageType>::
 106 SetInputTracheaLabelImage(const ImageType * image, ImagePixelType bg)
 107 {
 108   this->SetNthInput(3, const_cast<ImageType *>(image));
 109   m_BackgroundValueTrachea = bg;
 110 }
 111 //--------------------------------------------------------------------
 112
 113
 114 //--------------------------------------------------------------------
 115 template <class ImageType>
 116 template<class ArgsInfoType>
 117 void
 118 clitk::ExtractMediastinumFilter<ImageType>::
 119 SetArgsInfo(ArgsInfoType mArgsInfo)
 120 {
 121   SetVerboseOption_GGO(mArgsInfo);
 122   SetVerboseStep_GGO(mArgsInfo);
 123   SetWriteStep_GGO(mArgsInfo);
 124   SetVerboseWarningOff_GGO(mArgsInfo);
 125
 126   SetBackgroundValuePatient_GGO(mArgsInfo);
 127   SetBackgroundValueLung_GGO(mArgsInfo);
 128   SetBackgroundValueTrachea_GGO(mArgsInfo);
 129
 130   SetForegroundValueLeftLung_GGO(mArgsInfo);
 131   SetForegroundValueRightLung_GGO(mArgsInfo);
 132
 133   SetIntermediateSpacing_GGO(mArgsInfo);
 134   SetFuzzyThreshold1_GGO(mArgsInfo);
 135   SetFuzzyThreshold2_GGO(mArgsInfo);
 136 }
 137 //--------------------------------------------------------------------
 138
 139
 140 //--------------------------------------------------------------------
 141 template <class ImageType>
 142 void
 143 clitk::ExtractMediastinumFilter<ImageType>::
 144 GenerateOutputInformation() {
 145   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 146   ImagePointer outputImage = this->GetOutput(0);
 147   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
 148 }
 149 //--------------------------------------------------------------------
 150
 151
 152 //--------------------------------------------------------------------
 153 template <class ImageType>
 154 void
 155 clitk::ExtractMediastinumFilter<ImageType>::
 156 GenerateInputRequestedRegion()
 157 {
 158   // Call default
 159   Superclass::GenerateInputRequestedRegion();
 160   // Get input pointers
 161   ImagePointer patient = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 162   ImagePointer lung    = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 163   ImagePointer bones   = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
 164   ImagePointer trachea = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(3));
 165
 166   patient->SetRequestedRegion(patient->GetLargestPossibleRegion());
 167   lung->SetRequestedRegion(lung->GetLargestPossibleRegion());
 168   bones->SetRequestedRegion(bones->GetLargestPossibleRegion());
 169   trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
 170 }
 171 //--------------------------------------------------------------------
 172
 173
 174 //--------------------------------------------------------------------
 175 template <class ImageType>
 176 void
 177 clitk::ExtractMediastinumFilter<ImageType>::
 178 GenerateData()
 179 {
 180   // Get input pointers
 181   ImageConstPointer patient = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 182   ImageConstPointer lung    = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(1));
 183   ImageConstPointer bones   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(2));
 184   ImageConstPointer trachea = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(3));
 185
 186   // Get output pointer
 187   ImagePointer output;
 188
 189   // Step 1: patient minus lungs, minus bones
 190   StartNewStep("Patient contours minus lungs and minus bones");
 191   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 192   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 193   boolFilter->InPlaceOn();
 194   boolFilter->SetInput1(patient);
 195   boolFilter->SetInput2(lung);
 196   boolFilter->SetOperationType(BoolFilterType::AndNot);
 197   boolFilter->Update();
 198   boolFilter->SetInput1(boolFilter->GetOutput());
 199   boolFilter->SetInput2(bones);
 200   boolFilter->SetOperationType(BoolFilterType::AndNot);
 201   boolFilter->Update();
 202   output = boolFilter->GetOutput();
 203
 204   // Auto crop to gain support area
 205   output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue());
 206   this->template StopCurrentStep<ImageType>(output);
 207
 208   // Step 2: LR limits from lung (need separate lung ?)
 209   StartNewStep("Left/Right limits with lungs");
 210
 211   /* // WE DO NOT NEED THE FOLLOWING
 212      // Get separate lung images to get only the right and left lung
 213      // (label must be '1' because right is greater than left).
 214      ImagePointer right_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 2, 0);
 215      ImagePointer left_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 1, 0);
 216      writeImage<ImageType>(right_lung, "right.mhd");
 217      writeImage<ImageType>(left_lung, "left.mhd");
 218   */
 219
 220   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType> RelPosFilterType;
 221   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 222   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 223   relPosFilter->VerboseStepOff();
 224   relPosFilter->WriteStepOff();
 225   relPosFilter->SetInput(output);
 226   DD(output->GetLargestPossibleRegion().GetIndex());
 227   //  relPosFilter->SetInputObject(left_lung);
 228   relPosFilter->SetInputObject(lung);
 229   relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); // warning left lung is at right ;)
 230   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 231   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 232   relPosFilter->Update();
 233   output = relPosFilter->GetOutput();
 234   DD(output->GetLargestPossibleRegion());
 235
 236   relPosFilter->SetInput(output);
 237   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 238   relPosFilter->VerboseStepOff();
 239   relPosFilter->WriteStepOff();
 240   //  relPosFilter->SetInputObject(right_lung);
 241   relPosFilter->SetInputObject(lung);
 242   relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
 243   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 244   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 245   relPosFilter->Update();
 246   output = relPosFilter->GetOutput();
 247   DD(output->GetLargestPossibleRegion());
 248   this->template StopCurrentStep<ImageType>(output);
 249
 250   // Step 3: AP limits from bones
 251   StartNewStep("Ant/Post limits with bones");
 252   ImageConstPointer bones_ant;
 253   ImageConstPointer bones_post;
 254
 255   // Find ant-post coordinate of trachea (assume the first point is a
 256   // good estimation of the ant-post position of the trachea)
 257   typedef clitk::ExtractAirwayTreeInfoFilter<ImageType> AirwayFilter;
 258   typename AirwayFilter::Pointer airwayfilter = AirwayFilter::New();
 259   airwayfilter->SetInput(trachea);
 260   airwayfilter->Update();
 261   DD(airwayfilter->GetFirstTracheaPoint());
 262   ImagePointType point_trachea = airwayfilter->GetFirstTracheaPoint();
 263   ImageIndexType index_trachea;
 264   bones->TransformPhysicalPointToIndex(point_trachea, index_trachea);
 265   DD(index_trachea);
 266
 267   // Split bone image first into two parts (ant and post)
 268   typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 269   //  typedef itk::ExtractImageFilter<ImageType, ImageType> CropFilterType;
 270   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 271   ImageRegionType region = bones->GetLargestPossibleRegion();
 272   ImageSizeType size = region.GetSize();
 273   DD(size);
 274   size[1] =  index_trachea[1]; //size[1]/2.0;
 275   DD(size);
 276   region.SetSize(size);
 277   cropFilter->SetInput(bones);
 278   //  cropFilter->SetExtractionRegion(region);
 279   cropFilter->SetRegionOfInterest(region);
 280   cropFilter->ReleaseDataFlagOff();
 281   cropFilter->Update();
 282   bones_ant = cropFilter->GetOutput();
 283   writeImage<ImageType>(bones_ant, "b_ant.mhd");
 284
 285   //  cropFilter->ResetPipeline();// = CropFilterType::New();
 286   cropFilter = CropFilterType::New();
 287   ImageIndexType index = region.GetIndex();
 288   size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
 289   index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
 290   DD(size);
 291   region.SetIndex(index);
 292   region.SetSize(size);
 293   cropFilter->SetInput(bones);
 294   //  cropFilter->SetExtractionRegion(region);
 295   cropFilter->SetRegionOfInterest(region);
 296   cropFilter->ReleaseDataFlagOff();
 297   cropFilter->Update();
 298   bones_post = cropFilter->GetOutput();
 299   writeImage<ImageType>(bones_post, "b_post.mhd");
 300
 301   // Go !
 302   relPosFilter->SetCurrentStepNumber(0);
 303   relPosFilter->ResetPipeline();// = RelPosFilterType::New();
 304   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 305   relPosFilter->VerboseStepOff();
 306   relPosFilter->WriteStepOff();
 307   relPosFilter->SetInput(output);
 308   relPosFilter->SetInputObject(bones_post);
 309   relPosFilter->SetOrientationType(RelPosFilterType::AntTo);
 310   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 311   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 312   relPosFilter->Update();
 313   output = relPosFilter->GetOutput();
 314   writeImage<ImageType>(output, "post.mhd");
 315
 316   relPosFilter->SetInput(relPosFilter->GetOutput());
 317   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 318   relPosFilter->VerboseStepOff();
 319   relPosFilter->WriteStepOff();
 320   relPosFilter->SetInput(output);
 321   relPosFilter->SetInputObject(bones_ant);
 322   relPosFilter->SetOrientationType(RelPosFilterType::PostTo);
 323   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 324   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 325   relPosFilter->Update();
 326   output = relPosFilter->GetOutput();
 327   this->template StopCurrentStep<ImageType>(output);
 328   // Get CCL
 329   output = clitk::Labelize<ImageType>(output, GetBackgroundValue(), true, 100);
 330   // output = RemoveLabels<ImageType>(output, BG, param->GetLabelsToRemove());
 331   output = clitk::KeepLabels<ImageType>(output, GetBackgroundValue(),
 332                                         GetForegroundValue(), 1, 1, 0);
 333
 334   output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue());
 335   //  cropFilter = CropFilterType::New();
 336   //cropFilter->SetInput(output);
 337   //cropFilter->Update();
 338   //output = cropFilter->GetOutput();
 339
 340   // Final Step -> set output
 341   this->SetNthOutput(0, output);
 342 }
 343 //--------------------------------------------------------------------
 344
 345
 346 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX