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