X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkAddRelativePositionConstraintToLabelImageFilter.txx;h=e25f98b829fa3c65ed497f5fb5ffd73de863a066;hb=1d365ed1e37fdf249a317a721140288f4dc259cc;hp=1e06168877a91929bf2ed06444e4e4da9d06735e;hpb=e9c2298f9582aab29abe164c72caa37f089f3cd6;p=clitk.git

diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
index 1e06168..e25f98b 100644
--- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
+++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
@@ -3,7 +3,7 @@
 
   Authors belong to: 
   - University of LYON              http://www.universite-lyon.fr/
-  - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
+  - Léon Bérard cancer center       http://www.centreleonberard.fr
   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
 
   This software is distributed WITHOUT ANY WARRANTY; without even
@@ -14,7 +14,7 @@
 
   - BSD        See included LICENSE.txt file
   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-  ======================================================================-====*/
+  ===========================================================================**/
 
 // clitk
 #include "clitkCommon.h"
@@ -22,15 +22,22 @@
 #include "clitkAutoCropFilter.h"
 #include "clitkResampleImageWithOptionsFilter.h"
 #include "clitkBooleanOperatorLabelImageFilter.h"
+#include "clitkCropLikeImageFilter.h"
 
 // itk
 #include <deque>
-#include "itkStatisticsLabelMapFilter.h"
-#include "itkLabelImageToStatisticsLabelMapFilter.h"
-#include "itkRegionOfInterestImageFilter.h"
-#include "itkBinaryThresholdImageFilter.h"
-#include "itkBinaryErodeImageFilter.h"
-#include "itkBinaryBallStructuringElement.h"
+#include <itkStatisticsLabelMapFilter.h>
+#include <itkLabelImageToStatisticsLabelMapFilter.h>
+#include <itkRegionOfInterestImageFilter.h>
+#include <itkBinaryThresholdImageFilter.h>
+#include <itkBinaryErodeImageFilter.h>
+#include <itkBinaryBallStructuringElement.h>
+#include <itkAddImageFilter.h>
+#if ITK_VERSION_MAJOR >= 4
+  #include <itkDivideImageFilter.h>
+#else
+  #include <itkDivideByConstantImageFilter.h>
+#endif
 
 // itk [Bloch et al] 
 #include "RelativePositionPropImageFilter.h"
@@ -46,10 +53,16 @@ AddRelativePositionConstraintToLabelImageFilter():
   SetFuzzyThreshold(0.6);
   SetBackgroundValue(0);
   SetObjectBackgroundValue(0);
-  SetOrientationType(LeftTo);
-  ResampleBeforeRelativePositionFilterOn();
+  ClearOrientationType();
+  IntermediateSpacingFlagOn();
   SetIntermediateSpacing(10);
-  AutoCropOn();
+  AutoCropFlagOn();
+  InverseOrientationFlagOff();
+  RemoveObjectFlagOn();
+  CombineWithOrFlagOff();
+  VerboseStepFlagOff();
+  WriteStepFlagOff();
+  FuzzyMapOnlyFlagOff();
 }
 //--------------------------------------------------------------------
 
@@ -78,6 +91,60 @@ SetInputObject(const ImageType * image)
 //--------------------------------------------------------------------
   
 
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+ClearOrientationType() 
+{
+  m_OrientationTypeString.clear();
+  m_OrientationType.clear();
+  m_Angle1.clear();
+  m_Angle2.clear();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+int
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+GetNumberOfAngles()
+{
+  return m_OrientationType.size();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+AddOrientationTypeString(std::string t) 
+{
+  m_OrientationTypeString.push_back(t);
+
+  if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
+  if (t == "RightTo") { AddOrientationType(RightTo); return; }
+  if (t == "AntTo") { AddOrientationType(AntTo); return; }
+  if (t == "PostTo") { AddOrientationType(PostTo); return; }
+  if (t == "SupTo") { AddOrientationType(SupTo); return; }
+  if (t == "InfTo") { AddOrientationType(InfTo); return; }
+
+  if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
+  if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
+  if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
+  if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
+  if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
+  if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
+
+  if (t == "Angle") return;
+
+  clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
+}
+//--------------------------------------------------------------------
+  
+
 //--------------------------------------------------------------------
 template <class ImageType>
 void 
@@ -112,10 +179,12 @@ GenerateInputRequestedRegion()
 template <class ImageType>
 void 
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
-SetAngle1(double a) 
+AddAnglesInRad(double a, double b) 
 {
-  SetOrientationType(Angle);
-  m_Angle1 = a;
+  m_OrientationTypeString.push_back("Angle");
+  m_OrientationType.push_back(Angle);
+  m_Angle1.push_back(a);
+  m_Angle2.push_back(b);
 }
 //--------------------------------------------------------------------
 
@@ -124,10 +193,9 @@ SetAngle1(double a)
 template <class ImageType>
 void 
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
-SetAngle2(double a) 
+AddAnglesInDeg(double a, double b) 
 {
-  SetOrientationType(Angle);
-  m_Angle2 = a;
+  AddAnglesInRad(clitk::deg2rad(a), clitk::deg2rad(b));
 }
 //--------------------------------------------------------------------
 
@@ -136,17 +204,35 @@ SetAngle2(double a)
 template <class ImageType>
 void 
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
-SetOrientationType(OrientationTypeEnumeration orientation) 
+AddOrientationType(OrientationTypeEnumeration orientation) 
 {
-  m_OrientationType = orientation;
-  switch (m_OrientationType) {
-  case LeftTo:   m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(0);   break;
-  case RightTo:  m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0);   break;
-  case AntTo:    m_Angle1 = clitk::deg2rad(90);  m_Angle2 = clitk::deg2rad(0);   break;
-  case PostTo:   m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0);   break;
-  case InfTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(90);  break;
-  case SupTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(-90); break;
-  case Angle:  break;      
+  m_OrientationType.push_back(orientation);
+  switch (orientation) {
+  case RightTo:   
+    m_Angle1.push_back(clitk::deg2rad(0));   
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case LeftTo:  
+    m_Angle1.push_back(clitk::deg2rad(180)); 
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case AntTo:
+    m_Angle1.push_back(clitk::deg2rad(90));
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case PostTo:
+    m_Angle1.push_back(clitk::deg2rad(-90)); 
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case InfTo:    
+    m_Angle1.push_back(clitk::deg2rad(0));   
+    m_Angle2.push_back(clitk::deg2rad(90));
+    break;
+  case SupTo:    
+    m_Angle1.push_back(clitk::deg2rad(0));   
+    m_Angle2.push_back(clitk::deg2rad(-90));
+    break;
+  case Angle:  break;
   }
   /*         A1   A2
              Left      0    0
@@ -160,106 +246,187 @@ SetOrientationType(OrientationTypeEnumeration orientation)
 //--------------------------------------------------------------------
 
 
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+PrintOptions() 
+{
+  DD((int)this->GetBackgroundValue());
+  DD((int)this->GetObjectBackgroundValue());
+  DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
+  DD(this->GetIntermediateSpacingFlag());
+  DD(this->GetIntermediateSpacing());
+  DD(this->GetFuzzyThreshold());
+  DD(this->GetAutoCropFlag());
+  DD(this->GetInverseOrientationFlag());
+  DD(this->GetRemoveObjectFlag());
+  DD(this->GetCombineWithOrFlag());
+}
+//--------------------------------------------------------------------
+
+
 //--------------------------------------------------------------------
 template <class ImageType>
 void 
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 GenerateData() 
 {
+  if (GetNumberOfAngles() <1) {
+    clitkExceptionMacro("Add at least one orientation type");
+  }  
+
+  if (GetVerboseOptionFlag()) {
+    for(int i=0; i<GetNumberOfAngles(); i++) {
+        std::cout << "Orientation    \t:" << GetOrientationTypeString(i) << std::endl;
+    }
+    std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
+    std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
+    std::cout << "AutoCrop       \t:" << GetAutoCropFlag() << std::endl;
+    std::cout << "InverseOrient  \t:" << GetInverseOrientationFlag() << std::endl;
+    std::cout << "RemoveObject   \t:" << GetRemoveObjectFlag() << std::endl;
+    std::cout << "CombineWithOr  \t:" << GetCombineWithOrFlag() << std::endl;
+  }
+
   // Get input pointer
   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
-
-  //--------------------------------------------------------------------
-  //--------------------------------------------------------------------
   static const unsigned int dim = ImageType::ImageDimension;
-  StartNewStep("Initial resample");  
-  // Step 1 : resample
-  if (m_ResampleBeforeRelativePositionFilter) {
-    typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
-    typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
-    resampleFilter->SetInput(object);
-    resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
-    resampleFilter->SetGaussianFilteringEnabled(false);
-    //    resampleFilter->SetVerboseOptions(true);
-    resampleFilter->Update();
-    working_image = resampleFilter->GetOutput();
-  }
-  else {
-    working_image = object;
-  }
-  StopCurrentStep<ImageType>(working_image);
 
   // Step 2: object pad to input image -> we want to compute the
   // relative position for each point belonging to the input image
   // domain, so we have to extend (pad) the object image to fit the
   // domain size
+  working_image = object;
   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
-    StartNewStep("Pad object to image size");  
-    typename ImageType::Pointer output = ImageType::New();
-    SizeType size;
-    for(unsigned int i=0; i<dim; i++) {
-      size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
-    }
+    StartNewStep("Pad (resize) object to input size");  
+
+    if (0) { // OLD VERSION (TO REMOVE)
+      StartNewStep("Pad object to image size");  
+      typename ImageType::Pointer output = ImageType::New();
+      SizeType size;
+      for(unsigned int i=0; i<dim; i++) {
+        size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
+                         input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
+      }
 
-    // The index of the input is not necessarily zero, so we have to
-    // take it into account (not done)
-    RegionType region;
-    IndexType index = input->GetLargestPossibleRegion().GetIndex();
-    region.SetSize(size);
-    for(unsigned int i=0; i<dim; i++) {
-      if (index[i] != 0) {
-        std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
-        exit(0);
+      // The index of the input is not necessarily zero, so we have to
+      // take it into account (not done)
+      RegionType region;
+      IndexType index = input->GetLargestPossibleRegion().GetIndex();
+      region.SetSize(size);
+      for(unsigned int i=0; i<dim; i++) {
+        if (index[i] != 0) {
+          std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
+          exit(0);
+        }
       }
+      // output->SetLargestPossibleRegion(region);
+      output->SetRegions(region);
+      output->SetSpacing(working_image->GetSpacing());    
+      PointType origin = input->GetOrigin();
+      for(unsigned int i=0; i<dim; i++) {
+        origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
+      }
+      output->SetOrigin(origin);
+      //    output->SetOrigin(input->GetOrigin());
+
+      output->Allocate();
+      output->FillBuffer(m_BackgroundValue);
+      typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
+      // typename PasteFilterType::InputImageIndexType index;
+      for(unsigned int i=0; i<dim; i++) {
+        index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
+          + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
+      }
+      padFilter->SetSourceImage(working_image);
+      padFilter->SetDestinationImage(output);
+      padFilter->SetDestinationIndex(index);
+      padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
+      padFilter->Update();
+      working_image = padFilter->GetOutput();
     }
-    // output->SetLargestPossibleRegion(region);
-    output->SetRegions(region);
-    output->SetSpacing(working_image->GetSpacing());    
-    PointType origin = input->GetOrigin();
-    for(unsigned int i=0; i<dim; i++) {
-      origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
-    }
-    output->SetOrigin(origin);
-    //    output->SetOrigin(input->GetOrigin());
-
-    output->Allocate();
-    output->FillBuffer(m_BackgroundValue);
-    typename PadFilterType::Pointer padFilter = PadFilterType::New();
-    // typename PadFilterType::InputImageIndexType index;
-    for(unsigned int i=0; i<dim; i++) {
-      index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
-    + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
-    }
-    padFilter->SetSourceImage(working_image);
-    padFilter->SetDestinationImage(output);
-    padFilter->SetDestinationIndex(index);
-    padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
-    padFilter->Update();
-    working_image = padFilter->GetOutput();
+
+    // Resize object like input
+    working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
     StopCurrentStep<ImageType>(working_image);
   }
-  else {
-    // DD("[debug] RelPos : same size and spacing : no padding");
+
+  //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // Step 1 : resample
+  if (m_IntermediateSpacingFlag) {
+    StartNewStep("Resample object to intermediate spacing");  
+    typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
+    typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
+    resampleFilter->SetInput(working_image);
+    resampleFilter->SetDefaultPixelValue(0);
+    resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
+    resampleFilter->SetGaussianFilteringEnabled(false);
+    //    resampleFilter->SetVerboseOptions(true);
+    resampleFilter->Update();
+    working_image = resampleFilter->GetOutput();
+    StopCurrentStep<ImageType>(working_image);
   }
+
   // Keep object image (with resampline and pad)
   object_resampled = working_image;
- //  StopCurrentStep<ImageType>(working_image);
 
   // Step 3: compute rel pos in object
   StartNewStep("Relative Position Map");  
   typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->SetInput(working_image);
-  relPosFilter->SetAlpha1(m_Angle1); // xy plane
-  relPosFilter->SetAlpha2(m_Angle2);
-  relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
-  relPosFilter->SetFast(true);
-  relPosFilter->SetRadius(1); // seems sufficient in this cas
-  // relPosFilter->SetVerboseProgress(true);
-  relPosFilter->Update();
-  relPos = relPosFilter->GetOutput();
+  typename RelPosFilterType::Pointer relPosFilter;
+
+  for(int i=0; i<GetNumberOfAngles(); i++) {
+    // Compute fuzzy map
+    relPosFilter = RelPosFilterType::New();
+    relPosFilter->SetInput(working_image);
+    relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
+    relPosFilter->SetAlpha2(m_Angle2[i]);
+    relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
+    relPosFilter->SetFast(true);
+    relPosFilter->SetRadius(1); // seems sufficient in this case
+    // relPosFilter->SetVerboseProgress(true);
+    relPosFilter->Update();
+    relPos = relPosFilter->GetOutput();
+
+    if (GetNumberOfAngles() != 1) {
+      // Creation of the first m_FuzzyMap
+      if (i==0) {
+        m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
+        m_FuzzyMap->FillBuffer(0.0);
+      }
+      
+      // Add to current fuzzy map
+      typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
+      typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
+      addFilter->SetInput1(m_FuzzyMap);
+      addFilter->SetInput2(relPos);
+      addFilter->Update();
+      m_FuzzyMap = addFilter->GetOutput();
+    }
+    else m_FuzzyMap = relPos;
+  }
+
+  // Divide by the number of relpos
+  if (GetNumberOfAngles() != 1) {
+#if ITK_VERSION_MAJOR >= 4
+    typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
+    typename DivideFilter::Pointer divideFilter = DivideFilter::New();
+    divideFilter->SetConstant2(GetNumberOfAngles());
+#else
+    typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
+    typename DivideFilter::Pointer divideFilter = DivideFilter::New();
+    divideFilter->SetConstant(GetNumberOfAngles());
+#endif
+    divideFilter->SetInput(m_FuzzyMap);
+    divideFilter->Update();
+    m_FuzzyMap = divideFilter->GetOutput();
+  }
+
+  relPos = m_FuzzyMap;
   StopCurrentStep<FloatImageType>(relPos);
+  if (GetFuzzyMapOnlyFlag()) return;
                
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
@@ -297,8 +464,8 @@ GenerateData()
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   // Step 5: resample to initial spacing
-  if (m_ResampleBeforeRelativePositionFilter) {
-    StartNewStep("Resample to get the same sampling than input");
+  if (m_IntermediateSpacingFlag) {
+    StartNewStep("Resample to come back to the same sampling than input");
     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
     resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
@@ -308,12 +475,13 @@ GenerateData()
     // resampleFilter->SetVerboseOptions(true);
     resampleFilter->Update();
     working_image = resampleFilter->GetOutput();
-  StopCurrentStep<ImageType>(working_image);
+    StopCurrentStep<ImageType>(working_image);
   }
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
+  //if (!HaveSameSizeAndSpacing(working_image, input)) {
   if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
     StartNewStep("Pad to get the same size than input");
     typename ImageType::Pointer temp = ImageType::New();
@@ -321,7 +489,7 @@ GenerateData()
     temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
     temp->Allocate();
     temp->FillBuffer(m_BackgroundValue); 
-    typename PadFilterType::Pointer padFilter2 = PadFilterType::New();
+    typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
     padFilter2->SetSourceImage(working_image);
     padFilter2->SetDestinationImage(temp);
     padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
@@ -346,25 +514,53 @@ GenerateData()
   combineFilter->SetForegroundValue(m_BackgroundValue+1);
   combineFilter->SetInput1(input);
   combineFilter->SetInput2(working_image);
-  combineFilter->SetOperationType(BoolFilterType::And);
-  combineFilter->InPlaceOn();
+  if (GetInverseOrientationFlag())
+    combineFilter->SetOperationType(BoolFilterType::AndNot);
+  else {
+    if (GetCombineWithOrFlag())
+      combineFilter->SetOperationType(BoolFilterType::Or);
+    else
+      combineFilter->SetOperationType(BoolFilterType::And);
+  }
+  combineFilter->InPlaceOff(); // Do not modify initial input (!)
   combineFilter->Update(); 
   working_image = combineFilter->GetOutput();
- 
-  combineFilter = BoolFilterType::New();
-  combineFilter->SetInput1(working_image);
-  combineFilter->SetInput2(object);
-  combineFilter->SetOperationType(BoolFilterType::AndNot);
-  combineFilter->InPlaceOn();
-  combineFilter->Update(); 
 
-  working_image = combineFilter->GetOutput();
+  // Remove (if needed the object from the support)
+  if (GetRemoveObjectFlag()) {
+    combineFilter = BoolFilterType::New();
+    combineFilter->SetInput1(working_image);
+    combineFilter->SetInput2(object);
+    combineFilter->SetOperationType(BoolFilterType::AndNot);
+    combineFilter->InPlaceOn();
+    combineFilter->Update(); 
+    working_image = combineFilter->GetOutput();
+  }
+  // In the other case, we must *add* the initial object to keep it
+  // but not more than the initial support
+  else { 
+    combineFilter = BoolFilterType::New();
+    combineFilter->SetInput1(working_image);
+    combineFilter->SetInput2(object);
+    combineFilter->SetOperationType(BoolFilterType::Or);
+    combineFilter->InPlaceOn();
+    combineFilter->Update(); 
+    working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
+    combineFilter = BoolFilterType::New();
+    combineFilter->SetInput1(working_image);
+    combineFilter->SetInput2(input);
+    combineFilter->SetOperationType(BoolFilterType::And);
+    combineFilter->InPlaceOn();
+    combineFilter->Update(); 
+    working_image = combineFilter->GetOutput();
+  }
+
   StopCurrentStep<ImageType>(working_image);
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   // Step 7: autocrop
-  if (GetAutoCrop()) {
+  if (GetAutoCropFlag()) {
     StartNewStep("Final AutoCrop");
     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
@@ -380,7 +576,7 @@ GenerateData()
   
   // Final Step -> set output
   this->SetNthOutput(0, working_image);
-  return;
+  //  this->GraftOutput(working_image);
 }
 //--------------------------------------------------------------------