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