diff --git a/include/aidge/backend/cpu/data/Interpolation.hpp b/include/aidge/backend/cpu/data/Interpolation.hpp index 5909f02a190f4e10cdeb878505fdfea1a17e2d75..bfc32f1f5eccf2696ca64ca6b8d2c76468686d41 100644 --- a/include/aidge/backend/cpu/data/Interpolation.hpp +++ b/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 diff --git a/include/aidge/backend/cpu/operator/ResizeImpl_kernels.hpp b/include/aidge/backend/cpu/operator/ResizeImpl_kernels.hpp index 6449417baf855620669aba11ebca16d9384c4e7c..3adede478d39272470af37879d415ba80de9a0b0 100644 --- a/include/aidge/backend/cpu/operator/ResizeImpl_kernels.hpp +++ b/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, diff --git a/src/data/Interpolation.cpp b/src/data/Interpolation.cpp index fbf224d84f65c442e98967783d303605a177d390..bd1a90249b2c7ad2e385ab2313c9fe25098d7ba9 100644 --- a/src/data/Interpolation.cpp +++ b/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 diff --git a/unit_tests/operator/Test_ResizeImpl.cpp b/unit_tests/operator/Test_ResizeImpl.cpp index 6b3520fc88d36660ff44403bd41a47cd7ed96256..c56e72b6c3f1904311bbbb47c59d176463458021 100644 --- a/unit_tests/operator/Test_ResizeImpl.cpp +++ b/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