segmentation/clitkExtractLymphStationsFilter.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 CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
  20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkCommon.h"
  24 #include "clitkExtractLymphStationsFilter.h"
  25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
  26 #include "clitkSegmentationUtils.h"
  27 #include "clitkAutoCropFilter.h"
  28 #include "clitkSegmentationUtils.h"
  29 #include "clitkSliceBySliceRelativePositionFilter.h"
  30
  31 // itk
  32 #include <itkStatisticsLabelMapFilter.h>
  33 #include <itkLabelImageToStatisticsLabelMapFilter.h>
  34 #include <itkRegionOfInterestImageFilter.h>
  35 #include <itkBinaryThresholdImageFilter.h>
  36 #include <itkImageSliceConstIteratorWithIndex.h>
  37 #include <itkBinaryThinningImageFilter.h>
  38
  39 // itk ENST
  40 #include "RelativePositionPropImageFilter.h"
  41
  42 //--------------------------------------------------------------------
  43 template <class TImageType>
  44 clitk::ExtractLymphStationsFilter<TImageType>::
  45 ExtractLymphStationsFilter():
  46   clitk::FilterBase(),
  47   itk::ImageToImageFilter<TImageType, TImageType>()
  48 {
  49   this->SetNumberOfRequiredInputs(1);
  50   SetBackgroundValueMediastinum(0);
  51   SetBackgroundValueTrachea(0);
  52   SetBackgroundValue(0);
  53   SetForegroundValue(1);
  54
  55   SetIntermediateSpacing(6);
  56   SetFuzzyThreshold1(0.6);
  57 }
  58 //--------------------------------------------------------------------
  59
  60
  61 //--------------------------------------------------------------------
  62 template <class TImageType>
  63 void
  64 clitk::ExtractLymphStationsFilter<TImageType>::
  65 SetInputMediastinumLabelImage(const TImageType * image, ImagePixelType bg) {
  66   this->SetNthInput(0, const_cast<TImageType *>(image));
  67   m_BackgroundValueMediastinum = bg;
  68   SetCarenaZPositionInMM(image->GetOrigin()[2]+image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]);
  69   SetMiddleLobeBronchusZPositionInMM(image->GetOrigin()[2]);
  70   // DD(m_CarenaZPositionInMM);
  71 //   DD(m_MiddleLobeBronchusZPositionInMM);
  72 }
  73 //--------------------------------------------------------------------
  74
  75
  76 //--------------------------------------------------------------------
  77 template <class TImageType>
  78 void
  79 clitk::ExtractLymphStationsFilter<TImageType>::
  80 SetInputTracheaLabelImage(const TImageType * image, ImagePixelType bg) {
  81   this->SetNthInput(1, const_cast<TImageType *>(image));
  82   m_BackgroundValueTrachea = bg;
  83 }
  84 //--------------------------------------------------------------------
  85
  86
  87 //--------------------------------------------------------------------
  88 template <class TImageType>
  89 template<class ArgsInfoType>
  90 void
  91 clitk::ExtractLymphStationsFilter<TImageType>::
  92 SetArgsInfo(ArgsInfoType mArgsInfo)
  93 {
  94   SetVerboseOption_GGO(mArgsInfo);
  95   SetVerboseStep_GGO(mArgsInfo);
  96   SetWriteStep_GGO(mArgsInfo);
  97   SetVerboseWarningOff_GGO(mArgsInfo);
  98   SetCarenaZPositionInMM_GGO(mArgsInfo);
  99   SetMiddleLobeBronchusZPositionInMM_GGO(mArgsInfo);
 100   SetIntermediateSpacing_GGO(mArgsInfo);
 101   SetFuzzyThreshold1_GGO(mArgsInfo);
 102   //SetBackgroundValueMediastinum_GGO(mArgsInfo);
 103 }
 104 //--------------------------------------------------------------------
 105
 106
 107 //--------------------------------------------------------------------
 108 template <class TImageType>
 109 void
 110 clitk::ExtractLymphStationsFilter<TImageType>::
 111 GenerateOutputInformation() {
 112   //  Superclass::GenerateOutputInformation();
 113
 114   // Get input
 115   m_mediastinum = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(0));
 116   m_trachea = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(1));
 117
 118   // ----------------------------------------------------------------
 119   // ----------------------------------------------------------------
 120   // Superior limit = carena
 121   // Inferior limit = origine middle lobe bronchus
 122   StartNewStep("Inf/Sup mediastinum limits with carena/bronchus");
 123   ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
 124   ImageSizeType size = region.GetSize();
 125   ImageIndexType index = region.GetIndex();
 126   DD(m_CarenaZPositionInMM);
 127   DD(m_MiddleLobeBronchusZPositionInMM);
 128   index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
 129   size[2] = ceil((m_CarenaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]);
 130   region.SetSize(size);
 131   region.SetIndex(index);
 132   DD(region);
 133   typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 134   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 135   cropFilter->SetInput(m_mediastinum);
 136   cropFilter->SetRegionOfInterest(region);
 137   cropFilter->Update();
 138   m_working_image = cropFilter->GetOutput();
 139   // Auto Crop (because following rel pos is faster)
 140   m_working_image = clitk::AutoCrop<ImageType>(m_working_image, 0);
 141   StopCurrentStep<ImageType>(m_working_image);
 142   m_output = m_working_image;
 143
 144   // ----------------------------------------------------------------
 145   // ----------------------------------------------------------------
 146   // Separate trachea in two CCL
 147   StartNewStep("Separate trachea under carena");
 148   // DD(region);
 149   ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion();
 150   for(int i=0; i<3; i++) {
 151     index[i] = floor(((index[i]*m_mediastinum->GetSpacing()[i])+m_mediastinum->GetOrigin()[i]
 152                       -m_trachea->GetOrigin()[i])/m_trachea->GetSpacing()[i]);
 153     // DD(index[i]);
 154     size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]);
 155     //  DD(size[i]);
 156     if (index[i] < 0) {
 157       size[i] += index[i];
 158       index[i] = 0;
 159     }
 160     if (size[i]+index[i] > (trachea_region.GetSize()[i] + trachea_region.GetIndex()[i])) {
 161       size[i] = trachea_region.GetSize()[i] + trachea_region.GetIndex()[i] - index[i];
 162     }
 163   }
 164   // DD(index);
 165   //   DD(size);
 166   region.SetIndex(index);
 167   region.SetSize(size);
 168   //  typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 169   //  typename CropFilterType::Pointer
 170   cropFilter = CropFilterType::New();
 171  //  m_trachea.Print(std::cout);
 172   cropFilter->SetInput(m_trachea);
 173   cropFilter->SetRegionOfInterest(region);
 174   cropFilter->Update();
 175   m_working_trachea = cropFilter->GetOutput();
 176
 177   // Labelize and consider two main labels
 178   m_working_trachea = Labelize<ImageType>(m_working_trachea, 0, true, 1);
 179
 180   // Detect wich label is at Left
 181   typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
 182   SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
 183   iter.SetFirstDirection(0);
 184   iter.SetSecondDirection(1);
 185   iter.GoToBegin();
 186   bool stop = false;
 187   ImagePixelType leftLabel;
 188   ImagePixelType rightLabel;
 189   while (!stop) {
 190     if (iter.Get() != m_BackgroundValueTrachea) {
 191       //     DD(iter.GetIndex());
 192       //       DD((int)iter.Get());
 193       leftLabel = iter.Get();
 194       stop = true;
 195     }
 196     ++iter;
 197   }
 198   if (leftLabel == 1) rightLabel = 2;
 199   else rightLabel = 1;
 200   DD((int)leftLabel);
 201   DD((int)rightLabel);
 202
 203   // End step
 204   StopCurrentStep<ImageType>(m_working_trachea);
 205
 206   //-----------------------------------------------------
 207   /*  DD("TEST Skeleton");
 208   typedef itk::BinaryThinningImageFilter<ImageType, ImageType> SkeletonFilterType;
 209   typename SkeletonFilterType::Pointer skeletonFilter = SkeletonFilterType::New();
 210   skeletonFilter->SetInput(m_working_trachea);
 211   skeletonFilter->Update();
 212   writeImage<ImageType>(skeletonFilter->GetOutput(), "skel.mhd");
 213   writeImage<ImageType>(skeletonFilter->GetThinning(), "skel2.mhd");
 214   */
 215
 216   //-----------------------------------------------------
 217   StartNewStep("Left limits with bronchus (slice by slice)");
 218   // Select LeftLabel (set right label to 0)
 219   ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, rightLabel, 0);
 220   writeImage<ImageType>(temp, "temp1.mhd");
 221
 222   typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
 223   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 224   sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 225   sliceRelPosFilter->VerboseStepOn();
 226   sliceRelPosFilter->WriteStepOn();
 227   sliceRelPosFilter->SetInput(m_working_image);
 228   sliceRelPosFilter->SetInputObject(temp);
 229   sliceRelPosFilter->SetDirection(2);
 230   sliceRelPosFilter->SetFuzzyThreshold(0.5);
 231   sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::RightTo);
 232   sliceRelPosFilter->Update();
 233   m_working_image = sliceRelPosFilter->GetOutput();
 234   writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
 235
 236
 237   //-----------------------------------------------------
 238   StartNewStep("Right limits with bronchus (slice by slice)");
 239   // Select LeftLabel (set right label to 0)
 240   temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, leftLabel, 0);
 241   writeImage<ImageType>(temp, "temp2.mhd");
 242
 243   sliceRelPosFilter = SliceRelPosFilterType::New();
 244   sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 245   sliceRelPosFilter->VerboseStepOn();
 246   sliceRelPosFilter->WriteStepOn();
 247   sliceRelPosFilter->SetInput(m_working_image);
 248   sliceRelPosFilter->SetInputObject(temp);
 249   sliceRelPosFilter->SetDirection(2);
 250   sliceRelPosFilter->SetFuzzyThreshold(0.5);
 251   sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::LeftTo);
 252   sliceRelPosFilter->Update();
 253   m_working_image = sliceRelPosFilter->GetOutput();
 254   writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
 255
 256
 257   DD("end");
 258   m_output = m_working_image;
 259   StopCurrentStep<ImageType>(m_output);
 260
 261   // Set output image information (required)
 262   ImagePointer outputImage = this->GetOutput(0);
 263   outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
 264   outputImage->SetOrigin(m_working_image->GetOrigin());
 265   outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
 266   DD("end2");
 267 }
 268 //--------------------------------------------------------------------
 269
 270
 271 //--------------------------------------------------------------------
 272 template <class TImageType>
 273 void
 274 clitk::ExtractLymphStationsFilter<TImageType>::
 275 GenerateInputRequestedRegion() {
 276   DD("GenerateInputRequestedRegion");
 277   // Call default
 278   Superclass::GenerateInputRequestedRegion();
 279   // Following needed because output region can be greater than input (trachea)
 280   ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
 281   ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
 282   mediastinum->SetRequestedRegion(mediastinum->GetLargestPossibleRegion());
 283   trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
 284 }
 285 //--------------------------------------------------------------------
 286
 287
 288 //--------------------------------------------------------------------
 289 template <class TImageType>
 290 void
 291 clitk::ExtractLymphStationsFilter<TImageType>::
 292 GenerateData() {
 293   DD("GenerateData");
 294   // Final Step -> graft output (if SetNthOutput => redo)
 295   this->GraftOutput(m_output);
 296 }
 297 //--------------------------------------------------------------------
 298
 299
 300 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX