add cubic interpolation for Resize

include/aidge/backend/cpu/data/Interpolation.hpp

@@ -63,7 +63,8 @@ class InterpolationCPU : public Interpolation {
     template <typename T>
     static T interpolate(const std::vector<float> &coordsToInterpolate,
                          const std::set<Point<T>> &points,
-                         const Mode interpMode = Interpolation::Mode::Linear);
+                         const Mode interpMode = Interpolation::Mode::Linear,
+                         float cubic_coeff_a=-0.75f);
 
     /**
      * @brief performs linear interpolation on given points.
@@ -75,6 +76,18 @@ class InterpolationCPU : public Interpolation {
     static T linear(const std::vector<float> &originalCoords,
                     const std::set<Point<T>> &points);
 
+  template <typename T>
+  static T cubicInterpolate1D(T p0, T p1, T p2, T p3, float x, float cubic_coeff_a= -0.75f);
+  template <typename T>
+  static T cubicRecurse(const std::vector<float> &coordToInterpolate,
+                                   const std::set<Point<T>> &points,
+                                   DimIdx_t alongDim,
+                                   float cubic_coeff_a = -0.75f);
+
+  template <typename T>
+  static T cubic(const std::vector<float> &coordToInterpolate,
+      const std::set<Point<T>> &pointsToInterpolate,
+      float cubic_coeff_a = -0.75f);
     /**
      * @brief performs nearest interpolation on given points.
      * @note it is a wrapper for linearRecurse() private method

include/aidge/backend/cpu/operator/ResizeImpl_kernels.hpp

@@ -91,6 +91,7 @@ void ResizeImpl_cpu_forward_kernel(
                 InterpolationCPU::retrieveNeighbours(input,
                                                      inputDims,
                                                      coordInApprox,
+                                                     interpMode,
                                                      paddingMode);
             output[idxFlatOut] = InterpolationCPU::interpolate(coordInApprox,
                                                                neighbours,

src/data/Interpolation.cpp

@@ -286,17 +286,104 @@ T InterpolationCPU::nearest(const std::vector<float> &coordsToInterpolate,
     }
 }
 
+template <typename T>
+T cubicWeight(float x, float a) {
+    x = std::fabs(x);
+    if (x < 1.0f) {
+        return static_cast<T>((a + 2) * x * x * x - (a + 3) * x * x + 1);
+    } else if (x < 2.0f) {
+        return static_cast<T>(a * x * x * x - 5 * a * x * x + 8 * a * x - 4 * a);
+    } else {
+        return static_cast<T>(0);
+    }
+}
+
+template <typename T>
+T InterpolationCPU::cubicInterpolate1D(T p0, T p1, T p2, T p3, float x, float a) {
+    // Compute weights based on the Keys cubic kernel
+    float w0 = cubicWeight<T>(1 + x, a);
+    float w1 = cubicWeight<T>(x, a);
+    float w2 = cubicWeight<T>(1 - x, a);
+    float w3 = cubicWeight<T>(2 - x, a);
+
+    return static_cast<T>(p0 * w0 + p1 * w1 + p2 * w2 + p3 * w3);
+}
+
+template <typename T>
+T InterpolationCPU::cubicRecurse(const std::vector<float> &coordToInterpolate,
+                                 const std::set<Point<T>> &points,
+                                 DimIdx_t alongDim,
+                                 float a) {
+    AIDGE_ASSERT(points.size() >= 4,
+                 "Expected at least 4 points at dim {}, got {}", alongDim, points.size());
+
+    if (alongDim == coordToInterpolate.size() - 1) {
+        // Final dimension, extract 4 values in order
+        std::map<int, T> scalarByIndex;
+        for (const auto &pt : points) {
+            scalarByIndex[pt.first[alongDim]] = pt.second;
+        }
+
+        AIDGE_ASSERT(scalarByIndex.size() == 4,
+                     "Expected 4 scalar points along dim {}, got {}", alongDim, scalarByIndex.size());
+
+        auto it = scalarByIndex.begin();
+        float t = coordToInterpolate[alongDim] - static_cast<float>(std::next(it, 1)->first);
+        return cubicInterpolate1D(
+            it->second,
+            std::next(it, 1)->second,
+            std::next(it, 2)->second,
+            std::next(it, 3)->second,
+            t,
+            a);
+    }
+
+    // Non-final dimension: group points by current dimension's coord
+    std::map<int, std::set<Point<T>>> grouped;
+    for (const auto &pt : points) {
+        grouped[pt.first[alongDim]].insert(pt);
+    }
+
+    AIDGE_ASSERT(grouped.size() == 4,
+                 "Expected 4 groups along dim {}, got {}", alongDim, grouped.size());
+
+    std::vector<T> interpolatedValues;
+    for (const auto &[_, subgroup] : grouped) {
+        interpolatedValues.push_back(
+            cubicRecurse(coordToInterpolate, subgroup, alongDim + 1, a)
+        );
+    }
+
+    auto it = grouped.begin();
+    float t = coordToInterpolate[alongDim] - static_cast<float>(std::next(it, 1)->first);
+    return cubicInterpolate1D(
+        interpolatedValues[0],
+        interpolatedValues[1],
+        interpolatedValues[2],
+        interpolatedValues[3],
+        t,
+        a);
+}
+
+
+template <typename T>
+T InterpolationCPU::cubic(const std::vector<float> &coordToInterpolate,
+                          const std::set<Point<T>> &pointsToInterpolate,
+                          float cubic_coeff_a) {
+    return cubicRecurse(coordToInterpolate, pointsToInterpolate, 0, cubic_coeff_a);
+}
+
 template <typename T>
 T InterpolationCPU::interpolate(const std::vector<float> &coordsToInterpolate,
                                 const std::set<Point<T>> &points,
-                                const Mode interpMode) {
+                                const Mode interpMode,
+                                float cubic_coeff_a) {
 
+    // AIDGE_ASSERT(cubic_coeff_a <0, "a must be <0");
     T result{0};
     switch (interpMode) {
     case Interpolation::Mode::Cubic: {
-        AIDGE_THROW_OR_ABORT(
-            std::runtime_error,
-            "Unsupported interpolation mode selected : Cubic.");
+        return cubic(coordsToInterpolate, points, cubic_coeff_a);
         break;
     }
     case Interpolation::Mode::Linear: {
@@ -331,32 +418,39 @@ T InterpolationCPU::interpolate(const std::vector<float> &coordsToInterpolate,
 template int8_t InterpolationCPU::interpolate<int8_t>(
     const std::vector<float> &originalCoords,
     const std::set<Point<int8_t>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 template int16_t InterpolationCPU::interpolate<int16_t>(
     const std::vector<float> &originalCoords,
     const std::set<Point<int16_t>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 template int32_t InterpolationCPU::interpolate<int32_t>(
     const std::vector<float> &originalCoords,
     const std::set<Point<int32_t>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 template int64_t InterpolationCPU::interpolate<int64_t>(
     const std::vector<float> &originalCoords,
     const std::set<Point<int64_t>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 
 template half_float::half InterpolationCPU::interpolate<half_float::half>(
     const std::vector<float> &originalCoords,
     const std::set<Point<half_float::half>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 template float InterpolationCPU::interpolate<float>(
     const std::vector<float> &originalCoords,
     const std::set<Point<float>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 template double InterpolationCPU::interpolate<double>(
     const std::vector<float> &originalCoords,
     const std::set<Point<double>> &points,
-    const Mode interpMode);
+    const Mode interpMode,
+    float cubic_coeff_a);
 
 ////////////////////////////////////////////////////////////////////
 // INTERPOLATE LINEAR (& its associated recursive function)
@@ -432,5 +526,36 @@ template double
 InterpolationCPU::nearest(const std::vector<float> &originalCoords,
                           const std::set<Point<double>> &points,
                           const Interpolation::Mode nearestMode);
+////////////////////////////////////////////////////////////////////
+// INTERPOLATE CUBIC
+
+template int8_t InterpolationCPU::cubic<int8_t>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<int8_t>> &points,
+    float cubic_coeff_a);
+template int16_t InterpolationCPU::cubic<int16_t>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<int16_t>> &points,
+    float cubic_coeff_a);
+template int32_t InterpolationCPU::cubic<int32_t>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<int32_t>> &points,
+    float cubic_coeff_a);
+template int64_t InterpolationCPU::cubic<int64_t>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<int64_t>> &points,
+    float cubic_coeff_a);
+template half_float::half InterpolationCPU::cubic<half_float::half>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<half_float::half>> &points,
+    float cubic_coeff_a);
+template float InterpolationCPU::cubic<float>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<float>> &points,
+    float cubic_coeff_a);
+template double InterpolationCPU::cubic<double>(
+    const std::vector<float> &coordsToInterpolate,
+    const std::set<Point<double>> &points,
+    float cubic_coeff_a);
 
 } // namespace Aidge

unit_tests/operator/Test_ResizeImpl.cpp

@@ -244,6 +244,41 @@ TEST_CASE("[cpu/operator] Resize(forward)", "[Resize][CPU]") {
                               1e-5f,
                               1e-5f) == true);
     }
+    SECTION("Cubic Interpolation") {
+        std::shared_ptr<Tensor> input_tensor = std::make_shared<Tensor>(Array4D<float, 1, 1, 4, 4>{{{{
+            {1.0f, 2.0f, 3.0f, 4.0f},
+            {5.0f, 6.0f, 7.0f, 8.0f},
+            {9.0f, 10.0f, 11.0f, 12.0f},
+            {13.0f, 14.0f, 15.0f, 16.0f}}}
+        }});
+        Tensor expected_out_tensor = Tensor(Array4D<float, 1, 1, 8, 8>{{{{
+            { 0.47265625, 0.76953125, 1.24609375, 1.87500000, 2.28125000, 2.91015625, 3.38671875, 3.68359375 },
+            { 1.66015625, 1.95703125, 2.43359375, 3.06250000, 3.46875000, 4.09765625, 4.57421875, 4.87109375 },
+            { 3.56640625, 3.86328125, 4.33984375, 4.96875000, 5.37500000, 6.00390625, 6.48046875, 6.77734375 },
+            { 6.08203125, 6.37890625, 6.85546875, 7.48437500, 7.89062500, 8.51953125, 8.99609375, 9.29296875 },
+            { 7.70703125, 8.00390625, 8.48046875, 9.10937500, 9.51562500, 10.14453125, 10.62109375, 10.91796875 },
+            { 10.22265625, 10.51953125, 10.99609375, 11.62500000, 12.03125000, 12.66015625, 13.13671875, 13.43359375 },
+            { 12.12890625, 12.42578125, 12.90234375, 13.53125000, 13.93750000, 14.56640625, 15.04296875, 15.33984375 },
+            { 13.31640625, 13.61328125, 14.08984375, 14.71875000, 15.12500000, 15.75390625, 16.23046875, 16.52734375}
+        }}}});      
+
+        std::vector<float> scales = {1.0f, 1.0f, 2.0f, 2.0f}; // Adjust according to the expected interpolation scale
+        auto resize_node = Resize(scales, {}, Interpolation::CoordinateTransformation::HalfPixel, Interpolation::Mode::Cubic);
+        auto op = std::static_pointer_cast<Resize_Op>(resize_node->getOperator());
+        op->associateInput(0, input_tensor);
+    
+        op->setDataType(DataType::Float32);
+        op->setBackend("cpu");
+        op->forwardDims(true);
+        op->forward();
+    
+        op->getOutput(0)->print();
+        expected_out_tensor.print();
+    
+        CHECK(approxEq<float>(*op->getOutput(0), expected_out_tensor, 1e-5f, 1e-5f) == true);
+    }
+    
+    
 }
 
 } // namespace Aidge