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