segmentation/clitkExtractLymphStationsFilter.txx.save

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