diff --git a/app/otbClearCutsDetection.cxx b/app/otbClearCutsDetection.cxx
index 1ce07dbebc405727c39ab469cef9a2589ed715e7..935a8b627d34517e14dfa958d7f2b62df2484bd3 100644
--- a/app/otbClearCutsDetection.cxx
+++ b/app/otbClearCutsDetection.cxx
@@ -22,7 +22,7 @@
 
 // Filters
 #include "itkMaskImageFilter.h"
-#include "otbStatisticsImageCustomFilter.h" // TODO use the otb::StreamingStatisticsFilter
+#include "otbStreamingStatisticsImageFilter.h"
 #include "otbStreamingResampleImageFilter.h"
 #include "otbMultiChannelExtractROI.h"
 #include "otbDeltaNDVILabelerFilter.h"
@@ -72,7 +72,7 @@ public:
       FloatVectorImageType> MaskVectorImageFilterType;
   typedef otb::DeltaNDVILabelerFilter<FloatImageType,
       FloatImageType> NDVILabelImageFilterType;
-  typedef otb::StatisticsImageCustomFilter<FloatImageType> StatsFilterType;
+  typedef otb::StreamingStatisticsImageFilter<FloatImageType> StatsFilterType;
   typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType,
       FloatVectorImageType::InternalPixelType> ExtractROIFilterType;
 
@@ -122,6 +122,8 @@ public:
     MandatoryOff("inamask");
     AddParameter(ParameterType_InputImage,   "mask",  "Input mask (used for both input images)");
     MandatoryOff("mask");
+    AddParameter(ParameterType_Directory, "masksdir", "Vegetation masks directory");
+    MandatoryOff("masksdir");
 
     // Input images band indices
     AddParameter(ParameterType_Int, "nirb", "near infrared band index for input T0 image" );
@@ -145,16 +147,13 @@ public:
     SetDefaultParameterInt("reda",1);
 
     // Output image
-    AddParameter(ParameterType_OutputImage,  "out",   "Output d(NDVI) image");
-    SetDefaultOutputPixelType("out",ImagePixelType_float);
+    AddParameter(ParameterType_OutputImage,  "outdndvi",   "Output d(NDVI) image");
+    MandatoryOff("outdndvi");
+    SetDefaultOutputPixelType("outdndvi",ImagePixelType_float);
 
     AddParameter(ParameterType_OutputImage, "outlabel", "Output labeled image");
     SetDefaultOutputPixelType("outlabel",ImagePixelType_uint8);
 
-    // Vegetation masks directory
-    AddParameter(ParameterType_Directory, "masksdir", "Vegetation masks directory");
-    MandatoryOff("masksdir");
-
     AddRAMParameter();
 
   }
@@ -257,7 +256,8 @@ public:
 
     // Stats filter
     m_StatsFilter = StatsFilterType::New();
-    m_StatsFilter->SetIgnoredInputValue(m_DeltaNDVIFilter->GetFunctor().GetNoDataValue());
+    m_StatsFilter->SetIgnoreUserDefinedValue(true);
+    m_StatsFilter->SetUserIgnoredValue(m_DeltaNDVIFilter->GetFunctor().GetNoDataValue());
 
     // Mask directory
     if (HasValue("masksdir") || HasValue("mask"))
@@ -327,12 +327,19 @@ public:
       m_StatsFilter->SetInput(m_DeltaNDVIFilter->GetOutput());
       }
 
-    // Connect delta NDVI to output, through stats filter
-    SetParameterOutputImage("out", m_StatsFilter->GetOutput());
+    // Compute stats
+    AddProcess(m_StatsFilter->GetStreamer(),"Computing dNDVI statistics");
+    m_StatsFilter->Update();
+
+    if (HasValue("outdndvi"))
+      {
+      // Connect delta NDVI to output
+      SetParameterOutputImage("outdndvi", m_DeltaNDVIFilter->GetOutput());
+      }
 
     // Label image output
     m_NDVILabelFilter = NDVILabelImageFilterType::New();
-    m_NDVILabelFilter->SetInput(m_StatsFilter->GetOutput());
+    m_NDVILabelFilter->SetInput(m_DeltaNDVIFilter->GetOutput());
     m_NDVILabelFilter->SetInputMeanObject(m_StatsFilter->GetMeanOutput());
     m_NDVILabelFilter->SetInputSigmaObject(m_StatsFilter->GetSigmaOutput());
     m_NDVILabelFilter->SetInputNoDataValue(m_DeltaNDVIFilter->GetFunctor().GetNoDataValue());
@@ -341,18 +348,16 @@ public:
   }
 
   ResampleImageFilterType::Pointer m_ResampleFilter;
-  ExtractROIFilterType::Pointer m_ExtractROIFilter;
+  ExtractROIFilterType::Pointer    m_ExtractROIFilter;
   ResampleMaskImageFilterType::Pointer m_MaskResampleFilter;
   DeltaNDVIFilterType::Pointer m_DeltaNDVIFilter;
   MaskImageFilterType::Pointer m_MaskImageFilter;
   MaskVectorImageFilterType::Pointer m_MaskImageFilter0;
   MaskVectorImageFilterType::Pointer m_MaskImageFilter1;
-  NDVILabelImageFilterType::Pointer m_NDVILabelFilter;
+  NDVILabelImageFilterType::Pointer  m_NDVILabelFilter;
   StatsFilterType::Pointer m_StatsFilter;
-
   MaskHandlerType::Pointer m_MaskHandler;
-
-  AndFilterType::Pointer m_AndFilter;
+  AndFilterType::Pointer   m_AndFilter;
 
 };
 }
diff --git a/include/otbDeltaNDVIFunctor.h b/include/otbDeltaNDVIFunctor.h
index 887c9a6ec8bb5ef61dfe2b5149ac2160e0d83028..780f660e5a812b7e02e3c12a9deac5c1ca84831e 100644
--- a/include/otbDeltaNDVIFunctor.h
+++ b/include/otbDeltaNDVIFunctor.h
@@ -567,7 +567,7 @@ public:
     for (unsigned int i = 0 ; i < inputStackedLabelsPixel.Size() ; i++)
       {
       TLabelValue labelValue = inputStackedLabelsPixel[i];
-      if (labelValue != m_NoDataLabel && (label < labelValue || label == this->m_NoDataLabel))
+      if (labelValue != this->m_NoDataLabel && (label < labelValue || label == this->m_NoDataLabel))
         {
         label = labelValue;
         }
diff --git a/include/otbRegionComparator.h b/include/otbRegionComparator.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac8cabb2c68e403ff80351f0b78b4e33cdfb9b65
--- /dev/null
+++ b/include/otbRegionComparator.h
@@ -0,0 +1,229 @@
+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef otbRegionComparator_H_
+#define otbRegionComparator_H_
+
+#include "otbObjectList.h"
+#include "itkContinuousIndex.h"
+#include "itkMetaDataDictionary.h"
+#include "otbMetaDataKey.h"
+#include "itkMetaDataObject.h"
+#include "itkNumericTraits.h"
+
+#include "otbRemoteSensingRegion.h"
+
+namespace otb {
+
+template<class TInputImage1, class TInputImage2=TInputImage1> class RegionComparator{
+public:
+
+  typedef TInputImage1    InputImage1Type;
+  typedef typename InputImage1Type::Pointer   InputImage1PointerType;
+  typedef typename InputImage1Type::RegionType    InputImage1RegionType;
+
+  typedef TInputImage2    InputImage2Type;
+  typedef typename InputImage2Type::Pointer   InputImage2PointerType;
+  typedef typename InputImage2Type::RegionType    InputImage2RegionType;
+
+  typedef typename otb::RemoteSensingRegion<double> RSRegion;
+
+  void SetImage1(InputImage1PointerType im1) {m_InputImage1 = im1;}
+  void SetImage2(InputImage2PointerType im2) {m_InputImage2 = im2;}
+
+  InputImage2RegionType GetOverlapInImage1Indices()
+  {
+    // Largest region of im2 --> region in im1
+    InputImage1RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    InputImage1RegionType region2_in_1;
+    OutputRegionToInputRegion(region2, region2_in_1, m_InputImage2, m_InputImage1);
+
+    // Collision test
+    InputImage2RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    region2_in_1.Crop(region1);
+
+    // Return overlap
+    return region2_in_1;
+  }
+
+  InputImage2RegionType GetOverlapInImage2Indices()
+  {
+    // Largest region of im1 --> region in im2
+    InputImage1RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    InputImage1RegionType region1_in_2;
+    OutputRegionToInputRegion(region1, region1_in_2, m_InputImage1, m_InputImage2);
+
+    // Collision test
+    InputImage2RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    region1_in_2.Crop(region2);
+
+    // Return overlap
+    return region1_in_2;
+  }
+
+  bool DoesOverlap()
+  {
+    // Largest region of im1 --> region in im2
+    InputImage1RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    InputImage1RegionType region1_in_2;
+    OutputRegionToInputRegion(region1, region1_in_2, m_InputImage1, m_InputImage2);
+
+    // Collision test
+    InputImage2RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    return region1_in_2.Crop(region2);
+  }
+
+  /*
+   * Converts sourceRegion of im1 into targetRegion of im2
+   */
+  void OutputRegionToInputRegion(const InputImage1RegionType &sourceRegion, InputImage2RegionType &targetRegion,
+      const InputImage1PointerType &sourceImage, const InputImage2PointerType &targetImage)
+  {
+
+    // Source Region (source image indices)
+    typename InputImage1RegionType::IndexType sourceRegionIndexStart = sourceRegion.GetIndex();
+    typename InputImage1RegionType::IndexType sourceRegionIndexEnd;
+    for(unsigned int dim = 0; dim < InputImage1Type::ImageDimension; ++dim)
+      sourceRegionIndexEnd[dim]= sourceRegionIndexStart[dim] + sourceRegion.GetSize()[dim]-1;
+
+    // Source Region (Geo)
+    typename InputImage1Type::PointType sourceRegionIndexStartGeo, sourceRegionIndexEndGeo;
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexStart, sourceRegionIndexStartGeo);
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexEnd  , sourceRegionIndexEndGeo  );
+
+    // Source Region (target image indices)
+    typename InputImage2RegionType::IndexType targetIndexStart, targetIndexEnd;
+    targetImage->TransformPhysicalPointToIndex(sourceRegionIndexStartGeo, targetIndexStart);
+    targetImage->TransformPhysicalPointToIndex(sourceRegionIndexEndGeo  , targetIndexEnd);
+
+    // Target Region (target image indices)
+    typename InputImage2RegionType::IndexType targetRegionStart;
+    typename InputImage2RegionType::SizeType targetRegionSize;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      targetRegionStart[dim] = std::min(targetIndexStart[dim], targetIndexEnd[dim]);
+      targetRegionSize[dim]  = std::max(targetIndexStart[dim], targetIndexEnd[dim]) - targetRegionStart[dim] + 1;
+      }
+    InputImage2RegionType computedInputRegion(targetRegionStart, targetRegionSize);
+
+    // Avoid extrapolation
+    computedInputRegion.PadByRadius(1);
+
+    // Target Region
+    targetRegion = computedInputRegion;
+
+  }
+
+  /*
+   * Converts a RemoteSensingREgion into a ImageRegion (image1)
+   */
+  InputImage1RegionType RSRegionToImageRegion(const RSRegion rsRegion)
+  {
+    typename itk::ContinuousIndex<double> rsRegionOrigin = rsRegion.GetOrigin();
+    typename itk::ContinuousIndex<double> rsRegionSize = rsRegion.GetSize();
+    typename itk::ContinuousIndex<double> rsRegionEnd;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      rsRegionEnd[dim] = rsRegionOrigin[dim] + rsRegionSize[dim];
+      }
+    typename InputImage1RegionType::IndexType imageRegionStart, imageRegionEnd;
+    m_InputImage1->TransformPhysicalPointToIndex(rsRegionOrigin, imageRegionStart);
+    m_InputImage1->TransformPhysicalPointToIndex(rsRegionEnd, imageRegionEnd);
+
+    // Target Region (target image indices)
+    typename InputImage1RegionType::IndexType targetRegionStart;
+    typename InputImage1RegionType::SizeType targetRegionSize;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      targetRegionStart[dim] = std::min(imageRegionStart[dim], imageRegionEnd[dim]);
+      targetRegionSize[dim]  = std::max(imageRegionStart[dim], imageRegionEnd[dim]) - targetRegionStart[dim] + 1;
+      }
+
+    InputImage1RegionType region(targetRegionStart, targetRegionSize);
+    return region;
+
+  }
+
+  /*
+   * Converts a given region of an image into a remoste sensing region
+   */
+  RSRegion ImageRegionToRSRegion(const InputImage1RegionType &sourceRegion, InputImage1PointerType sourceImage)
+  {
+
+    // Source Region (source image indices)
+    typename InputImage1RegionType::IndexType sourceRegionIndexStart = sourceRegion.GetIndex();
+    typename InputImage1RegionType::IndexType sourceRegionIndexEnd;
+    for(unsigned int dim = 0; dim < InputImage1Type::ImageDimension; ++dim)
+      sourceRegionIndexEnd[dim]= sourceRegionIndexStart[dim] + sourceRegion.GetSize()[dim]-1;
+
+    // Source Region (Geo)
+    typename InputImage1Type::PointType sourceRegionIndexStartGeo, sourceRegionIndexEndGeo;
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexStart, sourceRegionIndexStartGeo);
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexEnd  , sourceRegionIndexEndGeo  );
+
+    // Source Region (Geo min & max)
+    typename itk::ContinuousIndex<double> sourceRegionMin, sourceRegionMax, sourceRegionSize;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      sourceRegionMin[dim] = std::min(sourceRegionIndexStartGeo[dim], sourceRegionIndexEndGeo[dim]);
+      sourceRegionMax[dim] = std::max(sourceRegionIndexStartGeo[dim], sourceRegionIndexEndGeo[dim]);
+      }
+    RSRegion region;
+    sourceRegionSize[0] = sourceRegionMax[0] - sourceRegionMin[0];
+    sourceRegionSize[1] = sourceRegionMax[1] - sourceRegionMin[1];
+    region.SetOrigin(sourceRegionMin);
+    region.SetSize(sourceRegionSize);
+    return region;
+
+  }
+  bool HaveSameProjection()
+  {
+    itk::MetaDataDictionary metaData1 = m_InputImage1 -> GetMetaDataDictionary();
+    itk::MetaDataDictionary metaData2 = m_InputImage2 -> GetMetaDataDictionary();
+
+    if (metaData1.HasKey(otb::MetaDataKey::ProjectionRefKey))
+      {
+      if (metaData2.HasKey(otb::MetaDataKey::ProjectionRefKey))
+        {
+        std::string proj1, proj2;
+        itk::ExposeMetaData<std::string>(metaData1, static_cast<std::string>(otb::MetaDataKey::ProjectionRefKey), proj1);
+        itk::ExposeMetaData<std::string>(metaData2, static_cast<std::string>(otb::MetaDataKey::ProjectionRefKey), proj2);
+        if (proj1.compare(proj2)==0)
+          return true;
+        }
+      }
+    return false;
+
+  }
+
+  bool HaveSameNbOfBands()
+  {
+    return ( (m_InputImage1->GetNumberOfComponentsPerPixel()) == (m_InputImage2->GetNumberOfComponentsPerPixel()) );
+  }
+
+  RSRegion ComputeIntersectingRemoteSensingRegion()
+  {
+    InputImage2RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    InputImage1RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    RSRegion rsr1 = ImageRegionToRSRegion(region1, m_InputImage1);
+    RSRegion rsr2 = ImageRegionToRSRegion(region2, m_InputImage2);
+    rsr1.Crop(rsr2);
+    return rsr1;
+  }
+private:
+
+  InputImage1PointerType m_InputImage1;
+  InputImage2PointerType m_InputImage2;
+
+};
+
+}
+
+#endif /* GTBRegionComparator_H_ */
diff --git a/otb-module.cmake b/otb-module.cmake
index d239b3070dcb4c77d5cb1a74169c8f5a05e02785..f627893b5b03dd881ef5b705a24be66f56a7b452 100644
--- a/otb-module.cmake
+++ b/otb-module.cmake
@@ -2,7 +2,6 @@ set(DOCUMENTATION "Clear Cuts Detection")
 
 otb_module(ClearCutsDetection
   DEPENDS
-	SimpleExtractionTools
 	Mosaic
 	OTBIndices
 	OTBStatistics