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