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