diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
index 6c06c2b347ba955df00afd5b1b54eab50f24db0f..69766a8b319d76704f5839c8e8749346cdca4c79 100644
--- a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
+++ b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.h
@@ -30,6 +30,7 @@
 #include "itkImageRegionConstIterator.h"
 #include "itkProcessObject.h"
 #include <type_traits>
+#include <utility>
 #include "otbFunctionTraits.h"
 #include "itkDefaultConvertPixelTraits.h"
 
@@ -174,7 +175,7 @@ template <class F, class T, size_t N> struct NumberOfOutputComponents<F,otb::Ima
     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;
+// 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)
@@ -276,24 +277,21 @@ template <typename C, typename R, typename... T> struct FunctorFilterSuperclassH
   using InputHasNeighborhood = std::tuple<typename IsNeighborhood<T>::ValueType...>;
 };
 
-// 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;
-// };
+
+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 
@@ -303,7 +301,10 @@ template <typename C, typename R, typename... T> struct FunctorFilterSuperclassH
  * \ingroup IntensityImageFilters   Multithreaded
  *
  * \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
@@ -324,7 +325,7 @@ public:
   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}});
+ // static Pointer New(const TFunction& f, itk::Size<2> radius = {{0,0}});
 
 /** Run-time type information (and related methods). */
   itkTypeMacro(FunctorImageFilter, ImageToImageFilter);
@@ -362,14 +363,14 @@ void SetFunctor(const FunctorType& functor)
 }
 
 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() {}
 
-protected:
-
 /** 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
@@ -394,6 +395,23 @@ protected:
   itk::Size<2> m_Radius;
 };
 
+template <typename Functor> auto NewFunctorFilter(const Functor& f, itk::Size<2> radius)
+{
+  using FilterType = FunctorImageFilter<Functor>;
+  using PointerType = typename FilterType::Pointer;
+
+  PointerType  p = new FilterType(f,radius);
+  return p;
+}
+
+template <typename Functor> auto NewFunctorFilter(const Functor& f, unsigned int numberOfOutputBands, itk::Size<2> radius)
+{
+  using FunctorType = NumberOfOutputBandsDecorator<Functor>;
+  FunctorType decoratedF(f,numberOfOutputBands);
+  return  NewFunctorFilter(decoratedF,radius);
+}
+
+
 }// namespace otb
 
 #ifndef OTB_MANUAL_INSTANTIATION
diff --git a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
index b163697616a893e96d99a4b8c344b1e4b5b1f79f..b22082af9817a4cc5bb54832d2cbdaf7fd92a284 100644
--- a/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
+++ b/Modules/Filtering/ImageManipulation/include/otbFunctorImageFilter.hxx
@@ -32,14 +32,6 @@ template <class TFunction>
 FunctorImageFilter<TFunction>::FunctorImageFilter()
 {}
 
-template <class TFunction>
-typename FunctorImageFilter<TFunction>::Pointer
-FunctorImageFilter<TFunction>::New(const TFunction& f, itk::Size<2> radius) 
-{
-  Pointer p = new FunctorImageFilter<TFunction>(f,radius);
-  return p;                                           
-}
-
 template <class TFunction>
 void
 FunctorImageFilter<TFunction>
diff --git a/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx b/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx
index 6416d222368027dcb2e168eac77e8165775f3b4e..4cb02775d5b147be0da775b3c7b8a877293067ff 100644
--- a/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx
+++ b/Modules/Filtering/ImageManipulation/test/otbFunctorImageFilter.cxx
@@ -42,26 +42,6 @@ template <typename O, typename ...T> struct VariadicAdd
   }
 };
 
-// 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)
@@ -170,15 +150,32 @@ template<typename TOut, typename TIn> struct Mean
   }
 };
 
+// 1 Image with neighborhood of VariableLengthVector -> 1 image with
+// VariableLengthVector
+// For each channel, returns the maximum value in neighborhood
+template<typename T> struct MaxInEachChannel
+{
+  auto operator()(const itk::Neighborhood<itk::VariableLengthVector<T>> & in) const
+  {
+    auto out = in.GetCenterValue();
+
+    for(auto it = in.Begin(); it!=in.End(); ++it)
+      {
+      for(auto band = 0u; band < out.Size();++band)
+        {
+        if((*it)[band]>out[band])
+          out[band] = (*it)[band];
+        }
+      }
+    return out;
+  }
 
+  size_t OutputSize(const std::array<size_t,1> & nbBands) const
+  {
+    return nbBands[0];
+  } 
+};
 
-// // 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>;
@@ -229,37 +226,60 @@ int otbFunctorImageFilter(int itkNotUsed(argc), char * itkNotUsed(argv) [])
                {
                  return scale*p;
                };
-  auto filterLambda = otb::FunctorImageFilter<decltype(Lambda1)>::New(Lambda1);
+  auto filterLambda = NewFunctorFilter(Lambda1);
   filterLambda->SetVInputs(image);
   filterLambda->Update();
 
+  // test FunctorImageFilter with a lambda that returns a
+  // VariableLengthVector
+  // Converts a neighborhood to a VariableLengthVector
+  auto Lambda2 = [](const itk::Neighborhood<double> & in)
+                 {
+                   itk::VariableLengthVector<double> out(in.Size());
+                   std::size_t idx{0};
+                   for(auto it = in.Begin(); it!=in.End();++it,++idx)
+                     {
+                     out[idx]=*it;
+                     }
+                   return out;
+                 };
+
+  // In this case, we use the helper function which allows to specify
+  // the number of outputs
+  auto filterLambda2  = NewFunctorFilter(Lambda2,vimage->GetNumberOfComponentsPerPixel(),{{3,3}});
+  filterLambda2->SetVInputs(image);
+  filterLambda2->Update();
+  
   // Test FunctorImageFilter with the VariadicConcatenate operator
   using ConcatFunctorType = VariadicConcatenate<double, double, itk::VariableLengthVector<double> >;
-  auto concatenate = otb::FunctorImageFilter<ConcatFunctorType>::New(ConcatFunctorType{});
+  auto concatenate = NewFunctorFilter(ConcatFunctorType{});
   concatenate->SetVInputs(image,vimage);
   concatenate->Update();
   
   // Test FunctorImageFilter With VariadicAdd functor
   using AddFunctorType = VariadicAdd<double, double, double>;
-  auto add = otb::FunctorImageFilter<AddFunctorType>::New(AddFunctorType{});
+  auto add = NewFunctorFilter(AddFunctorType{});
   add->SetVInputs(image,image);
   add->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);
+  auto extract = NewFunctorFilter(extractFunctor);
   extract->SetVInputs(vimage);
   extract->Update();
   
   // Test FunctorImageFilter With Mean functor
   using MeanFunctorType = Mean<double,double>;
-  auto median = otb::FunctorImageFilter<MeanFunctorType>::New(MeanFunctorType{},{{2,2}});
+  auto median = NewFunctorFilter(MeanFunctorType{},{{2,2}});
   median->SetVInputs(image);
   median->Update();
+
+  // Test FunctorImageFilter with MaxInEachChannel
+  using MaxInEachChannelType = MaxInEachChannel<double>;
+  auto maxInEachChannel = NewFunctorFilter(MaxInEachChannelType{},{{3,3}});
+  maxInEachChannel->SetVInputs(vimage);
+  maxInEachChannel->Update();
   
  return EXIT_SUCCESS;
 }