diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctionTraits.h b/Modules/Filtering/ImageManipulation/include/otbFunctionTraits.h
deleted file mode 100644
index c28d816833b53993b489eb68447e5ccfc24f7e6f..0000000000000000000000000000000000000000
--- a/Modules/Filtering/ImageManipulation/include/otbFunctionTraits.h
+++ /dev/null
@@ -1,312 +0,0 @@
-//--------------------------------------
-// utils/traits: Additional type traits
-//--------------------------------------
-//
-//          Copyright kennytm (auraHT Ltd.) 2011.
-// Distributed under the Boost Software License, Version 1.0.
-//    (See accompanying file doc/LICENSE_1_0.txt or copy at
-//          http://www.boost.org/LICENSE_1_0.txt)
-
-/**
-
-``<utils/traits.hpp>`` --- Additional type traits
-=================================================
-
-This module provides additional type traits and related functions, missing from
-the standard library.
-
-*/
-
-#ifndef TRAITS_HPP_9ALQFEFX7TO
-#define TRAITS_HPP_9ALQFEFX7TO 1
-
-#include <cstdlib>
-#include <tuple>
-#include <functional>
-#include <type_traits>
-
-namespace FunctionTraits {
-
-/**
-.. macro:: DECLARE_HAS_TYPE_MEMBER(member_name)
-
-    This macro declares a template ``has_member_name`` which will check whether
-    a type member ``member_name`` exists in a particular type.
-
-    Example::
-
-        DECLARE_HAS_TYPE_MEMBER(result_type)
-
-        ...
-
-        printf("%d\n", has_result_type< std::plus<int> >::value);
-        // ^ prints '1' (true)
-        printf("%d\n", has_result_type< double(*)() >::value);
-        // ^ prints '0' (false)
-*/
-#define DECLARE_HAS_TYPE_MEMBER(member_name)                            \
-  template <typename, typename = void>                                  \
-  struct has_##member_name                                              \
-  { enum { value = false }; };                                          \
-  template <typename T>                                                 \
-  struct has_##member_name<T, typename std::enable_if<sizeof(typename T::member_name)||true>::type> \
-  { enum { value = true }; };
-
-/**
-.. type:: struct utils::function_traits<F>
-
-    Obtain compile-time information about a function object *F*.
-
-    This template currently supports the following types:
-
-    * Normal function types (``R(T...)``), function pointers (``R(*)(T...)``)
-      and function references (``R(&)(T...)`` and ``R(&&)(T...)``).
-    * Member functions (``R(C::*)(T...)``)
-    * ``std::function<F>``
-    * Type of lambda functions, and any other types that has a unique
-      ``operator()``.
-    * Type of ``std::mem_fn`` (only for GCC's libstdc++ and LLVM's libc++).
-      Following the C++ spec, the first argument will be a raw pointer.
-*/
-template <typename T>
-struct function_traits
-  : public function_traits<decltype(&T::operator())>
-{};
-
-namespace xx_impl
-{
-template <typename C, typename R, typename... A>
-struct memfn_type
-{
-  typedef typename std::conditional<
-    std::is_const<C>::value,
-    typename std::conditional<
-      std::is_volatile<C>::value,
-      R (C::*)(A...) const volatile,
-      R (C::*)(A...) const
-      >::type,
-    typename std::conditional<
-      std::is_volatile<C>::value,
-      R (C::*)(A...) volatile,
-      R (C::*)(A...)
-      >::type
-    >::type type;
-};
-}
-
-template <typename ReturnType, typename... Args>
-struct function_traits<ReturnType(Args...)>
-{
-  /**
-    .. type:: type result_type
-
-        The type returned by calling an instance of the function object type *F*.
-    */
-  typedef ReturnType result_type;
-
-  /**
-    .. type:: type function_type
-
-        The function type (``R(T...)``).
-    */
-  typedef ReturnType function_type(Args...);
-
-  /**
-    .. type:: type member_function_type<OwnerType>
-
-        The member function type for an *OwnerType* (``R(OwnerType::*)(T...)``).
-    */
-  template <typename OwnerType>
-  using member_function_type = typename xx_impl::memfn_type<
-    typename std::remove_pointer<typename std::remove_reference<OwnerType>::type>::type,
-    ReturnType, Args...
-    >::type;
-
-  /**
-    .. data:: static const size_t arity
-
-        Number of arguments the function object will take.
-    */
-  enum { arity = sizeof...(Args) };
-
-  /**
-    .. type:: type arg<n>::type
-
-        The type of the *n*-th argument.
-    */
-  template <size_t i>
-  struct arg
-  {
-    typedef typename std::tuple_element<i, std::tuple<Args...>>::type type;
-  };
-};
-
-template <typename ReturnType, typename... Args>
-struct function_traits<ReturnType(*)(Args...)>
-  : public function_traits<ReturnType(Args...)>
-{};
-
-template <typename ClassType, typename ReturnType, typename... Args>
-struct function_traits<ReturnType(ClassType::*)(Args...)>
-  : public function_traits<ReturnType(Args...)>
-{
-  typedef ClassType& owner_type;
-};
-
-template <typename ClassType, typename ReturnType, typename... Args>
-struct function_traits<ReturnType(ClassType::*)(Args...) const>
-  : public function_traits<ReturnType(Args...)>
-{
-  typedef const ClassType& owner_type;
-};
-
-template <typename ClassType, typename ReturnType, typename... Args>
-struct function_traits<ReturnType(ClassType::*)(Args...) volatile>
-  : public function_traits<ReturnType(Args...)>
-{
-  typedef volatile ClassType& owner_type;
-};
-
-template <typename ClassType, typename ReturnType, typename... Args>
-struct function_traits<ReturnType(ClassType::*)(Args...) const volatile>
-  : public function_traits<ReturnType(Args...)>
-{
-  typedef const volatile ClassType& owner_type;
-};
-
-template <typename FunctionType>
-struct function_traits<std::function<FunctionType>>
-  : public function_traits<FunctionType>
-{};
-
-#if defined(_GLIBCXX_FUNCTIONAL)
-#define MEM_FN_SYMBOL_XX0SL7G4Z0J std::_Mem_fn
-#elif defined(_LIBCPP_FUNCTIONAL)
-#define MEM_FN_SYMBOL_XX0SL7G4Z0J std::__mem_fn
-#endif
-
-#ifdef MEM_FN_SYMBOL_XX0SL7G4Z0J
-
-template <typename R, typename C>
-struct function_traits<MEM_FN_SYMBOL_XX0SL7G4Z0J<R C::*>>
-  : public function_traits<R(C*)>
-{};
-template <typename R, typename C, typename... A>
-struct function_traits<MEM_FN_SYMBOL_XX0SL7G4Z0J<R(C::*)(A...)>>
-  : public function_traits<R(C*, A...)>
-{};
-template <typename R, typename C, typename... A>
-struct function_traits<MEM_FN_SYMBOL_XX0SL7G4Z0J<R(C::*)(A...) const>>
-  : public function_traits<R(const C*, A...)>
-{};
-template <typename R, typename C, typename... A>
-struct function_traits<MEM_FN_SYMBOL_XX0SL7G4Z0J<R(C::*)(A...) volatile>>
-  : public function_traits<R(volatile C*, A...)>
-{};
-template <typename R, typename C, typename... A>
-struct function_traits<MEM_FN_SYMBOL_XX0SL7G4Z0J<R(C::*)(A...) const volatile>>
-  : public function_traits<R(const volatile C*, A...)>
-{};
-
-#undef MEM_FN_SYMBOL_XX0SL7G4Z0J
-#endif
-
-template <typename T>
-struct function_traits<T&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<const T&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<volatile T&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<const volatile T&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<T&&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<const T&&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<volatile T&&> : public function_traits<T> {};
-template <typename T>
-struct function_traits<const volatile T&&> : public function_traits<T> {};
-
-
-#define FORWARD_RES_8QR485JMSBT                                         \
-  typename std::conditional<                                            \
-                                                                        std::is_lvalue_reference<R>::value, \
-                                                                        T&, \
-                                                                        typename std::remove_reference<T>::type&& \
-                                                                        >::type
-
-/**
-.. function:: auto utils::forward_like<Like, T>(T&& t) noexcept
-
-    Forward the reference *t* like the type of *Like*. That means, if *Like* is
-    an lvalue (reference), this function will return an lvalue reference of *t*.
-    Otherwise, if *Like* is an rvalue, this function will return an rvalue
-    reference of *t*.
-
-    This is mainly used to propagate the expression category (lvalue/rvalue) of
-    a member of *Like*, generalizing ``std::forward``.
-*/
-template <typename R, typename T>
-FORWARD_RES_8QR485JMSBT forward_like(T&& input) noexcept
-{
-  return static_cast<FORWARD_RES_8QR485JMSBT>(input);
-}
-
-#undef FORWARD_RES_8QR485JMSBT
-
-/**
-.. type:: struct utils::copy_cv<From, To>
-
-    Copy the CV qualifier between the two types. For example,
-    ``utils::copy_cv<const int, double>::type`` will become ``const double``.
-*/
-template <typename From, typename To>
-struct copy_cv
-{
-private:
-  typedef typename std::remove_cv<To>::type raw_To;
-  typedef typename std::conditional<std::is_const<From>::value,
-                                    const raw_To, raw_To>::type const_raw_To;
-public:
-  /**
-    .. type:: type type
-
-        Result of cv-copying.
-    */
-  typedef typename std::conditional<std::is_volatile<From>::value,
-                                    volatile const_raw_To, const_raw_To>::type type;
-};
-
-/**
-.. type:: struct utils::pointee<T>
-
-    Returns the type by derefering an instance of *T*. This is a generalization
-    of ``std::remove_pointer``, that it also works with iterators.
-*/
-template <typename T>
-struct pointee
-{
-  /**
-    .. type:: type type
-
-        Result of dereferencing.
-    */
-  typedef typename std::remove_reference<decltype(*std::declval<T>())>::type type;
-};
-
-/**
-.. function:: std::add_rvalue_reference<T>::type utils::rt_val<T>() noexcept
-
-    Returns a value of type *T*. It is guaranteed to do nothing and will not
-    throw a compile-time error, but using the returned result will cause
-    undefined behavior.
-*/
-template <typename T>
-typename std::add_rvalue_reference<T>::type rt_val() noexcept
-{
-  return std::move(*static_cast<T*>(nullptr));
-}
-
-}
-#endif
diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
index 69766a8b319d76704f5839c8e8749346cdca4c79..184f7c59e6b1d03fc1089292a176afef46a7044d 100644
--- a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
+++ b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
@@ -21,175 +21,27 @@
 #ifndef otbFunctorImageFilter_h
 #define otbFunctorImageFilter_h
 
-#include "itkImageToImageFilter.h"
 #include "otbVariadicInputsImageFilter.h"
 #include "otbImage.h"
 #include "otbVectorImage.h"
-#include "itkConstNeighborhoodIterator.h"
-#include "itkImageScanlineIterator.h"
-#include "itkImageRegionConstIterator.h"
-#include "itkProcessObject.h"
+
 #include <type_traits>
 #include <utility>
-#include "otbFunctionTraits.h"
-#include "itkDefaultConvertPixelTraits.h"
 
 namespace otb
 {
 
-namespace functor_filter_details
-{
-// Variadic SetRequestedRegion
-
-// This function sets the requested region for one image
-template<class T> int SetInputRequestedRegion(const T * img, const itk::ImageRegion<2> & region, const itk::Size<2>& radius)
-{
-  auto currentRegion = region;
-  currentRegion.PadByRadius(radius);
-
-  // The ugly cast in all ITK filters
-  T * nonConstImg = const_cast<T*>(img);
-  
-  if(currentRegion.Crop(img->GetLargestPossibleRegion()))
-    {
-    nonConstImg->SetRequestedRegion(currentRegion);
-    return 0;
-    }
-  else
-    {
-    nonConstImg->SetRequestedRegion(currentRegion);
-        
-    // build an exception
-    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
-    std::ostringstream msg;
-    msg << "::SetInputRequestedRegion<>()";
-    e.SetLocation(msg.str());
-    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
-    e.SetDataObject(nonConstImg);
-    throw e;
-    }
-
-  return 0;
-}
-
-template <class Tuple, size_t...Is> auto SetInputRequestedRegionsImpl(Tuple & t, const itk::ImageRegion<2> & region, std::index_sequence<Is...>,const itk::Size<2> & radius)
-{
-  return std::make_tuple(SetInputRequestedRegion(std::get<Is>(t),region,radius)...);
-}
-
-template <typename... T> auto SetInputRequestedRegions(std::tuple<T...> && t,const itk::ImageRegion<2> & region, const itk::Size<2> & radius)
-{
-  return SetInputRequestedRegionsImpl(t,region,std::make_index_sequence<sizeof...(T)>{},radius);
-}
-
-template <class Tuple, size_t...Is> auto GetNumberOfComponentsPerInputImpl(Tuple & t, std::index_sequence<Is...>)
-{
-  return std::array<size_t,sizeof...(Is)>{{std::get<Is>(t)->GetNumberOfComponentsPerPixel()...}};
-}
-
-template <typename ...T> auto GetNumberOfComponentsPerInput(std::tuple<T...> & t)
-{
-  return GetNumberOfComponentsPerInputImpl(t, std::make_index_sequence<sizeof...(T)>{});
-}
-
-template <typename N> struct MakeIterator {};
-
-template <> struct MakeIterator<std::false_type>
-{
-  template <class T> static auto Make(const T * img, const itk::ImageRegion<2> & region, const itk::Size<2>&)
-  {
-    itk::ImageRegionConstIterator<T> it(img,region);
-    return it;
-  }
-};
-
-template <> struct MakeIterator<std::true_type>
-{
-  template <class T> static auto Make(const T * img, const itk::ImageRegion<2> & region, const itk::Size<2>& radius)
-  {
-    itk::ConstNeighborhoodIterator<T> it(radius,img,region);
-    return it;
-  }
-};
-
-
-template <class TNeigh, class Tuple, size_t...Is> auto MakeIteratorsImpl(const Tuple& t, const itk::ImageRegion<2> & region, const itk::Size<2> & radius, std::index_sequence<Is...>, TNeigh)
-{
-  return std::make_tuple(MakeIterator<typename std::tuple_element<Is,TNeigh>::type >::Make(std::get<Is>(t),region,radius)...);
-}
-
-template<class TNeigh, typename... T> auto MakeIterators(std::tuple<T...> &&t,const itk::ImageRegion<2> & region, const itk::Size<2> & radius, TNeigh n)
-  {
-    return MakeIteratorsImpl(t,region,radius,std::make_index_sequence<sizeof...(T)>{},n);
-  }
-
-// Variadic call of operator from iterator tuple
-template <typename T> struct GetProxy{};
-
-
-template <typename T> struct GetProxy<itk::ImageRegionConstIterator<T> >
-{
-  static auto Get(const itk::ImageRegionConstIterator<T> & t)
-{
-  return t.Get();
-}
-};
-
-template <typename T> struct GetProxy<itk::ConstNeighborhoodIterator<T> >
-{
-  static auto Get(const itk::ConstNeighborhoodIterator<T> & t)
-{
-  return t.GetNeighborhood();
-}
-};
-
-template <class Tuple, class Oper, size_t...Is> auto CallOperatorImpl(Tuple& t, const Oper & oper,std::index_sequence<Is...>)
-{
-  return oper(GetProxy<typename std::remove_reference<decltype(std::get<Is>(t))>::type>::Get(std::get<Is>(t))...);
-}
-
-template <class Oper, typename ... Args> auto CallOperator(const Oper& oper, std::tuple<Args...> & t)
-{
-  return CallOperatorImpl(t,oper,std::make_index_sequence<sizeof...(Args)>{});
-}
-
-// Variadic move of iterators
-template<class Tuple,size_t...Is> auto MoveIteratorsImpl(Tuple & t, std::index_sequence<Is...>)
-{
-  return std::make_tuple(++(std::get<Is>(t) )...);
-}
-
-template<typename ... Args> void MoveIterators(std::tuple<Args...> & t)
-{
-  MoveIteratorsImpl(t,std::make_index_sequence<sizeof...(Args)>{});
-}
-
-
-// Default implementation does nothing
-template <class F, class O, size_t N> struct NumberOfOutputComponents
-{};
-
-template <class F, class T, size_t N> struct NumberOfOutputComponents<F,otb::Image<T>,N>
-  {
-  // We can not be here if output type is VectorImage
-    static void Set(const F&, otb::Image<T> *, std::array<size_t,N>){ std::cout<<"Default set"<<std::endl;}
-};
-
-// O is a VectorImage AND F has a fixed OuptutSize static constexrp size_t;
-template <class F, class T, size_t N> struct NumberOfOutputComponents<F,otb::VectorImage<T>,N>
-{
-  static void Set(const F & f, otb::VectorImage<T> * outputImage, std::array<size_t,N> inNbBands)
-  {
-    std::cout<<"Use OutputSize to set number of output components"<<std::endl;
-    outputImage->SetNumberOfComponentsPerPixel(f.OutputSize(inNbBands));
-  }
-};
-
-} // end namespace functor_filter_details
-
-
+/**  
+ * \struct IsNeighborhood 
+ * Struct testing if T is a neighborhood
+ * Provides:
+ * - ValueType type set to false_type or true_type
+ * - value set to true or false
+ * - PixelType type to the underlying pixel type
+ */
 template <class T, class Enable = void> struct IsNeighborhood{};
 
+/// Partial specialisation for scalar types
 template <class T> struct IsNeighborhood<T,typename std::enable_if<std::is_scalar<typename std::remove_reference<typename std::remove_cv<T>::type>::type>::value >::type>
 {
   using ValueType = std::false_type;
@@ -197,6 +49,7 @@ template <class T> struct IsNeighborhood<T,typename std::enable_if<std::is_scala
   using PixelType = T;
 };
 
+/// Partial specialisation for itk::Neighborhood<T>
 template <class T> struct IsNeighborhood<itk::Neighborhood<T>>
 {
   using ValueType = std::true_type;
@@ -204,6 +57,7 @@ template <class T> struct IsNeighborhood<itk::Neighborhood<T>>
   using PixelType = T;
 };
 
+/// Partial specialisation for const itk::Neighborhood<T> &
 template <class T> struct IsNeighborhood<const itk::Neighborhood<T>&>
 {
   using ValueType = std::true_type;
@@ -211,6 +65,7 @@ template <class T> struct IsNeighborhood<const itk::Neighborhood<T>&>
   using PixelType = T;
 };
 
+/// Partial specialisation for itk::VariableLengthVector<T>
 template <class T> struct IsNeighborhood<itk::VariableLengthVector<T>>
 {
   using ValueType = std::false_type;
@@ -218,7 +73,7 @@ template <class T> struct IsNeighborhood<itk::VariableLengthVector<T>>
   using PixelType = itk::VariableLengthVector<T>;
 };
 
-
+/// Partial specialisation for const itk::VariableLengthVector<T> &
 template <class T> struct IsNeighborhood<const itk::VariableLengthVector<T>&>
 {
   using ValueType = std::false_type;
@@ -226,13 +81,22 @@ template <class T> struct IsNeighborhood<const itk::VariableLengthVector<T>&>
   using PixelType = itk::VariableLengthVector<T>;
 };
 
-
-
-
+/**
+ * \struct InputImageTraits
+ * Struct allowing to derive input image types from operator()
+ * arguments
+ * Defines:
+ * - PixelType type to the underlying pixel type
+ * - ScalarType type to the underlying scalar type
+ * - ImageType type to the mocked up image type
+ *
+ * ImageType = Image<T> for T or itk::Neighborhood<T>
+ * ImageType = VectorImage<T> for itk::VariableLengthVector<T> or itk::Neighborhood<itk::VariableLengthVector<T>>
+ */
 template<typename T>
-struct TInputImage
+struct InputImageTraits
 {
-  using ArgumentType = T; //ArgType<T,N>;
+  using ArgumentType = T;
   using PixelType = typename IsNeighborhood<typename std::remove_cv<typename std::remove_reference< ArgumentType>::type >::type>::PixelType;
   using ScalarType = typename itk::DefaultConvertPixelTraits<PixelType>::ComponentType;
   using ImageType = typename std::conditional<std::is_scalar<PixelType>::value,
@@ -240,11 +104,19 @@ struct TInputImage
                                      otb::VectorImage< ScalarType > >::type;
 };
 
-
+/**
+ * \struct OutputImageTraits
+ * Struct allowing to derive output image type from operator()
+ * Defines:
+ * - ScalarType type to the underlying scalar type
+ * - ImageType type to the mocked up image type
+ * - ImageType = Image<T> for T
+ * - ImageType = VectorImage<T> for itk::VariableLengthVector<T>
+ */
 template<typename T>
-struct TOutputImage
+struct OutputImageTraits
 {
-  using ResultType = T; //FResultType<T>;
+  using ResultType = T;
   using ScalarType = typename itk::DefaultConvertPixelTraits<ResultType>::ComponentType;
   using ImageType =  typename std::conditional<std::is_scalar<ResultType>::value,
                                      typename otb::Image<ScalarType>,
@@ -252,59 +124,80 @@ struct TOutputImage
 };
 
 
-
+/**
+* \struct FunctorFilterSuperclassHelper 
+* Struct allowing to derive the superclass prototype for the
+* FunctorImageFilter class
+* Provides the following:
+* - OutputImageType : type of the output image
+* - FilterType : correct instanciation of VariadicInputsImageFilter from
+* - the operator() prototype
+* - InputHasNeighborhood a tuple of N false_type or true_type to denote
+* - if Ith arg of operator() expects a neighborhood.
+*/
 template <typename T> struct FunctorFilterSuperclassHelper : public FunctorFilterSuperclassHelper<decltype(&T::operator())> {};
 
-
+/// Partial specialisation for R(*)(T...)
 template <typename R, typename... T> struct FunctorFilterSuperclassHelper<R(*)(T...)>
 {
-  using OutputImageType = typename TOutputImage<R>::ImageType;
-  using FilterType = VariadicInputsImageFilter<OutputImageType,typename TInputImage<T>::ImageType...>;
+  using OutputImageType = typename OutputImageTraits<R>::ImageType;
+  using FilterType = VariadicInputsImageFilter<OutputImageType,typename InputImageTraits<T>::ImageType...>;
   using InputHasNeighborhood = std::tuple<typename IsNeighborhood<T>::ValueType...>;
 };
 
+// Partial specialisation for R(C::*)(T...) const
 template <typename C, typename R, typename... T> struct FunctorFilterSuperclassHelper<R(C::*)(T...) const>
 {
-  using OutputImageType = typename TOutputImage<R>::ImageType;
-  using FilterType = VariadicInputsImageFilter<OutputImageType,typename TInputImage<T>::ImageType...>;
+  using OutputImageType = typename OutputImageTraits<R>::ImageType;
+  using FilterType = VariadicInputsImageFilter<OutputImageType,typename InputImageTraits<T>::ImageType...>;
   using InputHasNeighborhood = std::tuple<typename IsNeighborhood<T>::ValueType...>;
 };
 
+// Partial specialisation for R(C::*)(T...)
 template <typename C, typename R, typename... T> struct FunctorFilterSuperclassHelper<R(C::*)(T...)>
 {
-  using OutputImageType = typename TOutputImage<R>::ImageType;
-  using FilterType = VariadicInputsImageFilter<OutputImageType,typename TInputImage<T>::ImageType...>;
+  using OutputImageType = typename OutputImageTraits<R>::ImageType;
+  using FilterType = VariadicInputsImageFilter<OutputImageType,typename InputImageTraits<T>::ImageType...>;
   using InputHasNeighborhood = std::tuple<typename IsNeighborhood<T>::ValueType...>;
 };
 
+/**
+ * This helper method builds a fully functional FunctorImageFilter from a functor instance
+ * 
+ * Functor can be any operator() that matches the following:
+ * - Accepts any number of arguments of T,
+ * (const) itk::VariableLengthVector<T> (&),(const)
+ * itk::Neighborhood<T> (&), (const)
+ * itk::Neighborhood<itk::VariableLengthVector<T>> (&) with T a scalar type
+ * - returns T or itk::VariableLengthVector<T>, with T a scalar type
+ *
+ * The returned filter is ready to use. Inputs can be set through the
+ * SetVInputs() method (see VariadicInputsImageFilter class for
+ * details)
+ * 
+ * \param f the Functor to build the filter from
+ * \param radius The size of neighborhood to use, if there is any
+ * itk::Neighborhood<T> in the operator() arguments.
+ * \return A ready to use OTB filter, which accepts n input image of
+ * type derived from the operator() arguments, and producing an image
+ * correpsonding to the operator() return type.
+ *
+ * Note that this function also works with a lambda as Functor,
+ * provided it returns a scalar type. If your lambda returns a
+ * VariableLengthVector, see the other NewFunctorFilter implementation.
+ */
+template <typename Functor> auto NewFunctorFilter(const Functor& f, itk::Size<2> radius = {{0,0}});
 
-template <typename F> struct NumberOfOutputBandsDecorator : F
-{
-public:
-  NumberOfOutputBandsDecorator(const F t, unsigned int nbComp) : F(t), m_NumberOfOutputComponents(nbComp) {}
-  using F::operator();
-
-  constexpr size_t OutputSize(...) const
-  {
-    return m_NumberOfOutputComponents;
-  }
-
-private:
-  unsigned int m_NumberOfOutputComponents;
-};
 
 /** \class FunctorImageFilter
- * \brief Implements 
+ * \brief TODO
  *
- * This class is 
+ * 
  *
- * \ingroup IntensityImageFilters   Multithreaded
+ * \ingroup IntensityImageFilters   Multithreaded Streamed
  *
  * \ingroup OTBImageManipulation
 */
-template <typename Functor> auto NewFunctorFilter(const Functor& f, itk::Size<2> radius = {{0,0}});
-template <typename Functor> auto NewFunctorFilter(const Functor& f, unsigned int numberOfOuptutBands, itk::Size<2> radius = {{0,0}});
-  
 template <class TFunction>
     class ITK_EXPORT FunctorImageFilter
   : public FunctorFilterSuperclassHelper<TFunction>::FilterType
@@ -323,14 +216,11 @@ public:
   using ProcessObjectType = itk::ProcessObject;
 
   using InputHasNeighborhood = typename FunctorFilterSuperclassHelper<TFunction>::InputHasNeighborhood;
-  
- /** Method for creation by passing the filter type. */
- // static Pointer New(const TFunction& f, itk::Size<2> radius = {{0,0}});
-
+ 
 /** Run-time type information (and related methods). */
   itkTypeMacro(FunctorImageFilter, ImageToImageFilter);
-
- 
+  
+  
 /** Get the functor object.  The functor is returned by reference.
    * (Functors do not have to derive from itk::LightObject, so they do
    * not necessarily have a reference count. So we cannot return a
@@ -350,48 +240,34 @@ public:
     return m_Functor;
   }
 
-/** Set the functor object.  This replaces the current Functor with a
-   * copy of the specified Functor. This allows the user to specify a
-   * functor that has ivars set differently than the default functor.
-   * This method requires an operator!=() be defined on the functor
-   * (or the compiler's default implementation of operator!=() being
-   * appropriate). */
-void SetFunctor(const FunctorType& functor)
-{
-  m_Functor = functor;
-  this->Modified();
-}
-
 private:
-friend auto NewFunctorFilter<TFunction>(const TFunction& f, itk::Size<2> radius);
-
-FunctorImageFilter();
-FunctorImageFilter(const FunctorType& f, itk::Size<2> radius) : m_Functor(f), m_Radius(radius) {};
-FunctorImageFilter(const Self &) ;
-void operator =(const Self&) ;
-~FunctorImageFilter() {}
-
-/** FunctorImageFilter can be implemented as a multithreaded filter.
-   * Therefore, this implementation provides a ThreadedGenerateData() routine
-   * which is called for each processing thread. The output image data is
-   * allocated automatically by the superclass prior to calling
-   * ThreadedGenerateData().  ThreadedGenerateData can only write to the
-   * portion of the output image specified by the parameter
-   * "outputRegionForThread"
-   *
-   * \sa ImageToImageFilter::ThreadedGenerateData(),
-   *     ImageToImageFilter::GenerateData()  */
-    virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId) override;
+  /// Actual creation of the filter is handled by this free function
+  friend auto NewFunctorFilter<TFunction>(const TFunction& f, itk::Size<2> radius);
+
+  /// Constructor of functor filter, will copy the functor
+  FunctorImageFilter(const FunctorType& f, itk::Size<2> radius) : m_Functor(f), m_Radius(radius) {};
+  FunctorImageFilter(const Self &) = delete;
+  void operator =(const Self&) = delete;
+  ~FunctorImageFilter() = default;
+
+/** Overload of ThreadedGenerateData  */
+virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId) override;
 
   /**
    * Pad the input requested region by radius
    */
   virtual void GenerateInputRequestedRegion(void) override;
 
+  /**
+   * Will use the OutputSize() method if
+   */
   virtual void GenerateOutputInformation() override;
-  
+
+
+  // The functor member
   FunctorType m_Functor;
 
+  // Radius if needed
   itk::Size<2> m_Radius;
 };
 
@@ -404,6 +280,23 @@ template <typename Functor> auto NewFunctorFilter(const Functor& f, itk::Size<2>
   return p;
 }
 
+
+template <typename F> struct NumberOfOutputBandsDecorator : F
+{
+public:
+  NumberOfOutputBandsDecorator(const F t, unsigned int nbComp) : F(t), m_NumberOfOutputComponents(nbComp) {}
+  using F::operator();
+
+  constexpr size_t OutputSize(...) const
+  {
+    return m_NumberOfOutputComponents;
+  }
+
+private:
+  unsigned int m_NumberOfOutputComponents;
+};
+
+
 template <typename Functor> auto NewFunctorFilter(const Functor& f, unsigned int numberOfOutputBands, itk::Size<2> radius)
 {
   using FunctorType = NumberOfOutputBandsDecorator<Functor>;
diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
index b22082af9817a4cc5bb54832d2cbdaf7fd92a284..f09f006430306cf7e6699ac209c5f58dfb57e8f5 100644
--- a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
+++ b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
@@ -22,11 +22,164 @@
 #define otbFunctorImageFilter_hxx
 
 #include "otbFunctorImageFilter.h"
-#include "itkImageRegionIterator.h"
 #include "itkProgressReporter.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkDefaultConvertPixelTraits.h"
 
 namespace otb
 {
+namespace functor_filter_details
+{
+// Variadic SetRequestedRegion
+
+// This function sets the requested region for one image
+template<class T> int SetInputRequestedRegion(const T * img, const itk::ImageRegion<2> & region, const itk::Size<2>& radius)
+{
+  auto currentRegion = region;
+  currentRegion.PadByRadius(radius);
+
+  // The ugly cast in all ITK filters
+  T * nonConstImg = const_cast<T*>(img);
+  
+  if(currentRegion.Crop(img->GetLargestPossibleRegion()))
+    {
+    nonConstImg->SetRequestedRegion(currentRegion);
+    return 0;
+    }
+  else
+    {
+    nonConstImg->SetRequestedRegion(currentRegion);
+        
+    // build an exception
+    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
+    std::ostringstream msg;
+    msg << "::SetInputRequestedRegion<>()";
+    e.SetLocation(msg.str());
+    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+    e.SetDataObject(nonConstImg);
+    throw e;
+    }
+
+  return 0;
+}
+
+template <class Tuple, size_t...Is> auto SetInputRequestedRegionsImpl(Tuple & t, const itk::ImageRegion<2> & region, std::index_sequence<Is...>,const itk::Size<2> & radius)
+{
+  return std::make_tuple(SetInputRequestedRegion(std::get<Is>(t),region,radius)...);
+}
+
+template <typename... T> auto SetInputRequestedRegions(std::tuple<T...> && t,const itk::ImageRegion<2> & region, const itk::Size<2> & radius)
+{
+  return SetInputRequestedRegionsImpl(t,region,std::make_index_sequence<sizeof...(T)>{},radius);
+}
+
+template <class Tuple, size_t...Is> auto GetNumberOfComponentsPerInputImpl(Tuple & t, std::index_sequence<Is...>)
+{
+  return std::array<size_t,sizeof...(Is)>{{std::get<Is>(t)->GetNumberOfComponentsPerPixel()...}};
+}
+
+template <typename ...T> auto GetNumberOfComponentsPerInput(std::tuple<T...> & t)
+{
+  return GetNumberOfComponentsPerInputImpl(t, std::make_index_sequence<sizeof...(T)>{});
+}
+
+template <typename N> struct MakeIterator {};
+
+template <> struct MakeIterator<std::false_type>
+{
+  template <class T> static auto Make(const T * img, const itk::ImageRegion<2> & region, const itk::Size<2>&)
+  {
+    itk::ImageRegionConstIterator<T> it(img,region);
+    return it;
+  }
+};
+
+template <> struct MakeIterator<std::true_type>
+{
+  template <class T> static auto Make(const T * img, const itk::ImageRegion<2> & region, const itk::Size<2>& radius)
+  {
+    itk::ConstNeighborhoodIterator<T> it(radius,img,region);
+    return it;
+  }
+};
+
+
+template <class TNeigh, class Tuple, size_t...Is> auto MakeIteratorsImpl(const Tuple& t, const itk::ImageRegion<2> & region, const itk::Size<2> & radius, std::index_sequence<Is...>, TNeigh)
+{
+  return std::make_tuple(MakeIterator<typename std::tuple_element<Is,TNeigh>::type >::Make(std::get<Is>(t),region,radius)...);
+}
+
+template<class TNeigh, typename... T> auto MakeIterators(std::tuple<T...> &&t,const itk::ImageRegion<2> & region, const itk::Size<2> & radius, TNeigh n)
+  {
+    return MakeIteratorsImpl(t,region,radius,std::make_index_sequence<sizeof...(T)>{},n);
+  }
+
+// Variadic call of operator from iterator tuple
+template <typename T> struct GetProxy{};
+
+
+template <typename T> struct GetProxy<itk::ImageRegionConstIterator<T> >
+{
+  static auto Get(const itk::ImageRegionConstIterator<T> & t)
+{
+  return t.Get();
+}
+};
+
+template <typename T> struct GetProxy<itk::ConstNeighborhoodIterator<T> >
+{
+  static auto Get(const itk::ConstNeighborhoodIterator<T> & t)
+{
+  return t.GetNeighborhood();
+}
+};
+
+template <class Tuple, class Oper, size_t...Is> auto CallOperatorImpl(Tuple& t, const Oper & oper,std::index_sequence<Is...>)
+{
+  return oper(GetProxy<typename std::remove_reference<decltype(std::get<Is>(t))>::type>::Get(std::get<Is>(t))...);
+}
+
+template <class Oper, typename ... Args> auto CallOperator(const Oper& oper, std::tuple<Args...> & t)
+{
+  return CallOperatorImpl(t,oper,std::make_index_sequence<sizeof...(Args)>{});
+}
+
+// Variadic move of iterators
+template<class Tuple,size_t...Is> auto MoveIteratorsImpl(Tuple & t, std::index_sequence<Is...>)
+{
+  return std::make_tuple(++(std::get<Is>(t) )...);
+}
+
+template<typename ... Args> void MoveIterators(std::tuple<Args...> & t)
+{
+  MoveIteratorsImpl(t,std::make_index_sequence<sizeof...(Args)>{});
+}
+
+
+// Default implementation does nothing
+template <class F, class O, size_t N> struct NumberOfOutputComponents
+{};
+
+template <class F, class T, size_t N> struct NumberOfOutputComponents<F,otb::Image<T>,N>
+  {
+  // We can not be here if output type is VectorImage
+    static void Set(const F&, otb::Image<T> *, std::array<size_t,N>){}
+};
+
+// O is a VectorImage AND F has a fixed OuptutSize static constexrp size_t;
+template <class F, class T, size_t N> struct NumberOfOutputComponents<F,otb::VectorImage<T>,N>
+{
+  static void Set(const F & f, otb::VectorImage<T> * outputImage, std::array<size_t,N> inNbBands)
+  {
+    outputImage->SetNumberOfComponentsPerPixel(f.OutputSize(inNbBands));
+  }
+};
+
+} // end namespace functor_filter_details
+
+
 
 template <class TFunction>
 FunctorImageFilter<TFunction>::FunctorImageFilter()