WIP: A really, really generic functor filter

double scale = 10.;
// dummy lambda that multiply by a scalar
auto lambda = [scale](double p)
{
  return scale*p;
};
// Create a filter that will accept an otb::Image<double> as input and produce an otb::Image<double>
auto filterLambda = NewFunctorFilter(lambda);
// Set the inputs using variadic SetVInputs
filterLambda->SetVInputs(image);
// All good, lets run it
filterLambda->Update();

// dummy lambda that sums 3 images with different types
auto lambda = [](double p, unsigned int p1, short p2)
{
  return static_cast<double>(p)+static_cast<double>(p1) + static_cast<double>(p2);
};
// Create a filter that will accept otb::Image<double>, otb::Image<unsigned int>,otb::Image<short>  as inputs and produce an otb::Image<double>
auto filterLambda = NewFunctorFilter(lambda);
// Set the inputs using variadic SetVinputs
filterLambda->SetVInputs(image, image0, image1); // image is otb::Image<double>, image0 is otb::Image<unsigned int> and image1 is otb::Image<short>
// All good, lets run it
filterLambda->Update();

// lambda that returns the norm of the input vector pixel
auto lambda = [](const itk::VariableLengthVector<double> & inv)
{
  return inv.GetNorm();
};
// Create a filter that will accept otb::VectorImage<double>, as inputs and produce an otb::Image<double>
auto filterLambda = NewFunctorFilter(lambda);
// Set the inputs using variadic SetVinputs
filterLambda->SetVInputs(vimage); // vimage is otb::VectorImage<double>
// All good, lets run it
filterLambda->Update();

// lambda that performs neighborhood averaging
auto lambda = [](const itk::Neighborhood<short> & in)
{
   double out(0);
   for(auto it = in.Begin(); it!=in.End();++it)
      out+=static_cast<TOut>(*it);

   out/=in.Size(); 
   return out;
};
// Create a filter that will accept otb::Image<short>, as inputs and produce an otb::Image<double>
auto filterLambda = NewFunctorFilter(lambda,{{3,3}});
// Set the inputs using variadic SetVinputs
filterLambda->SetVInputs(image); // image is otb::VectorImage<short>
// All good, lets run it
filterLambda->Update();

template <typename O, typename ...T> struct VariadicAdd
{
  auto operator()(T... ins) const
  {
    std::vector<O> outVector{static_cast<O>(ins)...};
    return std::accumulate(outVector.begin(), outVector.end(),0);
  }
};

using AddFunctorType = VariadicAdd<double, double, double>;
auto add = NewFunctorFilter(AddFunctorType{});
add->SetVInputs(image,image);
add->Update();

template<typename O, typename T> struct BandExtraction
{
  BandExtraction(unsigned int indices...) : m_Indices({indices}){}

  auto operator()(const itk::VariableLengthVector<T> & in) const
  {
    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;
  }

  // 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;
};

using ExtractFunctorType = BandExtraction<double,double>;
ExtractFunctorType extractFunctor{1,2};
auto extract = NewFunctorFilter(extractFunctor);
extract->SetVInputs(vimage);
extract->Update();

// 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)
{
  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<O> out;
    concatenateVectors(out, toVector(ins)...);
    
    return out;
  }

  // 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 ConcatFunctorType = VariadicConcatenate<double, double, itk::VariableLengthVector<double> >;
auto concatenate = NewFunctorFilter(ConcatFunctorType{});
concatenate->SetVInputs(image,vimage);
concatenate->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();