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