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 "clitkExtractAirwaysTreeInfoFilter.h"
  28 #include "clitkCropLikeImageFilter.h"
  29
  30 // std
  31 #include <deque>
  32
  33 // itk
  34 #include "itkStatisticsLabelMapFilter.h"
  35 #include "itkLabelImageToStatisticsLabelMapFilter.h"
  36 #include "itkRegionOfInterestImageFilter.h"
  37 #include "itkBinaryThresholdImageFilter.h"
  38
  39 // itk ENST
  40 #include "RelativePositionPropImageFilter.h"
  41
  42 //--------------------------------------------------------------------
  43 template <class ImageType>
  44 clitk::ExtractMediastinumFilter<ImageType>::
  45 ExtractMediastinumFilter():
  46   clitk::FilterBase(),
  47   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  48   itk::ImageToImageFilter<ImageType, ImageType>()
  49 {
  50   this->SetNumberOfRequiredInputs(4);
  51
  52   SetBackgroundValuePatient(0);
  53   SetBackgroundValueLung(0);
  54   SetBackgroundValueBones(0);
  55   SetForegroundValueLeftLung(1);
  56   SetForegroundValueRightLung(2);
  57   SetBackgroundValue(0);
  58   SetForegroundValue(1);
  59
  60   SetIntermediateSpacing(6);
  61   SetFuzzyThreshold1(0.4);
  62   SetFuzzyThreshold2(0.6);
  63   SetFuzzyThreshold3(0.2);
  64 }
  65 //--------------------------------------------------------------------
  66
  67
  68 //--------------------------------------------------------------------
  69 template <class ImageType>
  70 void
  71 clitk::ExtractMediastinumFilter<ImageType>::
  72 SetInputPatientLabelImage(const ImageType * image, ImagePixelType bg)
  73 {
  74   this->SetNthInput(0, const_cast<ImageType *>(image));
  75   m_BackgroundValuePatient = bg;
  76 }
  77 //--------------------------------------------------------------------
  78
  79
  80 //--------------------------------------------------------------------
  81 template <class ImageType>
  82 void
  83 clitk::ExtractMediastinumFilter<ImageType>::
  84 SetInputLungLabelImage(const ImageType * image, ImagePixelType bg,
  85                        ImagePixelType fgLeft, ImagePixelType fgRight)
  86 {
  87   this->SetNthInput(1, const_cast<ImageType *>(image));
  88   m_BackgroundValueLung = bg;
  89   m_ForegroundValueLeftLung = fgLeft;
  90   m_ForegroundValueRightLung = fgRight;
  91 }
  92 //--------------------------------------------------------------------
  93
  94
  95 //--------------------------------------------------------------------
  96 template <class ImageType>
  97 void
  98 clitk::ExtractMediastinumFilter<ImageType>::
  99 SetInputBonesLabelImage(const ImageType * image, ImagePixelType bg)
 100 {
 101   this->SetNthInput(2, const_cast<ImageType *>(image));
 102   m_BackgroundValueBones = bg;
 103 }
 104 //--------------------------------------------------------------------
 105
 106
 107 //--------------------------------------------------------------------
 108 template <class ImageType>
 109 void
 110 clitk::ExtractMediastinumFilter<ImageType>::
 111 SetInputTracheaLabelImage(const ImageType * image, ImagePixelType bg)
 112 {
 113   this->SetNthInput(3, const_cast<ImageType *>(image));
 114   m_BackgroundValueTrachea = bg;
 115 }
 116 //--------------------------------------------------------------------
 117
 118
 119 //--------------------------------------------------------------------
 120 template <class ImageType>
 121 template<class ArgsInfoType>
 122 void
 123 clitk::ExtractMediastinumFilter<ImageType>::
 124 SetArgsInfo(ArgsInfoType mArgsInfo)
 125 {
 126   SetVerboseOption_GGO(mArgsInfo);
 127   SetVerboseStep_GGO(mArgsInfo);
 128   SetWriteStep_GGO(mArgsInfo);
 129   SetVerboseWarningOff_GGO(mArgsInfo);
 130
 131   SetBackgroundValuePatient_GGO(mArgsInfo);
 132   SetBackgroundValueLung_GGO(mArgsInfo);
 133   SetBackgroundValueTrachea_GGO(mArgsInfo);
 134
 135   SetForegroundValueLeftLung_GGO(mArgsInfo);
 136   SetForegroundValueRightLung_GGO(mArgsInfo);
 137
 138   SetIntermediateSpacing_GGO(mArgsInfo);
 139   SetFuzzyThreshold1_GGO(mArgsInfo);
 140   SetFuzzyThreshold2_GGO(mArgsInfo);
 141   SetFuzzyThreshold3_GGO(mArgsInfo);
 142
 143   SetAFDBFilename_GGO(mArgsInfo);
 144 }
 145 //--------------------------------------------------------------------
 146
 147
 148 //--------------------------------------------------------------------
 149 template <class ImageType>
 150 void
 151 clitk::ExtractMediastinumFilter<ImageType>::
 152 GenerateOutputInformation() {
 153   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 154   ImagePointer outputImage = this->GetOutput(0);
 155   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
 156 }
 157 //--------------------------------------------------------------------
 158
 159
 160 //--------------------------------------------------------------------
 161 template <class ImageType>
 162 void
 163 clitk::ExtractMediastinumFilter<ImageType>::
 164 GenerateInputRequestedRegion()
 165 {
 166   // Call default
 167   Superclass::GenerateInputRequestedRegion();
 168
 169   // Get input pointers
 170   LoadAFDB();
 171   ImagePointer patient = GetAFDB()->template GetImage <ImageType>("patient");
 172   ImagePointer lung = GetAFDB()->template GetImage <ImageType>("lungs");
 173   ImagePointer bones = GetAFDB()->template GetImage <ImageType>("bones");
 174   ImagePointer trachea = GetAFDB()->template GetImage <ImageType>("trachea");
 175
 176   patient->SetRequestedRegion(patient->GetLargestPossibleRegion());
 177   lung->SetRequestedRegion(lung->GetLargestPossibleRegion());
 178   bones->SetRequestedRegion(bones->GetLargestPossibleRegion());
 179   trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
 180 }
 181 //--------------------------------------------------------------------
 182
 183
 184 //--------------------------------------------------------------------
 185 template <class ImageType>
 186 void
 187 clitk::ExtractMediastinumFilter<ImageType>::
 188 GenerateData()
 189 {
 190   // Get input pointers
 191   ImagePointer patient = GetAFDB()->template GetImage <ImageType>("patient");
 192   ImagePointer lung = GetAFDB()->template GetImage <ImageType>("lungs");
 193   ImagePointer bones = GetAFDB()->template GetImage <ImageType>("bones");
 194   ImagePointer trachea = GetAFDB()->template GetImage <ImageType>("trachea");
 195
 196   // Get output pointer
 197   ImagePointer output;
 198
 199   // Step 0: Crop support (patient) to lung extend in RL
 200   StartNewStep("Crop support like lungs along LR");
 201   typedef clitk::CropLikeImageFilter<ImageType> CropFilterType;
 202   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 203   cropFilter->SetInput(patient);
 204   cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
 205   cropFilter->Update();
 206   output = cropFilter->GetOutput();
 207   this->template StopCurrentStep<ImageType>(output);
 208
 209   // Step 0: Crop support (previous) to bones extend in AP
 210   StartNewStep("Crop support like bones along AP");
 211   cropFilter = CropFilterType::New();
 212   cropFilter->SetInput(output);
 213   cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
 214   cropFilter->Update();
 215   output = cropFilter->GetOutput();
 216   this->template StopCurrentStep<ImageType>(output);
 217
 218   // Step 1: patient minus lungs, minus bones
 219   StartNewStep("Patient contours minus lungs and minus bones");
 220   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 221   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 222   boolFilter->InPlaceOn();
 223   boolFilter->SetInput1(output);
 224   boolFilter->SetInput2(lung);
 225   boolFilter->SetOperationType(BoolFilterType::AndNot);
 226   boolFilter->Update();
 227   boolFilter->SetInput1(boolFilter->GetOutput());
 228   boolFilter->SetInput2(bones);
 229   boolFilter->SetOperationType(BoolFilterType::AndNot);
 230   boolFilter->Update();
 231   output = boolFilter->GetOutput();
 232
 233   // Auto crop to gain support area
 234   output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue());
 235   this->template StopCurrentStep<ImageType>(output);
 236
 237   // Step 2: LR limits from lung (need separate lung ?)
 238   StartNewStep("Left/Right limits with lungs");
 239
 240   /*
 241   // WE DO NOT NEED THE FOLLOWING ?
 242   // Get separate lung images to get only the right and left lung (because RelativePositionPropImageFilter only consider fg=1);
 243   // (label must be '1' because right is greater than left).
 244   ImagePointer right_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 2, 0);
 245   ImagePointer left_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 1, 0);
 246   writeImage<ImageType>(right_lung, "right.mhd");
 247   writeImage<ImageType>(left_lung, "left.mhd");
 248   */
 249
 250   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType> RelPosFilterType;
 251   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 252   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 253   relPosFilter->VerboseStepOff();
 254   relPosFilter->WriteStepOff();
 255   relPosFilter->SetInput(output);
 256   //relPosFilter->SetInputObject(left_lung);
 257   relPosFilter->SetInputObject(lung);
 258   relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); // warning left lung is at right ;)
 259   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 260   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 261   relPosFilter->Update();
 262   output = relPosFilter->GetOutput();
 263   // writeImage<ImageType>(right_lung, "step4-left.mhd");
 264
 265   relPosFilter->SetInput(output);
 266   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 267   relPosFilter->VerboseStepOff();
 268   relPosFilter->WriteStepOff();
 269   //relPosFilter->SetInputObject(right_lung);
 270   relPosFilter->SetInputObject(lung);
 271   relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
 272   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 273   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 274   relPosFilter->Update();
 275   output = relPosFilter->GetOutput();
 276   this->template StopCurrentStep<ImageType>(output);
 277
 278   // Step 3: AP limits from bones
 279   StartNewStep("Ant/Post limits with bones");
 280   ImageConstPointer bones_ant;
 281   ImageConstPointer bones_post;
 282
 283   // Find ant-post coordinate of trachea (assume the carena position is a
 284   // good estimation of the ant-post position of the trachea)
 285   ImagePointType carina;
 286   LoadAFDB();
 287   GetAFDB()->GetPoint3D("Carina", carina);
 288   DD(carina);
 289   ImageIndexType index_trachea;
 290   bones->TransformPhysicalPointToIndex(carina, index_trachea);
 291   DD(index_trachea);
 292
 293   // Split bone image first into two parts (ant and post)
 294   typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> ROIFilterType;
 295   //  typedef itk::ExtractImageFilter<ImageType, ImageType> ROIFilterType;
 296   typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
 297   ImageRegionType region = bones->GetLargestPossibleRegion();
 298   ImageSizeType size = region.GetSize();
 299   DD(size);
 300   size[1] =  index_trachea[1]; //size[1]/2.0;
 301   DD(size);
 302   region.SetSize(size);
 303   roiFilter->SetInput(bones);
 304   //  roiFilter->SetExtractionRegion(region);
 305   roiFilter->SetRegionOfInterest(region);
 306   roiFilter->ReleaseDataFlagOff();
 307   roiFilter->Update();
 308   bones_ant = roiFilter->GetOutput();
 309   writeImage<ImageType>(bones_ant, "b_ant.mhd");
 310
 311   //  roiFilter->ResetPipeline();// = ROIFilterType::New();
 312   roiFilter = ROIFilterType::New();
 313   ImageIndexType index = region.GetIndex();
 314   index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
 315   size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
 316   DD(size);
 317   region.SetIndex(index);
 318   region.SetSize(size);
 319   roiFilter->SetInput(bones);
 320   //  roiFilter->SetExtractionRegion(region);
 321   roiFilter->SetRegionOfInterest(region);
 322   roiFilter->ReleaseDataFlagOff();
 323   roiFilter->Update();
 324   bones_post = roiFilter->GetOutput();
 325   writeImage<ImageType>(bones_post, "b_post.mhd");
 326
 327   // Go !
 328   relPosFilter->SetCurrentStepNumber(0);
 329   relPosFilter->ResetPipeline();// = RelPosFilterType::New();
 330   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 331   relPosFilter->VerboseStepOff();
 332   relPosFilter->WriteStepOff();
 333   relPosFilter->SetInput(output);
 334   relPosFilter->SetInputObject(bones_post);
 335   relPosFilter->SetOrientationType(RelPosFilterType::AntTo);
 336   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 337   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 338   relPosFilter->Update();
 339   output = relPosFilter->GetOutput();
 340   writeImage<ImageType>(output, "post.mhd");
 341
 342   relPosFilter->SetInput(relPosFilter->GetOutput());
 343   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 344   relPosFilter->VerboseStepOff();
 345   relPosFilter->WriteStepOff();
 346   relPosFilter->SetInput(output);
 347   relPosFilter->SetInputObject(bones_ant);
 348   relPosFilter->SetOrientationType(RelPosFilterType::PostTo);
 349   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 350   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 351   relPosFilter->Update();
 352   output = relPosFilter->GetOutput();
 353   this->template StopCurrentStep<ImageType>(output);
 354
 355   // Get CCL
 356   StartNewStep("Keep main connected component");
 357   output = clitk::Labelize<ImageType>(output, GetBackgroundValue(), true, 500);
 358   // output = RemoveLabels<ImageType>(output, BG, param->GetLabelsToRemove());
 359   output = clitk::KeepLabels<ImageType>(output, GetBackgroundValue(),
 360                                         GetForegroundValue(), 1, 1, 0);
 361   this->template StopCurrentStep<ImageType>(output);
 362
 363   // Step : Lower limits from lung (need separate lung ?)
 364   StartNewStep("Lower limits with lungs");
 365   relPosFilter = RelPosFilterType::New();
 366   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 367   relPosFilter->VerboseStepOff();
 368   relPosFilter->WriteStepOff();
 369   relPosFilter->SetInput(output);
 370   DD(output->GetLargestPossibleRegion().GetIndex());
 371   //  relPosFilter->SetInputObject(left_lung);
 372   relPosFilter->SetInputObject(lung);
 373   relPosFilter->SetOrientationType(RelPosFilterType::SupTo);
 374   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 375   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
 376   relPosFilter->Update();
 377   output = relPosFilter->GetOutput();
 378   DD(output->GetLargestPossibleRegion());
 379
 380   output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue());
 381   //  roiFilter = ROIFilterType::New();
 382   //roiFilter->SetInput(output);
 383   //roiFilter->Update();
 384   //output = roiFilter->GetOutput();
 385
 386   // Final Step -> set output
 387   this->SetNthOutput(0, output);
 388 }
 389 //--------------------------------------------------------------------
 390
 391
 392 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX