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 <deque>
  33 #include <itkStatisticsLabelMapFilter.h>
  34 #include <itkLabelImageToStatisticsLabelMapFilter.h>
  35 #include <itkRegionOfInterestImageFilter.h>
  36 #include <itkBinaryThresholdImageFilter.h>
  37 #include <itkImageSliceConstIteratorWithIndex.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 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");
 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   StopCurrentStep<ImageType>(m_working_trachea);
 180
 181   // Detect wich label is at Left
 182   typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
 183   SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
 184   iter.SetFirstDirection(0);
 185   iter.SetSecondDirection(1);
 186   iter.GoToBegin();
 187   bool stop = false;
 188   ImagePixelType leftLabel;
 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   DD((int)leftLabel);
 199
 200   /*
 201   // Relative position
 202   StartNewStep("Left/Right limits with trachea");
 203
 204   // Select LeftLabel (set label 1 to 0)
 205   ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 1, 0);
 206   writeImage<ImageType>(temp, "temp1.mhd");
 207
 208   // Left relative position
 209   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType> RelPosFilterType;
 210   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 211   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 212   relPosFilter->VerboseStepOff();
 213   relPosFilter->WriteStepOff();
 214   relPosFilter->SetInput(m_working_image);
 215   relPosFilter->SetInputObject(temp);
 216   relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
 217   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 218   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 219   relPosFilter->Update();
 220   m_working_image = relPosFilter->GetOutput();
 221
 222   // Select RightLabel (set label 2 to 0)
 223   temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 2, 0);
 224   writeImage<ImageType>(temp, "temp2.mhd");
 225
 226   // Left relative position
 227   relPosFilter = RelPosFilterType::New();
 228   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 229   relPosFilter->VerboseStepOn();
 230   relPosFilter->WriteStepOn();
 231   relPosFilter->SetInput(m_working_image);
 232   relPosFilter->SetInputObject(temp);
 233   relPosFilter->SetOrientationType(RelPosFilterType::LeftTo);
 234   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 235   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 236   relPosFilter->Update();
 237   m_working_image = relPosFilter->GetOutput();
 238   */
 239
 240   //-----------------------------------------------------
 241   StartNewStep("Left/Right limits with trachea (slice by slice");
 242
 243   // Select LeftLabel (set label 1 to 0)
 244   ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 1, 0);
 245   writeImage<ImageType>(temp, "temp1.mhd");
 246
 247   /*
 248     - écrire utilisation de filter SliceBySliceRelPosFilter (combine in a 3D)
 249     - filter SliceBySliceRelPosFilter + combine in a 3D
 250           - resampling, affine transfo to pair the slices
 251           - extract list of images (ecrire utilisation de ExtractSliceFilter)
 252    */
 253
 254   typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
 255   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 256   sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 257   sliceRelPosFilter->VerboseStepOn();
 258   sliceRelPosFilter->WriteStepOn();
 259   sliceRelPosFilter->SetInput(m_working_image);
 260   sliceRelPosFilter->SetInputObject(temp);
 261   sliceRelPosFilter->SetDirection(2);
 262   sliceRelPosFilter->Update();
 263   m_working_image = sliceRelPosFilter->GetOutput();
 264   writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
 265
 266   /*
 267
 268
 269   itk::ImageSliceConstIteratorWithIndex<ImageType> it(temp, temp->GetRequestedRegion());
 270   itk::ImageRegionIterator<SliceType> * ot;
 271   typename SliceType::Pointer slice;
 272   it.SetFirstDirection(0);
 273   it.SetSecondDirection(1);
 274   it.GoToBegin();
 275   typename SliceType::RegionType mSliceRegion;
 276   typename SliceType::IndexType mSliceIndex;
 277   typename SliceType::SizeType mSliceSize;
 278   typename SliceType::SpacingType mSliceSpacing;
 279   typename SliceType::PointType mSliceOrigin;
 280   mSliceIndex[0] = temp->GetLargestPossibleRegion().GetIndex()[0];
 281   mSliceIndex[1] = temp->GetLargestPossibleRegion().GetIndex()[1];
 282   mSliceSize[0] = temp->GetLargestPossibleRegion().GetSize()[0];
 283   mSliceSize[1] = temp->GetLargestPossibleRegion().GetSize()[1];
 284   mSliceSpacing[0] = temp->GetSpacing()[0];
 285   mSliceSpacing[1] = temp->GetSpacing()[1];
 286   mSliceOrigin[0] = temp->GetOrigin()[0];
 287   mSliceOrigin[1] = temp->GetOrigin()[1];
 288   mSliceRegion.SetIndex(mSliceIndex);
 289   mSliceRegion.SetSize(mSliceSize);
 290   int i=0;
 291   while( !it.IsAtEnd() ) {
 292     // DD(i);
 293     slice = SliceType::New();
 294     slice->SetRegions(mSliceRegion);
 295     slice->SetOrigin(mSliceOrigin);
 296     slice->SetSpacing(mSliceSpacing);
 297     slice->Allocate();
 298     ot = new itk::ImageRegionIterator<SliceType>(slice, slice->GetLargestPossibleRegion());
 299     ot->GoToBegin();
 300     // DD(it.GetIndex());
 301     while(!it.IsAtEndOfSlice()) {
 302       // DD(ot->GetIndex());
 303       while(!it.IsAtEndOfLine()) {
 304         ot->Set(it.Get());
 305         ++it;
 306         ++(*ot);
 307       }
 308       it.NextLine();
 309     }
 310     mImageSlices.push_back(slice);
 311     it.NextSlice();
 312     ++i;
 313   }
 314   writeImage<SliceType>(mImageSlices[10], "slice.mhd");
 315
 316   // Perform RelPos by slice
 317   for(int i=0; i<mImageSlices.size(); i++) {
 318     DD(i);
 319     typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
 320     typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 321     //    relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 322     relPosFilter->VerboseStepOff();
 323     relPosFilter->WriteStepOff();
 324     relPosFilter->SetInput(m_working_image);
 325     relPosFilter->SetInputObject(temp);
 326     relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
 327     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 328     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 329     relPosFilter->Update();
 330     m_working_image = relPosFilter->GetOutput();
 331   }
 332
 333   */
 334
 335   DD("end");
 336   m_output = m_working_image;
 337   StopCurrentStep<ImageType>(m_output);
 338
 339   // Set output image information (required)
 340   ImagePointer outputImage = this->GetOutput(0);
 341   outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
 342   outputImage->SetOrigin(m_working_image->GetOrigin());
 343   outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
 344   DD("end2");
 345 }
 346 //--------------------------------------------------------------------
 347
 348
 349 //--------------------------------------------------------------------
 350 template <class TImageType>
 351 void
 352 clitk::ExtractLymphStationsFilter<TImageType>::
 353 GenerateInputRequestedRegion() {
 354   DD("GenerateInputRequestedRegion");
 355   // Call default
 356   Superclass::GenerateInputRequestedRegion();
 357   // Following needed because output region can be greater than input (trachea)
 358   ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
 359   ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
 360   mediastinum->SetRequestedRegion(mediastinum->GetLargestPossibleRegion());
 361   trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
 362 }
 363 //--------------------------------------------------------------------
 364
 365
 366 //--------------------------------------------------------------------
 367 template <class TImageType>
 368 void
 369 clitk::ExtractLymphStationsFilter<TImageType>::
 370 GenerateData() {
 371   DD("GenerateData");
 372   // Final Step -> graft output (if SetNthOutput => redo)
 373   this->GraftOutput(m_output);
 374 }
 375 //--------------------------------------------------------------------
 376
 377
 378 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX