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