// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// resize.hpp — 画像リサイズ (複数補間方式)
//
// nearest-neighbor, bilinear (既存 Image::resize と同等),
// bicubic, Lanczos3 補間。

#ifndef CALX_IMAGE_IMGPROC_RESIZE_HPP
#define CALX_IMAGE_IMGPROC_RESIZE_HPP

#include <image/Image.hpp>
#include <algorithm>
#include <cmath>

namespace calx::image {

    enum class Interpolation { Nearest, Bilinear, Bicubic, Lanczos };

    namespace detail_resize {
        // Bicubic 重み関数 (Catmull-Rom, a=-0.5)
        template<typename T>
        T cubicWeight(T x) {
            x = std::abs(x);
            if (x <= T{1}) return (T{1.5}*x - T{2.5})*x*x + T{1};
            if (x <  T{2}) return ((-T{0.5}*x + T{2.5})*x - T{4})*x + T{2};
            return T{0};
        }

        // Lanczos 重み関数
        template<typename T>
        T lanczosWeight(T x, int a = 3) {
            x = std::abs(x);
            if (x < T{1e-10}) return T{1};
            if (x >= static_cast<T>(a)) return T{0};
            T px = T{3.14159265358979323846} * x;
            T pxa = px / static_cast<T>(a);
            return std::sin(px) / px * std::sin(pxa) / pxa;
        }
    }

    /// 画像リサイズ
    template<typename T>
    Image<T> resize(const Image<T>& src, size_t newH, size_t newW,
                     Interpolation interp = Interpolation::Bilinear)
    {
        size_t srcH = src.height(), srcW = src.width();
        int nch = src.numChannels();
        Image<T> dst(newH, newW, nch);

        double sy = static_cast<double>(srcH) / newH;
        double sx = static_cast<double>(srcW) / newW;

        for (int ch = 0; ch < nch; ++ch) {
            for (size_t y = 0; y < newH; ++y) {
                for (size_t x = 0; x < newW; ++x) {
                    double fy = y * sy, fx = x * sx;

                    switch (interp) {
                    case Interpolation::Nearest: {
                        int iy = static_cast<int>(fy + 0.5);
                        int ix = static_cast<int>(fx + 0.5);
                        iy = std::clamp(iy, 0, static_cast<int>(srcH) - 1);
                        ix = std::clamp(ix, 0, static_cast<int>(srcW) - 1);
                        dst.at(y, x, ch) = src.at(iy, ix, ch);
                        break;
                    }
                    case Interpolation::Bilinear: {
                        int y0 = static_cast<int>(fy);
                        int x0 = static_cast<int>(fx);
                        int y1 = std::min(y0 + 1, static_cast<int>(srcH) - 1);
                        int x1 = std::min(x0 + 1, static_cast<int>(srcW) - 1);
                        double dy = fy - y0, dx = fx - x0;
                        dst.at(y, x, ch) = static_cast<T>(
                            (1-dy)*(1-dx)*src.at(y0,x0,ch) + (1-dy)*dx*src.at(y0,x1,ch) +
                            dy*(1-dx)*src.at(y1,x0,ch) + dy*dx*src.at(y1,x1,ch));
                        break;
                    }
                    case Interpolation::Bicubic: {
                        int iy = static_cast<int>(fy);
                        int ix = static_cast<int>(fx);
                        double dy = fy - iy, dx = fx - ix;
                        double val = 0;
                        for (int m = -1; m <= 2; ++m) {
                            int py = std::clamp(iy + m, 0, static_cast<int>(srcH) - 1);
                            double wy = detail_resize::cubicWeight<double>(dy - m);
                            for (int n = -1; n <= 2; ++n) {
                                int px = std::clamp(ix + n, 0, static_cast<int>(srcW) - 1);
                                double wx = detail_resize::cubicWeight<double>(dx - n);
                                val += wy * wx * static_cast<double>(src.at(py, px, ch));
                            }
                        }
                        dst.at(y, x, ch) = static_cast<T>(val);
                        break;
                    }
                    case Interpolation::Lanczos: {
                        int iy = static_cast<int>(fy);
                        int ix = static_cast<int>(fx);
                        double dy = fy - iy, dx = fx - ix;
                        double val = 0, wsum = 0;
                        for (int m = -2; m <= 3; ++m) {
                            int py = std::clamp(iy + m, 0, static_cast<int>(srcH) - 1);
                            double wy = detail_resize::lanczosWeight<double>(dy - m);
                            for (int n = -2; n <= 3; ++n) {
                                int px = std::clamp(ix + n, 0, static_cast<int>(srcW) - 1);
                                double wx = detail_resize::lanczosWeight<double>(dx - n);
                                double w = wy * wx;
                                val += w * static_cast<double>(src.at(py, px, ch));
                                wsum += w;
                            }
                        }
                        if (wsum > 1e-10) val /= wsum;
                        dst.at(y, x, ch) = static_cast<T>(val);
                        break;
                    }
                    }
                }
            }
        }
        return dst;
    }

} // namespace calx::image

#endif // CALX_IMAGE_IMGPROC_RESIZE_HPP
