itk/clitkAddRelativePositionConstraintToLabelImageFilter.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 // clitk
  20 #include "clitkCommon.h"
  21 #include "clitkBooleanOperatorLabelImageFilter.h"
  22 #include "clitkAutoCropFilter.h"
  23 #include "clitkResampleImageWithOptionsFilter.h"
  24 #include "clitkBooleanOperatorLabelImageFilter.h"
  25
  26 // itk
  27 #include <deque>
  28 #include "itkStatisticsLabelMapFilter.h"
  29 #include "itkLabelImageToStatisticsLabelMapFilter.h"
  30 #include "itkRegionOfInterestImageFilter.h"
  31 #include "itkBinaryThresholdImageFilter.h"
  32 #include "itkBinaryErodeImageFilter.h"
  33 #include "itkBinaryBallStructuringElement.h"
  34
  35 // itk [Bloch et al]
  36 #include "RelativePositionPropImageFilter.h"
  37
  38 //--------------------------------------------------------------------
  39 template <class TImageType>
  40 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
  41 AddRelativePositionConstraintToLabelImageFilter():
  42   clitk::FilterBase(),
  43   itk::ImageToImageFilter<TImageType, TImageType>()
  44 {
  45   this->SetNumberOfRequiredInputs(2);
  46   SetFuzzyThreshold(0.6);
  47   SetBackgroundValue(0);
  48   SetObjectBackgroundValue(0);
  49   SetOrientationType(LeftTo);
  50   ResampleBeforeRelativePositionFilterOn();
  51   SetIntermediateSpacing(10);
  52
  53   // Step 1 : resample (option=sampling)
  54   // Step 2 : Pad (no)
  55   // Step 3 : relative position  (option = angle)
  56   // Step 4 : Threshold
  57   // Step 5 : Erode for boundary
  58   // Step 6 : resample to initial spacing
  59   // Step 7 : pad if not the same size : it can occur when downsample and upsample
  60   // Step 6: combine input+thresholded relpos
  61   // Step 7: autocrop
  62
  63   // Step 1 : resample
  64   ResampleBeforeRelativePositionFilterOn();
  65   SetIntermediateSpacing(10);
  66   SetOrientationType(LeftTo);
  67   // SegmentationStep * step = new SegmentationStep;
  68   //   step->SetName("Initial resampling and relative position map");
  69   //   SetStep(step, &Self::ResampleAndRelativePositionMap);
  70
  71   // Step 3 : threshold + postprocess
  72   // SetFuzzyThreshold(0.4);
  73   //  step = new SegmentationStep;
  74   //   step->SetName("Map threshold and post-processing");
  75   //   SetStep(step, &Self::MapThresholdAndPostProcessing);
  76 }
  77 //--------------------------------------------------------------------
  78
  79
  80 //--------------------------------------------------------------------
  81 template <class TImageType>
  82 void
  83 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
  84 SetInput(const ImageType * image) {
  85   // Process object is not const-correct so the const casting is required.
  86   this->SetNthInput(0, const_cast<ImageType *>(image));
  87 }
  88 //--------------------------------------------------------------------
  89
  90
  91 //--------------------------------------------------------------------
  92 template <class TImageType>
  93 void
  94 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
  95 SetInputObject(const ImageType * image) {
  96   // Process object is not const-correct so the const casting is required.
  97   this->SetNthInput(1, const_cast<ImageType *>(image));
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 template <class TImageType>
 104 void
 105 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
 106 GenerateOutputInformation() {
 107   ImagePointer input = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
 108   ImagePointer outputImage = this->GetOutput(0);
 109   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
 110 }
 111 //--------------------------------------------------------------------
 112
 113
 114 //--------------------------------------------------------------------
 115 template <class TImageType>
 116 void
 117 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
 118 GenerateInputRequestedRegion() {
 119   // Call default
 120   itk::ImageToImageFilter<TImageType, TImageType>::GenerateInputRequestedRegion();
 121   // Get input pointers and set requested region to common region
 122   ImagePointer input1 = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
 123   ImagePointer input2 = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
 124   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
 125   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
 126 }
 127 //--------------------------------------------------------------------
 128
 129
 130 //--------------------------------------------------------------------
 131 template <class TImageType>
 132 void
 133 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
 134 SetAngle1(double a) {
 135   SetOrientationType(Angle);
 136   m_Angle1 = a;
 137 }
 138 //--------------------------------------------------------------------
 139
 140
 141 //--------------------------------------------------------------------
 142 template <class TImageType>
 143 void
 144 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
 145 SetAngle2(double a) {
 146   SetOrientationType(Angle);
 147   m_Angle2 = a;
 148 }
 149 //--------------------------------------------------------------------
 150
 151
 152 //--------------------------------------------------------------------
 153 template <class TImageType>
 154 void
 155 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
 156 SetOrientationType(OrientationTypeEnumeration orientation) {
 157   m_OrientationType = orientation;
 158   switch (m_OrientationType) {
 159   case LeftTo:   m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(0);   break;
 160   case RightTo:  m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0);   break;
 161   case AntTo:    m_Angle1 = clitk::deg2rad(90);  m_Angle2 = clitk::deg2rad(0);   break;
 162   case PostTo:   m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0);   break;
 163   case InfTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(90);  break;
 164   case SupTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(-90); break;
 165   case Angle:  break;
 166   }
 167   /*         A1   A2
 168              Left      0    0
 169              Right   180    0
 170              Ant      90    0
 171              Post    -90    0
 172              Inf       0   90
 173              Sup       0  -90
 174   */
 175 }
 176 //--------------------------------------------------------------------
 177
 178
 179 //--------------------------------------------------------------------
 180 template <class TImageType>
 181 void
 182 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
 183 GenerateData() {
 184   // Get input pointer
 185   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 186   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 187
 188   //--------------------------------------------------------------------
 189   //--------------------------------------------------------------------
 190   StartNewStep("Resample And Relative Position Map");
 191
 192   static const unsigned int dim = ImageType::ImageDimension;
 193   // Step 1 : resample
 194   if (m_ResampleBeforeRelativePositionFilter) {
 195     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
 196     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
 197     resampleFilter->SetInput(object);
 198     resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
 199     resampleFilter->SetGaussianFilteringEnabled(false);
 200     // resampleFilter->SetVerboseOptions(true);
 201     resampleFilter->Update();
 202     working_image = resampleFilter->GetOutput();
 203   }
 204   else {
 205     working_image = object;
 206   }
 207   // writeImage<ImageType>(working_image, "res.mhd");
 208
 209   // Step 2: object pad to input image -> we want to compute the
 210   // relative position for each point belonging to the input image
 211   // domain, so we have to extend (pad) the object image to fit the
 212   // domain size
 213   typename ImageType::Pointer output = ImageType::New();
 214   SizeType size;
 215   for(unsigned int i=0; i<dim; i++) {
 216     size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
 217   }
 218   RegionType region;
 219   region.SetSize(size);
 220   // output->SetLargestPossibleRegion(region);
 221   output->SetRegions(region);
 222   output->SetSpacing(working_image->GetSpacing());
 223   output->SetOrigin(input->GetOrigin());
 224   output->Allocate();
 225   output->FillBuffer(m_BackgroundValue);
 226   typename PadFilterType::Pointer padFilter = PadFilterType::New();
 227   typename PadFilterType::InputImageIndexType index;
 228   for(unsigned int i=0; i<dim; i++) {
 229     index[i] = lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/(double)m_IntermediateSpacing);
 230   }
 231   padFilter->SetSourceImage(working_image);
 232   padFilter->SetDestinationImage(output);
 233   padFilter->SetDestinationIndex(index);
 234   padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
 235   padFilter->Update();
 236   working_image = padFilter->GetOutput();
 237
 238   // Keep object image (with resampline and pad)
 239   object_resampled = working_image;
 240  //  writeImage<ImageType>(working_image, "pad.mhd");
 241
 242   // Step 3: compute rel pos in object
 243   typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
 244   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 245   relPosFilter->SetInput(working_image);
 246   relPosFilter->SetAlpha1(m_Angle1); // xy plane
 247   relPosFilter->SetAlpha2(m_Angle2);
 248   relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
 249   relPosFilter->SetFast(true);
 250   relPosFilter->SetRadius(1); // seems sufficient in this cas
 251   // relPosFilter->SetVerboseProgress(true);
 252   relPosFilter->Update();
 253   relPos = relPosFilter->GetOutput();
 254   this->template StopCurrentStep<FloatImageType>(relPos);
 255
 256   //--------------------------------------------------------------------
 257   //--------------------------------------------------------------------
 258   StartNewStep("Map Threshold And Post Processing");
 259
 260   // Step 1: threshold
 261   typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
 262   typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
 263   thresholdFilter->SetInput(relPos);
 264   thresholdFilter->SetOutsideValue(m_BackgroundValue);
 265   thresholdFilter->SetInsideValue(m_BackgroundValue+1);
 266   thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
 267   thresholdFilter->Update();
 268   working_image = thresholdFilter->GetOutput();
 269
 270   // Step 2 : erosion with initial mask to exclude pixels that were
 271   // inside the resampled version and outside the original mask
 272   typedef itk::BinaryBallStructuringElement<unsigned int, 3> StructuringElementType;
 273   StructuringElementType kernel;
 274   kernel.SetRadius(1);
 275   kernel.CreateStructuringElement();
 276   typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
 277   typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
 278   erodeFilter->SetInput(working_image);
 279   erodeFilter->SetKernel(kernel);
 280   erodeFilter->SetBackgroundValue(m_BackgroundValue);
 281   erodeFilter->SetErodeValue(m_BackgroundValue+1);
 282   erodeFilter->Update();
 283   working_image = erodeFilter->GetOutput();
 284
 285   // Step 5: resample to initial spacing
 286   if (m_ResampleBeforeRelativePositionFilter) {
 287     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
 288     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
 289     resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
 290     resampleFilter->SetInput(working_image);
 291     resampleFilter->SetOutputSpacing(input->GetSpacing());
 292     resampleFilter->SetGaussianFilteringEnabled(false);
 293     // resampleFilter->SetVerboseOptions(true);
 294     resampleFilter->Update();
 295     working_image = resampleFilter->GetOutput();
 296   }
 297
 298   // writeImage<ImageType>(working_image, "resinitial.mhd");
 299
 300   // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
 301   if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
 302     typename ImageType::Pointer temp = ImageType::New();
 303     temp->CopyInformation(input);
 304     temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
 305     temp->Allocate();
 306     temp->FillBuffer(m_BackgroundValue);
 307     typename PadFilterType::Pointer padFilter2 = PadFilterType::New(); // if yes : redo relpos
 308     padFilter2->SetSourceImage(working_image);
 309     padFilter2->SetDestinationImage(temp);
 310     padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
 311     padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
 312     padFilter2->Update();
 313     working_image = padFilter2->GetOutput();
 314   }
 315   // writeImage<ImageType>(working_image, "pad2.mhd");
 316
 317   // Step 6: combine input+thresholded relpos
 318   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 319   typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
 320   combineFilter->SetBackgroundValue(m_BackgroundValue);
 321   combineFilter->SetBackgroundValue1(m_BackgroundValue);
 322   combineFilter->SetBackgroundValue2(m_BackgroundValue);
 323   combineFilter->SetForegroundValue(m_BackgroundValue+1);
 324   combineFilter->SetInput1(input);
 325   combineFilter->SetInput2(working_image);
 326   combineFilter->SetOperationType(BoolFilterType::And);
 327   combineFilter->InPlaceOn();
 328   combineFilter->Update();
 329   working_image = combineFilter->GetOutput();
 330
 331   combineFilter = BoolFilterType::New();
 332   combineFilter->SetInput1(working_image);
 333   combineFilter->SetInput2(object);
 334   combineFilter->SetOperationType(BoolFilterType::AndNot);
 335   combineFilter->InPlaceOn();
 336   combineFilter->Update();
 337
 338   working_image = combineFilter->GetOutput();
 339   // writeImage<ImageType>(working_image, "combine.mhd");
 340
 341   // Step 7: autocrop
 342   typedef clitk::AutoCropFilter<ImageType> CropFilterType;
 343   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 344   cropFilter->SetInput(working_image);
 345   cropFilter->ReleaseDataFlagOff();
 346   cropFilter->Update();
 347   working_image = cropFilter->GetOutput();
 348   this->template StopCurrentStep<ImageType>(working_image);
 349
 350   //--------------------------------------------------------------------
 351   //--------------------------------------------------------------------
 352
 353   // Final Step -> set output
 354   this->SetNthOutput(0, working_image);
 355   return;
 356 }
 357 //--------------------------------------------------------------------
 358