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