diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
index 9125660e60424c20423c71106a5a91a0aac35d62..13dbb661b3795072fe599aeb89e764e251238bb1 100644
--- a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
+++ b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
@@ -39,6 +39,8 @@ 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;
@@ -46,7 +48,7 @@ template<class T> int SetInputRequestedRegion(const T * img, const itk::ImageReg
 
   // The ugly cast in all ITK filters
   T * nonConstImg = const_cast<T*>(img);
-
+  
   if(currentRegion.Crop(img->GetLargestPossibleRegion()))
     {
     nonConstImg->SetRequestedRegion(currentRegion);
@@ -79,12 +81,14 @@ template <typename... T> auto SetInputRequestedRegions(std::tuple<T...> && t,con
   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()...}};
+}
 
-// Variadic creation of iterator tuple
-template <class T> auto MakeIterator(const T * img, const itk::ImageRegion<2> & region)
+template <typename ...T> auto GetNumberOfComponentsPerInput(std::tuple<T...> & t)
 {
-  itk::ImageRegionConstIterator<T> it(img,region);
-  return it;
+  return GetNumberOfComponentsPerInputImpl(t, std::make_index_sequence<sizeof...(T)>{});
 }
 
 template <class T> auto MakeIterator(itk::SmartPointer<T> img, const itk::ImageRegion<2> & region)
@@ -132,6 +136,27 @@ 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;}
+};
+
+// Case 1: 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
 
 
@@ -195,24 +220,24 @@ template <typename C, typename R, typename... T> struct FunctorFilterSuperclassH
   using FilterType = VariadicInputsImageFilter<OutputImageType,typename TInputImage<T>::ImageType...>;
 };
 
-// Default implementation does nothing
-template <class F, class O, class cond = void> struct NumberOfOutputComponents
-{
-  // We can not be here if output type is VectorImage
-  //static_assert(std::is_same<O, otb::VectorImage<typename O::InternalPixelType> >::value,"Return type of Functor is a VariableLenghtVector, add a constexpr size_t OutputSize member");
-  static void Set(const F&, O *){}
-};
-
-// Case 1: O is a VectorImage AND F has a fixed OuptutSize unsigned
-// int constexpr
-template <class F, class T> struct NumberOfOutputComponents<F,T,typename std::enable_if<std::is_same<size_t, decltype(F::OutputSize)>::value>::type >
-{
-static void Set(const F &, otb::VectorImage<T> * outputImage)
-  {
-    std::cout<<"Use OutputSize to set number of output components"<<std::endl;
-    outputImage->SetNumberOfComponentsPerPixel(F::OutputSize);
-  }
-};
+// template <class F> class NumberOfOutputDecorator
+// {
+// public:
+//   NumberOfOutputDecorator(const F & f, size_t n) : m_Functor(f), m_NumberOfOutputs(n) {}
+  
+//   template<typename ...T> auto operator()(T...args)
+//   {
+//     return m_Functor(args...);
+//   }
+  
+//   constexpr size_t OutputSize(...) const
+//   {
+//     return m_NumberOfOutputs;
+//   }
+// private:
+//   F      m_Functor;
+//   size_t m_NumberOfOutputs;
+// };
 
 /** \class FunctorImageFilter
  * \brief Implements 
@@ -241,9 +266,7 @@ public:
   using ProcessObjectType = itk::ProcessObject;
 
  /** Method for creation by passing the filter type. */
-
-
-  static Pointer New(const FunctorType& f) ;
+  static Pointer New(const TFunction& f, itk::Size<2> radius = {{0,0}});
 
 /** Run-time type information (and related methods). */
   itkTypeMacro(FunctorImageFilter, ImageToImageFilter);
@@ -282,7 +305,7 @@ void SetFunctor(const FunctorType& functor)
 
 private:
 FunctorImageFilter();
-FunctorImageFilter(const FunctorType& f) : m_Functor(f) {};
+FunctorImageFilter(const FunctorType& f, itk::Size<2> radius) : m_Functor(f), m_Radius(radius) {};
 FunctorImageFilter(const Self &) ;
 void operator =(const Self&) ;
 ~FunctorImageFilter() {}
@@ -309,6 +332,8 @@ protected:
   virtual void GenerateOutputInformation() override;
   
   FunctorType m_Functor;
+
+  itk::Size<2> m_Radius;
 };
 
 }// namespace otb
diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
index e8b7808c82d640b3d06de7e015ed50ade77fc877..7ac39b47d1c3aed0b29d2aa1d2fda24215485cc3 100644
--- a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
+++ b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
@@ -30,16 +30,13 @@ namespace otb
 
 template <class TFunction>
 FunctorImageFilter<TFunction>::FunctorImageFilter()
-{
-  //  this->SetNumberOfRequiredInputs(m_Functor.GetNumberOfInputs());
-  //  this->InPlaceOff();
-}
+{}
 
 template <class TFunction>
-typename FunctorImageFilter<TFunction>::Pointer 
-FunctorImageFilter<TFunction>::New(const TFunction& f) 
+typename FunctorImageFilter<TFunction>::Pointer
+FunctorImageFilter<TFunction>::New(const TFunction& f, itk::Size<2> radius) 
 {
-  Pointer p = new FunctorImageFilter<TFunction>(f);
+  Pointer p = new FunctorImageFilter<TFunction>(f,radius);
   return p;                                           
 }
 
@@ -63,7 +60,17 @@ template <class TFunction>
 void
 FunctorImageFilter<TFunction>::GenerateOutputInformation()
 {
-  NumberOfOutputComponents<TFunction,OutputImageType>::Set(m_Functor,this->GetOutput());
+  // Call Superclass implementation
+  Superclass::GenerateOutputInformation();
+
+  // Get All variadic inputs
+  auto inputs = this->GetVInputs();
+
+  // Retrieve an array of number of components per input
+  auto inputNbComps = functor_filter_details::GetNumberOfComponentsPerInput(inputs);
+
+  // Call the helper to set the number of components for the output image
+  functor_filter_details::NumberOfOutputComponents<TFunction,OutputImageType,inputNbComps.size()>::Set(m_Functor,this->GetOutput(),inputNbComps);
 }
 
 /**
@@ -74,8 +81,7 @@ void
 FunctorImageFilter<TFunction>
 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
 {
-
-  
+  // Build output iterator
   itk::ImageScanlineIterator<OutputImageType> outIt(this->GetOutput(),outputRegionForThread);
   itk::ProgressReporter p(this,threadId,outputRegionForThread.GetNumberOfPixels());
 
diff --git a/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx b/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx
index 466fc7d0c85353c541ba98e05e6c0a2781cc3904..729155aa0fa1862b19908411f01713e75fa8181a 100644
--- a/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx
+++ b/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx
@@ -20,60 +20,166 @@
 
 #include "itkMacro.h"
 #include "otbFunctorImageFilter.h"
-#include "otbVariadicInputsImageFilter.h"
-
 #include "otbImage.h"
 #include "otbVectorImage.h"
 #include "itkNeighborhood.h"
-#include "itkImageRegionIterator.h"
-#include "itkImageSource.h"
 #include <tuple>
 
+#include <numeric>
+
+// Example functors
 
-struct Funct1
+// N scalar images -> image
+// This functor takes N scalar image (variadic N) and returns the sum
+// of all pixels
+template <typename O, typename ...T> struct VariadicAdd
 {
-  double operator()(double p) const
+  auto operator()(T... ins) const
   {
-    return p;
+    std::vector<O> outVector{static_cast<O>(ins)...};
+
+    return std::accumulate(outVector.begin(), outVector.end(),0);
   }
 };
 
-struct Funct2
+// Redefinition of AddImageFilter, fully variadic
+template<typename TOutputImage,
+         typename ... TInputImage>
+using AddImageFilter= otb::FunctorImageFilter<
+  VariadicAdd<typename TOutputImage::PixelType,
+              typename TInputImage::PixelType...> >;
+
+
+// Temptative make_add_filter function to make things even simpler
+template<typename O, typename ...T> auto make_add_filter(typename otb::Image<T>::Pointer... imgs)
+{
+  using Filter = AddImageFilter<otb::Image<O>,otb::Image<T>...>;
+
+  auto filter =  Filter::New(typename Filter::FunctorType{});
+
+  filter->SetVInputs(imgs...);
+
+  return filter;
+}
+
+// helper function to implement next functor (convert a scalar value
+// to a VariableLengthVector)
+template <typename T> itk::VariableLengthVector<T> toVector(const T & in)
+{
+  itk::VariableLengthVector<T> out;
+  out.SetSize(1);
+  out[0] = in;
+  return out;
+}
+
+// helper function to implement next functor, VariableLengthVectorVersion (returns in)
+template <typename  T> const itk::VariableLengthVector<T> & toVector(const itk::VariableLengthVector<T> & in)
 {
-  itk::VariableLengthVector<double> operator()(double p) const
+  return in;
+}
+
+// helper function to implement next functor, Merge two VariableLengthVector in-place
+template <typename v1, typename v2> void concatenateVectors(v1 & a, const v2 & b)
+{
+  const size_t previousSizeOfA = a.GetSize();
+  
+  a.SetSize(previousSizeOfA+b.GetSize());
+  
+  for(size_t it = 0; it<b.Size();++it)
+    {
+    a[previousSizeOfA+it] = static_cast<typename v1::ValueType>(b[it]);
+    }
+}
+
+// helper function to implement next functor, Merge N VariableLengthVector in-place
+template <typename v1, typename v2, typename ...vn> void concatenateVectors(v1 & a, const v2 & b, const vn&... z)
+{
+  concatenateVectors(a,b);
+  concatenateVectors(a,z...);
+}
+
+// N  images (all types) -> vector image
+// This functor concatenates N images (N = variadic) of type
+// VectorImage and or Image, into a single VectorImage
+template<typename O, typename ...T> struct VariadicConcatenate
+{
+  auto operator()(const T &...  ins) const
   {
-    itk::VariableLengthVector<double> result(OutputSize);
-    result[0] = result[1] = p;
-    return result;
+    itk::VariableLengthVector<O> out;
+    concatenateVectors(out, toVector(ins)...);
+    
+    return out;
   }
 
-  static constexpr size_t OutputSize = 2;
+  // Must define OutputSize because output pixel is vector
+  constexpr size_t OutputSize(const std::array<size_t, sizeof...(T)> inputsNbBands) const
+  {
+    return std::accumulate(inputsNbBands.begin(),inputsNbBands.end(),0);
+  }
 };
 
-using IntImageNeighborhood = itk::Neighborhood<int>;
-struct Funct3
+
+// 1 VectorImage -> 1 VectorImage with a different size depending on a
+// parameter of the functor
+// This Functor 
+template<typename O, typename T> struct BandExtraction
 {
-  double operator()(const IntImageNeighborhood & p) const
+  BandExtraction(unsigned int indices...) : m_Indices({indices})
+  {}
+  
+  // TODO define a constructor to initialize m_Indices  
+  auto operator()(const itk::VariableLengthVector<T> & in) const
   {
-    return static_cast<double>(p.GetCenterValue());
+    itk::VariableLengthVector<O> out(m_Indices.size());
+
+    size_t idx = 0;
+    for(auto v: m_Indices)
+      {
+      out[idx] = static_cast<O>(in[v]);
+      ++idx;
+      }
+
+    return out;
   }
-};
 
-auto Lambda1 = [](double p)
-               {
-                 return p;
-               };
+  // This time OutputSize does not depend on input image size, hence
+  // the ...
+  size_t OutputSize(...) const
+  {
+    return m_Indices.size();
+  }
+  
+  // set of band indices to extract
+  std::set<unsigned int> m_Indices;
+};
 
-struct Funct4
-{
-  itk::VariableLengthVector<double> operator()(double p, const itk::VariableLengthVector<double> & vp) const
+// 1 Image with neighborhood -> 1 Image
+// This Functor computes the mean in neighborhood
+template<typename TOut, typename TIn> struct Median
+{  
+  auto operator()(const itk::Neighborhood<TIn> & in) const
   {
-    itk::VariableLengthVector<double> result(2);
-    result.Fill(p);
-    return result;
+    TOut out(0);
+
+    for(auto it = in.Begin(); it!=in.End();++it)
+      out+=static_cast<TOut>(*it);
+
+    out/=in.Size();
+    
+    return out;
   }
 };
 
+
+
+// // A lambda
+// auto Lambda2 = [](double p)
+//                {
+//                  itk::VariableLengthVector<double> ret(3);
+//                  ret.Fill(p);
+//                  return ret;
+//                };
+
 using OImage = typename otb::Image<int>;
 using Neig = typename itk::Neighborhood<int>;
 
@@ -92,11 +198,7 @@ int otbFunctorImageFilter(int itkNotUsed(argc), char * itkNotUsed(argv) [])
   using ImageType       = Image<double>;
   using RegionType      = typename ImageType::RegionType;
   using SizeType        = typename RegionType::SizeType;
-  using IndexType       = typename RegionType::IndexType;
   
-  using NeighborhoodType = itk::Neighborhood<double,2>;
-  using VectorPixelType  = typename VectorImageType::PixelType;
-
   auto vimage = VectorImageType::New();
   auto image  = ImageType::New();
 
@@ -113,45 +215,51 @@ int otbFunctorImageFilter(int itkNotUsed(argc), char * itkNotUsed(argv) [])
   image->Allocate();
   image->FillBuffer(0);
 
-  auto iterators = MakeIterators(std::make_tuple(image,vimage),image->GetBufferedRegion());
-
-  MoveIterators(iterators);
-
-  Funct4 f{};
-  
-  auto res = CallOperator(f,iterators);
-
   // Test VariadicInputsImageFilter
-
   auto filter = otb::VariadicInputsImageFilter<VectorImageType,VectorImageType,ImageType>::New();
-
   filter->SetVInput<0>(vimage);
   filter->SetVInput<1>(image);
-
   filter->SetVInputs(vimage,image);
-
   std::cout<<filter->GetVInput<0>()<< filter->GetVInput<1>()<<std::endl;
 
   
-  // test FunctorImageFilter
-  auto filter1 = otb::FunctorImageFilter<Funct1>::New(Funct1{});
-  auto filter2 = otb::FunctorImageFilter<Funct2>::New(Funct2{});
-  //auto filter3 = otb::FunctorImageFilter<Funct3>::New(Funct3{});
-  auto filter4 = otb::FunctorImageFilter<Funct4>::New(Funct4{});
+  // test FunctorImageFilter with a lambda
+  double scale = 10.;  
+  auto Lambda1 = [scale](double p)
+               {
+                 return scale*p;
+               };
   auto filterLambda = otb::FunctorImageFilter<decltype(Lambda1)>::New(Lambda1);
+  filterLambda->SetVInputs(image);
+  filterLambda->Update();
 
-  filter1->SetVInputs(image);
-  filter1->Update();
+  // Test FunctorImageFilter with the VariadicConcatenate operator
+  using ConcatFunctorType = VariadicConcatenate<double, double, itk::VariableLengthVector<double> >;
+  auto concatenate = otb::FunctorImageFilter<ConcatFunctorType>::New(ConcatFunctorType{});
+  concatenate->SetVInputs(image,vimage);
+  concatenate->Update();
   
-  filter2->SetVInputs(image);
-  filter2->Update();
+  // Test FunctorImageFilter With VariadicAdd functor
+  using AddFunctorType = VariadicAdd<double, double, double>;
+  auto add = otb::FunctorImageFilter<AddFunctorType>::New(AddFunctorType{});
+  add->SetVInputs(image,image);
+  add->Update();
 
-  filter4->SetVInputs(image,vimage);
-  filter4->Update();
-
-  filterLambda->SetVInputs(image);
-  filterLambda->Update();
+  // Example of simple helper function to make things easier
+  auto add2 = make_add_filter<double,double,double>(image,image);
 
+  // Test FunctorImageFilter with BandExtraction functor
+  using ExtractFunctorType = BandExtraction<double,double>;
+  ExtractFunctorType extractFunctor{1,2};
+  auto extract = otb::FunctorImageFilter<ExtractFunctorType>::New(extractFunctor);
+  extract->SetVInputs(vimage);
+  extract->Update();
+  
+  // Test FunctorImageFilter With Median functor
+  using MedianFunctorType = Median<double,double>;
+  auto median = otb::FunctorImageFilter<MedianFunctorType>::New(MedianFunctorType{});
+  median->SetVInputs(image);
+  median->Update();
   
  return EXIT_SUCCESS;
 }