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

// bilateral.hpp — Bilateral フィルタ
//
// エッジ保存平滑化フィルタ。空間的近さと輝度値の類似性の両方で
// 重み付けを行い、エッジをぼかさずにノイズを除去する。

#ifndef CALX_IMAGE_FILTERS_BILATERAL_HPP
#define CALX_IMAGE_FILTERS_BILATERAL_HPP

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

namespace calx::image {

    /// Bilateral フィルタ
    /// @param src          入力画像 (任意チャンネル数)
    /// @param d            カーネル直径 (奇数推奨, 0 の場合 sigmaSpace から自動計算)
    /// @param sigmaColor   色空間のガウシアン σ (値域は [0,1] 正規化前提)
    /// @param sigmaSpace   座標空間のガウシアン σ (ピクセル単位)
    /// @return             フィルタ適用後の画像
    template<typename T>
    Image<T> bilateralFilter(const Image<T>& src, int d = 9,
                             T sigmaColor = T{75.0 / 255},
                             T sigmaSpace = T{75.0}) {
        int H = static_cast<int>(src.height());
        int W = static_cast<int>(src.width());
        int nch = src.numChannels();
        if (H == 0 || W == 0) return src;

        // カーネル半径
        if (d <= 0) d = static_cast<int>(std::ceil(sigmaSpace * 3)) * 2 + 1;
        int radius = d / 2;

        double invSigmaColor2 = -0.5 / (static_cast<double>(sigmaColor) * static_cast<double>(sigmaColor));
        double invSigmaSpace2 = -0.5 / (static_cast<double>(sigmaSpace) * static_cast<double>(sigmaSpace));

        // 空間ガウシアンのルックアップテーブル
        std::vector<std::vector<double>> spatialLUT(2 * radius + 1, std::vector<double>(2 * radius + 1));
        for (int dy = -radius; dy <= radius; ++dy)
            for (int dx = -radius; dx <= radius; ++dx)
                spatialLUT[dy + radius][dx + radius] =
                    std::exp((dy * dy + dx * dx) * invSigmaSpace2);

        Image<T> dst(H, W, nch, T{0});

        // 各チャンネルを独立にフィルタリング
        for (int c = 0; c < nch; ++c) {
            const auto& srcCh = src.channel(c);
            auto& dstCh = dst.channel(c);

            for (int y = 0; y < H; ++y) {
                for (int x = 0; x < W; ++x) {
                    double centerVal = static_cast<double>(srcCh(y, x));
                    double sumW = 0;
                    double sumV = 0;

                    for (int dy = -radius; dy <= radius; ++dy) {
                        int ny = std::clamp(y + dy, 0, H - 1);
                        for (int dx = -radius; dx <= radius; ++dx) {
                            int nx = std::clamp(x + dx, 0, W - 1);
                            double neighborVal = static_cast<double>(srcCh(ny, nx));

                            double colorDiff = neighborVal - centerVal;
                            double wColor = std::exp(colorDiff * colorDiff * invSigmaColor2);
                            double wSpace = spatialLUT[dy + radius][dx + radius];
                            double w = wColor * wSpace;

                            sumW += w;
                            sumV += w * neighborVal;
                        }
                    }

                    dstCh(y, x) = static_cast<T>(sumV / sumW);
                }
            }
        }
        return dst;
    }

} // namespace calx::image

#endif // CALX_IMAGE_FILTERS_BILATERAL_HPP
