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 ImageType>
  40 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  41 AddRelativePositionConstraintToLabelImageFilter():
  42   clitk::FilterBase(),
  43   itk::ImageToImageFilter<ImageType, ImageType>()
  44 {
  45   this->SetNumberOfRequiredInputs(2);
  46   SetFuzzyThreshold(0.6);
  47   SetBackgroundValue(0);
  48   SetObjectBackgroundValue(0);
  49   SetOrientationType(LeftTo);
  50   ResampleBeforeRelativePositionFilterOn();
  51   SetIntermediateSpacing(10);
  52   AutoCropFlagOn();
  53   NotFlagOff();
  54 }
  55 //--------------------------------------------------------------------
  56
  57
  58 //--------------------------------------------------------------------
  59 template <class ImageType>
  60 void
  61 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  62 SetInput(const ImageType * image)
  63 {
  64   // Process object is not const-correct so the const casting is required.
  65   this->SetNthInput(0, const_cast<ImageType *>(image));
  66 }
  67 //--------------------------------------------------------------------
  68
  69
  70 //--------------------------------------------------------------------
  71 template <class ImageType>
  72 void
  73 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  74 SetInputObject(const ImageType * image)
  75 {
  76   // Process object is not const-correct so the const casting is required.
  77   this->SetNthInput(1, const_cast<ImageType *>(image));
  78 }
  79 //--------------------------------------------------------------------
  80
  81
  82 //--------------------------------------------------------------------
  83 template <class ImageType>
  84 void
  85 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  86 SetOrientationTypeString(std::string t)
  87 {
  88   SetOrientationType(Angle);
  89   if (t[0] == 'L') SetOrientationType(LeftTo);
  90   if (t[0] == 'R') SetOrientationType(RightTo);
  91   if (t[0] == 'A') SetOrientationType(AntTo);
  92   if (t[0] == 'P') SetOrientationType(PostTo);
  93   if (t[0] == 'S') SetOrientationType(SupTo);
  94   if (t[0] == 'I') SetOrientationType(InfTo);
  95 }
  96 //--------------------------------------------------------------------
  97
  98
  99 //--------------------------------------------------------------------
 100 template <class ImageType>
 101 void
 102 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 103 GenerateOutputInformation()
 104 {
 105   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 106   ImagePointer outputImage = this->GetOutput(0);
 107   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
 108 }
 109 //--------------------------------------------------------------------
 110
 111
 112 //--------------------------------------------------------------------
 113 template <class ImageType>
 114 void
 115 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 116 GenerateInputRequestedRegion()
 117 {
 118   // Call default
 119   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
 120   // Get input pointers and set requested region to common region
 121   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 122   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 123   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
 124   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
 125 }
 126 //--------------------------------------------------------------------
 127
 128
 129 //--------------------------------------------------------------------
 130 template <class ImageType>
 131 void
 132 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 133 SetAngle1(double a)
 134 {
 135   SetOrientationType(Angle);
 136   m_Angle1 = a;
 137 }
 138 //--------------------------------------------------------------------
 139
 140
 141 //--------------------------------------------------------------------
 142 template <class ImageType>
 143 void
 144 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 145 SetAngle2(double a)
 146 {
 147   SetOrientationType(Angle);
 148   m_Angle2 = a;
 149 }
 150 //--------------------------------------------------------------------
 151
 152
 153 //--------------------------------------------------------------------
 154 template <class ImageType>
 155 void
 156 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 157 SetOrientationType(OrientationTypeEnumeration orientation)
 158 {
 159   m_OrientationType = orientation;
 160   switch (m_OrientationType) {
 161   case LeftTo:   m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(0);   break;
 162   case RightTo:  m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0);   break;
 163   case AntTo:    m_Angle1 = clitk::deg2rad(90);  m_Angle2 = clitk::deg2rad(0);   break;
 164   case PostTo:   m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0);   break;
 165   case InfTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(90);  break;
 166   case SupTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(-90); break;
 167   case Angle:  break;
 168   }
 169   /*         A1   A2
 170              Left      0    0
 171              Right   180    0
 172              Ant      90    0
 173              Post    -90    0
 174              Inf       0   90
 175              Sup       0  -90
 176   */
 177 }
 178 //--------------------------------------------------------------------
 179
 180
 181 //--------------------------------------------------------------------
 182 template <class ImageType>
 183 void
 184 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 185 GenerateData()
 186 {
 187   // Print Option
 188   /*
 189     DD(GetFuzzyThreshold());
 190     DD((int)GetBackgroundValue());
 191     DD((int)GetObjectBackgroundValue());
 192     DD(GetOrientationType());
 193     DD(GetResampleBeforeRelativePositionFilter());
 194     DD(GetIntermediateSpacing());
 195     DD(GetAutoCropFlag());
 196     DD(GetNotFlag());
 197   */
 198
 199   // Get input pointer
 200   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 201   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 202
 203   //--------------------------------------------------------------------
 204   //--------------------------------------------------------------------
 205   static const unsigned int dim = ImageType::ImageDimension;
 206   StartNewStep("Initial resample");
 207   // Step 1 : resample
 208   if (m_ResampleBeforeRelativePositionFilter) {
 209     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
 210     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
 211     resampleFilter->SetInput(object);
 212     resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
 213     resampleFilter->SetGaussianFilteringEnabled(false);
 214     //    resampleFilter->SetVerboseOptions(true);
 215     resampleFilter->Update();
 216     working_image = resampleFilter->GetOutput();
 217   }
 218   else {
 219     working_image = object;
 220   }
 221   StopCurrentStep<ImageType>(working_image);
 222
 223   // Step 2: object pad to input image -> we want to compute the
 224   // relative position for each point belonging to the input image
 225   // domain, so we have to extend (pad) the object image to fit the
 226   // domain size
 227   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
 228     StartNewStep("Pad object to image size");
 229     typename ImageType::Pointer output = ImageType::New();
 230     SizeType size;
 231     for(unsigned int i=0; i<dim; i++) {
 232       size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
 233     }
 234
 235     // The index of the input is not necessarily zero, so we have to
 236     // take it into account (not done)
 237     RegionType region;
 238     IndexType index = input->GetLargestPossibleRegion().GetIndex();
 239     region.SetSize(size);
 240     for(unsigned int i=0; i<dim; i++) {
 241       if (index[i] != 0) {
 242         std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
 243         exit(0);
 244       }
 245     }
 246     // output->SetLargestPossibleRegion(region);
 247     output->SetRegions(region);
 248     output->SetSpacing(working_image->GetSpacing());
 249     PointType origin = input->GetOrigin();
 250     for(unsigned int i=0; i<dim; i++) {
 251       origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
 252     }
 253     output->SetOrigin(origin);
 254     //    output->SetOrigin(input->GetOrigin());
 255
 256     output->Allocate();
 257     output->FillBuffer(m_BackgroundValue);
 258     typename PadFilterType::Pointer padFilter = PadFilterType::New();
 259     // typename PadFilterType::InputImageIndexType index;
 260     for(unsigned int i=0; i<dim; i++) {
 261       index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
 262     + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
 263     }
 264     padFilter->SetSourceImage(working_image);
 265     padFilter->SetDestinationImage(output);
 266     padFilter->SetDestinationIndex(index);
 267     padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
 268     padFilter->Update();
 269     working_image = padFilter->GetOutput();
 270     StopCurrentStep<ImageType>(working_image);
 271   }
 272   else {
 273     // DD("[debug] RelPos : same size and spacing : no padding");
 274   }
 275   // Keep object image (with resampline and pad)
 276   object_resampled = working_image;
 277  //  StopCurrentStep<ImageType>(working_image);
 278
 279   // Step 3: compute rel pos in object
 280   StartNewStep("Relative Position Map");
 281   typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
 282   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 283   relPosFilter->SetInput(working_image);
 284   relPosFilter->SetAlpha1(m_Angle1); // xy plane
 285   relPosFilter->SetAlpha2(m_Angle2);
 286   relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
 287   relPosFilter->SetFast(true);
 288   relPosFilter->SetRadius(1); // seems sufficient in this cas
 289   // relPosFilter->SetVerboseProgress(true);
 290   relPosFilter->Update();
 291   relPos = relPosFilter->GetOutput();
 292   StopCurrentStep<FloatImageType>(relPos);
 293
 294   //--------------------------------------------------------------------
 295   //--------------------------------------------------------------------
 296   StartNewStep("Map Threshold");
 297   // Step 1: threshold
 298   typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
 299   typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
 300   thresholdFilter->SetInput(relPos);
 301   thresholdFilter->SetOutsideValue(m_BackgroundValue);
 302   thresholdFilter->SetInsideValue(m_BackgroundValue+1);
 303   thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
 304   thresholdFilter->Update();
 305   working_image = thresholdFilter->GetOutput();
 306   StopCurrentStep<ImageType>(working_image);
 307
 308   //--------------------------------------------------------------------
 309   //--------------------------------------------------------------------
 310   StartNewStep("Post Processing: erosion with initial mask");
 311   // Step 2 : erosion with initial mask to exclude pixels that were
 312   // inside the resampled version and outside the original mask
 313   typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
 314   StructuringElementType kernel;
 315   kernel.SetRadius(1);
 316   kernel.CreateStructuringElement();
 317   typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
 318   typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
 319   erodeFilter->SetInput(working_image);
 320   erodeFilter->SetKernel(kernel);
 321   erodeFilter->SetBackgroundValue(m_BackgroundValue);
 322   erodeFilter->SetErodeValue(m_BackgroundValue+1);
 323   erodeFilter->Update();
 324   working_image = erodeFilter->GetOutput();
 325   StopCurrentStep<ImageType>(working_image);
 326
 327   //--------------------------------------------------------------------
 328   //--------------------------------------------------------------------
 329   // Step 5: resample to initial spacing
 330   if (m_ResampleBeforeRelativePositionFilter) {
 331     StartNewStep("Resample to get the same sampling than input");
 332     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
 333     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
 334     resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
 335     resampleFilter->SetInput(working_image);
 336     resampleFilter->SetOutputSpacing(input->GetSpacing());
 337     resampleFilter->SetGaussianFilteringEnabled(false);
 338     // resampleFilter->SetVerboseOptions(true);
 339     resampleFilter->Update();
 340     working_image = resampleFilter->GetOutput();
 341     StopCurrentStep<ImageType>(working_image);
 342   }
 343
 344   //--------------------------------------------------------------------
 345   //--------------------------------------------------------------------
 346   // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
 347   // DD(working_image->GetLargestPossibleRegion());
 348   // DD(input->GetLargestPossibleRegion());
 349   //if (!HaveSameSizeAndSpacing(working_image, input)) {
 350   if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
 351     StartNewStep("Pad to get the same size than input");
 352     typename ImageType::Pointer temp = ImageType::New();
 353     temp->CopyInformation(input);
 354     temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
 355     temp->Allocate();
 356     temp->FillBuffer(m_BackgroundValue);
 357     typename PadFilterType::Pointer padFilter2 = PadFilterType::New();
 358     padFilter2->SetSourceImage(working_image);
 359     padFilter2->SetDestinationImage(temp);
 360     // DD(input->GetLargestPossibleRegion().GetIndex());
 361     padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
 362     padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
 363     padFilter2->Update();
 364     working_image = padFilter2->GetOutput();
 365     StopCurrentStep<ImageType>(working_image);
 366   }
 367   else {
 368     //DD("[debug] Rel Pos : no padding after");
 369   }
 370
 371   //--------------------------------------------------------------------
 372   //--------------------------------------------------------------------
 373   // Step 6: combine input+thresholded relpos
 374   StartNewStep("Combine with initial input (boolean And)");
 375   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 376   typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
 377   combineFilter->SetBackgroundValue(m_BackgroundValue);
 378   combineFilter->SetBackgroundValue1(m_BackgroundValue);
 379   combineFilter->SetBackgroundValue2(m_BackgroundValue);
 380   combineFilter->SetForegroundValue(m_BackgroundValue+1);
 381   combineFilter->SetInput1(input);
 382   combineFilter->SetInput2(working_image);
 383   if (GetNotFlag())
 384     combineFilter->SetOperationType(BoolFilterType::AndNot);
 385   else
 386     combineFilter->SetOperationType(BoolFilterType::And);
 387   combineFilter->InPlaceOn();
 388   combineFilter->Update();
 389   working_image = combineFilter->GetOutput();
 390
 391   combineFilter = BoolFilterType::New();
 392   combineFilter->SetInput1(working_image);
 393   combineFilter->SetInput2(object);
 394   combineFilter->SetOperationType(BoolFilterType::AndNot);
 395   combineFilter->InPlaceOn();
 396   combineFilter->Update();
 397
 398   working_image = combineFilter->GetOutput();
 399   StopCurrentStep<ImageType>(working_image);
 400
 401   //--------------------------------------------------------------------
 402   //--------------------------------------------------------------------
 403   // Step 7: autocrop
 404   if (GetAutoCropFlag()) {
 405     StartNewStep("Final AutoCrop");
 406     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
 407     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 408     cropFilter->SetInput(working_image);
 409     cropFilter->ReleaseDataFlagOff();
 410     cropFilter->Update();
 411     working_image = cropFilter->GetOutput();
 412     StopCurrentStep<ImageType>(working_image);
 413   }
 414
 415   //--------------------------------------------------------------------
 416   //--------------------------------------------------------------------
 417
 418   // Final Step -> set output
 419   this->SetNthOutput(0, working_image);
 420   //  this->GraftOutput(working_image);
 421   return;
 422 }
 423 //--------------------------------------------------------------------
 424