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