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 #include <itkAddImageFilter.h>
  35 #include <itkDivideByConstantImageFilter.h>
  36
  37 // itk [Bloch et al]
  38 #include "RelativePositionPropImageFilter.h"
  39
  40 //--------------------------------------------------------------------
  41 template <class ImageType>
  42 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  43 AddRelativePositionConstraintToLabelImageFilter():
  44   clitk::FilterBase(),
  45   itk::ImageToImageFilter<ImageType, ImageType>()
  46 {
  47   this->SetNumberOfRequiredInputs(2);
  48   SetFuzzyThreshold(0.6);
  49   SetBackgroundValue(0);
  50   SetObjectBackgroundValue(0);
  51   ClearOrientationType();
  52   IntermediateSpacingFlagOn();
  53   SetIntermediateSpacing(10);
  54   AutoCropFlagOn();
  55   InverseOrientationFlagOff();
  56   RemoveObjectFlagOn();
  57   CombineWithOrFlagOff();
  58 }
  59 //--------------------------------------------------------------------
  60
  61
  62 //--------------------------------------------------------------------
  63 template <class ImageType>
  64 void
  65 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  66 SetInput(const ImageType * image)
  67 {
  68   // Process object is not const-correct so the const casting is required.
  69   this->SetNthInput(0, const_cast<ImageType *>(image));
  70 }
  71 //--------------------------------------------------------------------
  72
  73
  74 //--------------------------------------------------------------------
  75 template <class ImageType>
  76 void
  77 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  78 SetInputObject(const ImageType * image)
  79 {
  80   // Process object is not const-correct so the const casting is required.
  81   this->SetNthInput(1, const_cast<ImageType *>(image));
  82 }
  83 //--------------------------------------------------------------------
  84
  85
  86 //--------------------------------------------------------------------
  87 template <class ImageType>
  88 void
  89 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
  90 ClearOrientationType()
  91 {
  92   m_OrientationTypeString.clear();
  93   m_OrientationType.clear();
  94   m_Angle1.clear();
  95   m_Angle2.clear();
  96 }
  97 //--------------------------------------------------------------------
  98
  99
 100 //--------------------------------------------------------------------
 101 template <class ImageType>
 102 int
 103 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 104 GetNumberOfAngles()
 105 {
 106   return m_OrientationType.size();
 107 }
 108 //--------------------------------------------------------------------
 109
 110
 111 //--------------------------------------------------------------------
 112 template <class ImageType>
 113 void
 114 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 115 AddOrientationTypeString(std::string t)
 116 {
 117   m_OrientationTypeString.push_back(t);
 118   switch (t[0]) {
 119   case 'L' : AddOrientationType(AtLeftTo); break;
 120   case 'R' : AddOrientationType(AtRightTo);break;
 121   case 'A' : AddOrientationType(AntTo);break;
 122   case 'P' : AddOrientationType(PostTo);break;
 123   case 'S' : AddOrientationType(SupTo);break;
 124   case 'I' : AddOrientationType(InfTo);break;
 125   default: clitkExceptionMacro("Error, you must provide L,R or A,P or S,I");
 126   }
 127 }
 128 //--------------------------------------------------------------------
 129
 130
 131 //--------------------------------------------------------------------
 132 template <class ImageType>
 133 void
 134 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 135 GenerateOutputInformation()
 136 {
 137   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 138   ImagePointer outputImage = this->GetOutput(0);
 139   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
 140 }
 141 //--------------------------------------------------------------------
 142
 143
 144 //--------------------------------------------------------------------
 145 template <class ImageType>
 146 void
 147 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 148 GenerateInputRequestedRegion()
 149 {
 150   // Call default
 151   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
 152   // Get input pointers and set requested region to common region
 153   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 154   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 155   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
 156   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
 157 }
 158 //--------------------------------------------------------------------
 159
 160
 161 //--------------------------------------------------------------------
 162 template <class ImageType>
 163 void
 164 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 165 AddAngles(double a, double b)
 166 {
 167   AddOrientationTypeString("Angle");
 168   m_Angle1.push_back(a);
 169   m_Angle2.push_back(b);
 170 }
 171 //--------------------------------------------------------------------
 172
 173
 174 //--------------------------------------------------------------------
 175 template <class ImageType>
 176 void
 177 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 178 AddOrientationType(OrientationTypeEnumeration orientation)
 179 {
 180   m_OrientationType.push_back(orientation);
 181   switch (orientation) {
 182   case AtRightTo:
 183     m_Angle1.push_back(clitk::deg2rad(0));
 184     m_Angle2.push_back(clitk::deg2rad(0));
 185     break;
 186   case AtLeftTo:
 187     m_Angle1.push_back(clitk::deg2rad(180));
 188     m_Angle2.push_back(clitk::deg2rad(0));
 189     break;
 190   case AntTo:
 191     m_Angle1.push_back(clitk::deg2rad(90));
 192     m_Angle2.push_back(clitk::deg2rad(0));
 193     break;
 194   case PostTo:
 195     m_Angle1.push_back(clitk::deg2rad(-90));
 196     m_Angle2.push_back(clitk::deg2rad(0));
 197     break;
 198   case InfTo:
 199     m_Angle1.push_back(clitk::deg2rad(0));
 200     m_Angle2.push_back(clitk::deg2rad(90));
 201     break;
 202   case SupTo:
 203     m_Angle1.push_back(clitk::deg2rad(0));
 204     m_Angle2.push_back(clitk::deg2rad(-90));
 205     break;
 206   case Angle:  break;
 207   }
 208   /*         A1   A2
 209              Left      0    0
 210              Right   180    0
 211              Ant      90    0
 212              Post    -90    0
 213              Inf       0   90
 214              Sup       0  -90
 215   */
 216 }
 217 //--------------------------------------------------------------------
 218
 219
 220 //--------------------------------------------------------------------
 221 template <class ImageType>
 222 void
 223 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 224 GenerateData()
 225 {
 226   if (GetNumberOfAngles() <1) {
 227     clitkExceptionMacro("Add at least one orientation type");
 228   }
 229
 230   // Get input pointer
 231   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 232   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 233   static const unsigned int dim = ImageType::ImageDimension;
 234
 235   // Step 2: object pad to input image -> we want to compute the
 236   // relative position for each point belonging to the input image
 237   // domain, so we have to extend (pad) the object image to fit the
 238   // domain size
 239   working_image = object;
 240   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
 241     StartNewStep("Pad (resize) object to input size");
 242
 243     if (0) { // OLD VERSION (TO REMOVE)
 244       StartNewStep("Pad object to image size");
 245       typename ImageType::Pointer output = ImageType::New();
 246       SizeType size;
 247       for(unsigned int i=0; i<dim; i++) {
 248         size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
 249                          input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
 250       }
 251
 252       // The index of the input is not necessarily zero, so we have to
 253       // take it into account (not done)
 254       RegionType region;
 255       IndexType index = input->GetLargestPossibleRegion().GetIndex();
 256       region.SetSize(size);
 257       for(unsigned int i=0; i<dim; i++) {
 258         if (index[i] != 0) {
 259           std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
 260           exit(0);
 261         }
 262       }
 263       // output->SetLargestPossibleRegion(region);
 264       output->SetRegions(region);
 265       output->SetSpacing(working_image->GetSpacing());
 266       PointType origin = input->GetOrigin();
 267       for(unsigned int i=0; i<dim; i++) {
 268         origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
 269       }
 270       output->SetOrigin(origin);
 271       //    output->SetOrigin(input->GetOrigin());
 272
 273       output->Allocate();
 274       output->FillBuffer(m_BackgroundValue);
 275       typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
 276       // typename PasteFilterType::InputImageIndexType index;
 277       for(unsigned int i=0; i<dim; i++) {
 278         index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
 279           + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
 280       }
 281       padFilter->SetSourceImage(working_image);
 282       padFilter->SetDestinationImage(output);
 283       padFilter->SetDestinationIndex(index);
 284       padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
 285       padFilter->Update();
 286       working_image = padFilter->GetOutput();
 287     }
 288
 289     // Resize object like input
 290     working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
 291     StopCurrentStep<ImageType>(working_image);
 292   }
 293
 294   //--------------------------------------------------------------------
 295   //--------------------------------------------------------------------
 296   // Step 1 : resample
 297   if (m_IntermediateSpacingFlag) {
 298     StartNewStep("Resample object to intermediate spacing");
 299     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
 300     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
 301     resampleFilter->SetInput(working_image);
 302     resampleFilter->SetDefaultPixelValue(0);
 303     resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
 304     resampleFilter->SetGaussianFilteringEnabled(false);
 305     //    resampleFilter->SetVerboseOptions(true);
 306     resampleFilter->Update();
 307     working_image = resampleFilter->GetOutput();
 308     StopCurrentStep<ImageType>(working_image);
 309   }
 310
 311   // Keep object image (with resampline and pad)
 312   object_resampled = working_image;
 313
 314   // Step 3: compute rel pos in object
 315   StartNewStep("Relative Position Map");
 316   typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
 317   typename RelPosFilterType::Pointer relPosFilter;
 318
 319   typename FloatImageType::Pointer m_FuzzyMap;
 320   for(int i=0; i<GetNumberOfAngles(); i++) {
 321     // Compute fuzzy map
 322     relPosFilter = RelPosFilterType::New();
 323     relPosFilter->SetInput(working_image);
 324     relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
 325     relPosFilter->SetAlpha2(m_Angle2[i]);
 326     relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
 327     relPosFilter->SetFast(true);
 328     relPosFilter->SetRadius(1); // seems sufficient in this case
 329     // relPosFilter->SetVerboseProgress(true);
 330     relPosFilter->Update();
 331     relPos = relPosFilter->GetOutput();
 332
 333     if (GetNumberOfAngles() != 1) {
 334       // Creation of the first m_FuzzyMap
 335       if (i==0) {
 336         m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
 337         m_FuzzyMap->FillBuffer(0.0);
 338       }
 339
 340       // Add to current fuzzy map
 341       typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
 342       typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
 343       addFilter->SetInput1(m_FuzzyMap);
 344       addFilter->SetInput2(relPos);
 345       addFilter->Update();
 346       m_FuzzyMap = addFilter->GetOutput();
 347     }
 348     else m_FuzzyMap = relPos;
 349   }
 350
 351   // Divide by the number of relpos
 352   if (GetNumberOfAngles() != 1) {
 353     typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
 354     typename DivideFilter::Pointer divideFilter = DivideFilter::New();
 355     divideFilter->SetInput(m_FuzzyMap);
 356     divideFilter->SetConstant(GetNumberOfAngles());
 357     divideFilter->Update();
 358     m_FuzzyMap = divideFilter->GetOutput();
 359   }
 360
 361   relPos = m_FuzzyMap;
 362   StopCurrentStep<FloatImageType>(relPos);
 363
 364   //--------------------------------------------------------------------
 365   //--------------------------------------------------------------------
 366   StartNewStep("Map Threshold");
 367   // Step 1: threshold
 368   typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
 369   typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
 370   thresholdFilter->SetInput(relPos);
 371   thresholdFilter->SetOutsideValue(m_BackgroundValue);
 372   thresholdFilter->SetInsideValue(m_BackgroundValue+1);
 373   thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
 374   thresholdFilter->Update();
 375   working_image = thresholdFilter->GetOutput();
 376   StopCurrentStep<ImageType>(working_image);
 377
 378   //--------------------------------------------------------------------
 379   //--------------------------------------------------------------------
 380   StartNewStep("Post Processing: erosion with initial mask");
 381   // Step 2 : erosion with initial mask to exclude pixels that were
 382   // inside the resampled version and outside the original mask
 383   typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
 384   StructuringElementType kernel;
 385   kernel.SetRadius(1);
 386   kernel.CreateStructuringElement();
 387   typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
 388   typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
 389   erodeFilter->SetInput(working_image);
 390   erodeFilter->SetKernel(kernel);
 391   erodeFilter->SetBackgroundValue(m_BackgroundValue);
 392   erodeFilter->SetErodeValue(m_BackgroundValue+1);
 393   erodeFilter->Update();
 394   working_image = erodeFilter->GetOutput();
 395   StopCurrentStep<ImageType>(working_image);
 396
 397   //--------------------------------------------------------------------
 398   //--------------------------------------------------------------------
 399   // Step 5: resample to initial spacing
 400   if (m_IntermediateSpacingFlag) {
 401     StartNewStep("Resample to come back to the same sampling than input");
 402     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
 403     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
 404     resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
 405     resampleFilter->SetInput(working_image);
 406     resampleFilter->SetOutputSpacing(input->GetSpacing());
 407     resampleFilter->SetGaussianFilteringEnabled(false);
 408     // resampleFilter->SetVerboseOptions(true);
 409     resampleFilter->Update();
 410     working_image = resampleFilter->GetOutput();
 411     StopCurrentStep<ImageType>(working_image);
 412   }
 413
 414   //--------------------------------------------------------------------
 415   //--------------------------------------------------------------------
 416   // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
 417   //if (!HaveSameSizeAndSpacing(working_image, input)) {
 418   if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
 419     StartNewStep("Pad to get the same size than input");
 420     typename ImageType::Pointer temp = ImageType::New();
 421     temp->CopyInformation(input);
 422     temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
 423     temp->Allocate();
 424     temp->FillBuffer(m_BackgroundValue);
 425     typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
 426     padFilter2->SetSourceImage(working_image);
 427     padFilter2->SetDestinationImage(temp);
 428     // DD(input->GetLargestPossibleRegion().GetIndex());
 429     padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
 430     padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
 431     padFilter2->Update();
 432     working_image = padFilter2->GetOutput();
 433     StopCurrentStep<ImageType>(working_image);
 434   }
 435   else {
 436     //DD("[debug] Rel Pos : no padding after");
 437   }
 438
 439   //--------------------------------------------------------------------
 440   //--------------------------------------------------------------------
 441   // Step 6: combine input+thresholded relpos
 442   StartNewStep("Combine with initial input (boolean And)");
 443   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 444   typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
 445   combineFilter->SetBackgroundValue(m_BackgroundValue);
 446   combineFilter->SetBackgroundValue1(m_BackgroundValue);
 447   combineFilter->SetBackgroundValue2(m_BackgroundValue);
 448   combineFilter->SetForegroundValue(m_BackgroundValue+1);
 449   combineFilter->SetInput1(input);
 450   combineFilter->SetInput2(working_image);
 451   if (GetInverseOrientationFlag())
 452     combineFilter->SetOperationType(BoolFilterType::AndNot);
 453   else {
 454     if (GetCombineWithOrFlag())
 455       combineFilter->SetOperationType(BoolFilterType::Or);
 456     else
 457       combineFilter->SetOperationType(BoolFilterType::And);
 458   }
 459   combineFilter->InPlaceOff(); // Do not modify initial input (!)
 460   combineFilter->Update();
 461   working_image = combineFilter->GetOutput();
 462
 463   // Remove (if needed the object from the support)
 464   if (GetRemoveObjectFlag()) {
 465     combineFilter = BoolFilterType::New();
 466     combineFilter->SetInput1(working_image);
 467     combineFilter->SetInput2(object);
 468     combineFilter->SetOperationType(BoolFilterType::AndNot);
 469     combineFilter->InPlaceOn();
 470     combineFilter->Update();
 471     working_image = combineFilter->GetOutput();
 472   }
 473   // In the other case, we must *add* the initial object to keep it
 474   // but not more than the initial support
 475   else {
 476     combineFilter = BoolFilterType::New();
 477     combineFilter->SetInput1(working_image);
 478     combineFilter->SetInput2(object);
 479     combineFilter->SetOperationType(BoolFilterType::Or);
 480     combineFilter->InPlaceOn();
 481     combineFilter->Update();
 482     working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
 483     combineFilter = BoolFilterType::New();
 484     combineFilter->SetInput1(working_image);
 485     combineFilter->SetInput2(input);
 486     combineFilter->SetOperationType(BoolFilterType::And);
 487     combineFilter->InPlaceOn();
 488     combineFilter->Update();
 489     working_image = combineFilter->GetOutput();
 490   }
 491
 492   StopCurrentStep<ImageType>(working_image);
 493
 494   //--------------------------------------------------------------------
 495   //--------------------------------------------------------------------
 496   // Step 7: autocrop
 497   if (GetAutoCropFlag()) {
 498     StartNewStep("Final AutoCrop");
 499     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
 500     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 501     cropFilter->SetInput(working_image);
 502     cropFilter->ReleaseDataFlagOff();
 503     cropFilter->Update();
 504     working_image = cropFilter->GetOutput();
 505     StopCurrentStep<ImageType>(working_image);
 506   }
 507
 508   //--------------------------------------------------------------------
 509   //--------------------------------------------------------------------
 510
 511   // Final Step -> set output
 512   this->SetNthOutput(0, working_image);
 513   //  this->GraftOutput(working_image);
 514 }
 515 //--------------------------------------------------------------------
 516