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