diff --git a/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.h b/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.h
index b1a925b458c77fe571a20e03449239fd52fb3aa8..f70d93c6b6d7fc3222c125e4c2dd166ebc8a3817 100644
--- a/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.h
+++ b/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.h
@@ -21,8 +21,14 @@
 #ifndef otbBCOInterpolateImageFunction_h
 #define otbBCOInterpolateImageFunction_h
 
-#include "itkInterpolateImageFunction.h"
+#include <boost/version.hpp>
+#if BOOST_VERSION >= 105800
+#include <boost/container/small_vector.hpp>
+#else
 #include "vnl/vnl_vector.h"
+#endif
+
+#include "itkInterpolateImageFunction.h"
 #include "otbMath.h"
 
 #include "otbVectorImage.h"
@@ -57,7 +63,7 @@ class ITK_EXPORT BCOInterpolateImageFunctionBase : public itk::InterpolateImageF
 {
 public:
   /** Standard class typedefs. */
-  typedef BCOInterpolateImageFunctionBase Self;
+  typedef BCOInterpolateImageFunctionBase                       Self;
   typedef itk::InterpolateImageFunction<TInputImage, TCoordRep> Superclass;
 
   /** Run-time type information (and related methods). */
@@ -89,15 +95,21 @@ public:
   typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
   typedef TCoordRep                                ContinuousIndexValueType;
 
-  /** Coeficients container type.*/
+#if BOOST_VERSION >= 105800
+  // Faster path for small radii.
+  /** Coeficients container type. */
+  typedef boost::container::small_vector<double, 7> CoefContainerType;
+#else
+  /** Coeficients container type. */
   typedef vnl_vector<double> CoefContainerType;
+#endif
 
   /** Set/Get the window radius */
-  virtual void SetRadius(unsigned int radius);
+  virtual void         SetRadius(unsigned int radius);
   virtual unsigned int GetRadius() const;
 
   /** Set/Get the optimisation coefficient (Common values are -0.5, -0.75 or -1.0) */
-  virtual void SetAlpha(double alpha);
+  virtual void   SetAlpha(double alpha);
   virtual double GetAlpha() const;
 
   /** Evaluate the function at a ContinuousIndex position
@@ -115,7 +127,7 @@ protected:
   ~BCOInterpolateImageFunctionBase() override{};
   void PrintSelf(std::ostream& os, itk::Indent indent) const override;
   /** Compute the BCO coefficients. */
-  virtual CoefContainerType EvaluateCoef(const ContinuousIndexValueType& indexValue) const;
+  CoefContainerType EvaluateCoef(const ContinuousIndexValueType& indexValue) const;
 
   /** Used radius for the BCO */
   unsigned int m_Radius;
@@ -135,10 +147,10 @@ class ITK_EXPORT BCOInterpolateImageFunction : public otb::BCOInterpolateImageFu
 {
 public:
   /** Standard class typedefs. */
-  typedef BCOInterpolateImageFunction Self;
+  typedef BCOInterpolateImageFunction                             Self;
   typedef BCOInterpolateImageFunctionBase<TInputImage, TCoordRep> Superclass;
-  typedef itk::SmartPointer<Self>       Pointer;
-  typedef itk::SmartPointer<const Self> ConstPointer;
+  typedef itk::SmartPointer<Self>                                 Pointer;
+  typedef itk::SmartPointer<const Self>                           ConstPointer;
 
   itkTypeMacro(BCOInterpolateImageFunction, BCOInterpolateImageFunctionBase);
   itkNewMacro(Self);
@@ -169,14 +181,14 @@ private:
 
 template <typename TPixel, unsigned int VImageDimension, class TCoordRep>
 class ITK_EXPORT BCOInterpolateImageFunction<otb::VectorImage<TPixel, VImageDimension>, TCoordRep>
-    : public otb::BCOInterpolateImageFunctionBase<otb::VectorImage<TPixel, VImageDimension>, TCoordRep>
+  : public otb::BCOInterpolateImageFunctionBase<otb::VectorImage<TPixel, VImageDimension>, TCoordRep>
 {
 public:
   /** Standard class typedefs.*/
-  typedef BCOInterpolateImageFunction Self;
+  typedef BCOInterpolateImageFunction                                                           Self;
   typedef BCOInterpolateImageFunctionBase<otb::VectorImage<TPixel, VImageDimension>, TCoordRep> Superclass;
-  typedef itk::SmartPointer<Self>       Pointer;
-  typedef itk::SmartPointer<const Self> ConstPointer;
+  typedef itk::SmartPointer<Self>                                                               Pointer;
+  typedef itk::SmartPointer<const Self>                                                         ConstPointer;
 
   itkTypeMacro(BCOInterpolateImageFunction, BCOInterpolateImageFunctionBase);
   itkNewMacro(Self);
diff --git a/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.hxx b/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.hxx
index a1c27a37b0cacc207b122a46f20aad4716ef1570..9070621aee0aaf1c59d55e896d0d98b51ddc1fea 100644
--- a/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.hxx
+++ b/Modules/Core/Interpolation/include/otbBCOInterpolateImageFunction.hxx
@@ -72,10 +72,9 @@ template <class TInputImage, class TCoordRep>
 typename BCOInterpolateImageFunctionBase<TInputImage, TCoordRep>::CoefContainerType
 BCOInterpolateImageFunctionBase<TInputImage, TCoordRep>::EvaluateCoef(const ContinuousIndexValueType& indexValue) const
 {
-  // Init BCO coefficient container
+  double offset, dist, position, step;
 
-  CoefContainerType BCOCoef(m_WinSize, 0.);
-  double            offset, dist, position, step;
+  typename BCOInterpolateImageFunctionBase<TInputImage, TCoordRep>::CoefContainerType bcoCoef(this->m_WinSize);
 
   offset = indexValue - itk::Math::Floor<IndexValueType>(indexValue + 0.5);
 
@@ -83,38 +82,41 @@ BCOInterpolateImageFunctionBase<TInputImage, TCoordRep>::EvaluateCoef(const Cont
   step     = 4. / static_cast<double>(2 * m_Radius);
   position = -static_cast<double>(m_Radius) * step;
 
-  double sum = 0.0;
+  double sum   = 0.0;
+  auto   alpha = this->m_Alpha;
 
   for (unsigned int i = 0; i < m_WinSize; ++i)
   {
-
     // Compute the BCO coefficients according to alpha.
     dist = std::abs(position - offset * step);
 
+    double coef;
     if (dist <= 2.)
     {
       if (dist <= 1.)
       {
-        BCOCoef[i] = (m_Alpha + 2.) * std::abs(dist * dist * dist) - (m_Alpha + 3.) * dist * dist + 1;
+        coef = (alpha + 2.) * dist * dist * dist - (alpha + 3.) * dist * dist + 1;
       }
       else
       {
-        BCOCoef[i] = m_Alpha * std::abs(dist * dist * dist) - 5 * m_Alpha * dist * dist + 8 * m_Alpha * std::abs(dist) - 4 * m_Alpha;
+        coef = alpha * dist * dist * dist - 5 * alpha * dist * dist + 8 * alpha * dist - 4 * alpha;
       }
     }
     else
     {
-      BCOCoef[i] = 0;
+      coef = 0.0;
     }
 
-    sum += BCOCoef[i];
+    sum += coef;
     position += step;
+
+    bcoCoef[i] = coef;
   }
 
   for (unsigned int i = 0; i < m_WinSize; ++i)
-    BCOCoef[i]        = BCOCoef[i] / sum;
+    bcoCoef[i] /= sum;
 
-  return BCOCoef;
+  return bcoCoef;
 }
 
 template <class TInputImage, class TCoordRep>
@@ -127,7 +129,6 @@ template <class TInputImage, class TCoordRep>
 typename BCOInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
 BCOInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateAtContinuousIndex(const ContinuousIndexType& index) const
 {
-
   unsigned int dim;
 
   IndexType baseIndex;
@@ -135,8 +136,8 @@ BCOInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateAtContinuousIndex(c
 
   RealType value = itk::NumericTraits<RealType>::Zero;
 
-  CoefContainerType BCOCoefX = this->EvaluateCoef(index[0]);
-  CoefContainerType BCOCoefY = this->EvaluateCoef(index[1]);
+  const auto& BCOCoefX = this->EvaluateCoef(index[0]);
+  const auto& BCOCoefY = this->EvaluateCoef(index[1]);
 
   // Compute base index = closet index
   for (dim = 0; dim < ImageDimension; dim++)
@@ -198,12 +199,18 @@ BCOInterpolateImageFunction<otb::VectorImage<TPixel, VImageDimension>, TCoordRep
   IndexType neighIndex;
 
 
+#if BOOST_VERSION >= 105800
+  // faster path for <= 8 components
+  boost::container::small_vector<ScalarRealType, 8> lineRes(componentNumber);
+#else
   std::vector<ScalarRealType> lineRes(componentNumber);
-  OutputType                  output(componentNumber);
+#endif
+
+  OutputType output(componentNumber);
   output.Fill(itk::NumericTraits<ScalarRealType>::Zero);
 
-  CoefContainerType BCOCoefX = this->EvaluateCoef(index[0]);
-  CoefContainerType BCOCoefY = this->EvaluateCoef(index[1]);
+  const auto& BCOCoefX = this->EvaluateCoef(index[0]);
+  const auto& BCOCoefY = this->EvaluateCoef(index[1]);
 
   // Compute base index = closet index
   for (dim = 0; dim < ImageDimension; dim++)