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

// canny.hpp — Canny エッジ検出
//
// 1. ガウシアンぼかし
// 2. Sobel 勾配 (方向 + 強度)
// 3. 非最大値抑制 (NMS)
// 4. ヒステリシス閾値処理

#ifndef CALX_IMAGE_IMGPROC_CANNY_HPP
#define CALX_IMAGE_IMGPROC_CANNY_HPP

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

namespace calx::image {

    /// Canny エッジ検出
    /// @param src         入力グレースケール画像 (1ch)
    /// @param lowThresh   低閾値
    /// @param highThresh  高閾値
    /// @param sigma       ガウシアンぼかしの σ (0 = ぼかしなし)
    /// @return エッジ画像 (1ch, 0 or 1)
    template<typename T>
    Image<T> cannyEdge(const Image<T>& src, T lowThresh, T highThresh,
                        T sigma = T{1.4})
    {
        assert(src.numChannels() == 1);
        size_t H = src.height(), W = src.width();

        // Step 1: ガウシアンぼかし
        std::vector<std::vector<double>> blurred(H, std::vector<double>(W, 0));
        if (sigma > T{0.01}) {
            int ksize = static_cast<int>(std::ceil(sigma * 3)) * 2 + 1;
            int half = ksize / 2;
            std::vector<double> kernel(ksize);
            double sum = 0;
            for (int i = 0; i < ksize; ++i) {
                double x = i - half;
                kernel[i] = std::exp(-x * x / (2 * static_cast<double>(sigma * sigma)));
                sum += kernel[i];
            }
            for (auto& k : kernel) k /= sum;

            // 水平 pass
            std::vector<std::vector<double>> tmp(H, std::vector<double>(W, 0));
            for (size_t y = 0; y < H; ++y)
                for (size_t x = 0; x < W; ++x) {
                    double v = 0;
                    for (int k = -half; k <= half; ++k) {
                        int xx = std::clamp(static_cast<int>(x) + k, 0, static_cast<int>(W) - 1);
                        v += kernel[k + half] * static_cast<double>(src.at(y, xx, 0));
                    }
                    tmp[y][x] = v;
                }
            // 垂直 pass
            for (size_t y = 0; y < H; ++y)
                for (size_t x = 0; x < W; ++x) {
                    double v = 0;
                    for (int k = -half; k <= half; ++k) {
                        int yy = std::clamp(static_cast<int>(y) + k, 0, static_cast<int>(H) - 1);
                        v += kernel[k + half] * tmp[yy][x];
                    }
                    blurred[y][x] = v;
                }
        } else {
            for (size_t y = 0; y < H; ++y)
                for (size_t x = 0; x < W; ++x)
                    blurred[y][x] = static_cast<double>(src.at(y, x, 0));
        }

        // Step 2: Sobel 勾配
        std::vector<std::vector<double>> Gx(H, std::vector<double>(W, 0));
        std::vector<std::vector<double>> Gy(H, std::vector<double>(W, 0));
        std::vector<std::vector<double>> mag(H, std::vector<double>(W, 0));
        std::vector<std::vector<double>> dir(H, std::vector<double>(W, 0));

        for (size_t y = 1; y + 1 < H; ++y)
            for (size_t x = 1; x + 1 < W; ++x) {
                Gx[y][x] = -blurred[y-1][x-1] + blurred[y-1][x+1]
                           -2*blurred[y][x-1] + 2*blurred[y][x+1]
                           -blurred[y+1][x-1] + blurred[y+1][x+1];
                Gy[y][x] = -blurred[y-1][x-1] - 2*blurred[y-1][x]
                           -blurred[y-1][x+1]
                           +blurred[y+1][x-1] + 2*blurred[y+1][x]
                           +blurred[y+1][x+1];
                mag[y][x] = std::sqrt(Gx[y][x]*Gx[y][x] + Gy[y][x]*Gy[y][x]);
                dir[y][x] = std::atan2(Gy[y][x], Gx[y][x]);
            }

        // Step 3: 非最大値抑制
        std::vector<std::vector<double>> nms(H, std::vector<double>(W, 0));
        for (size_t y = 1; y + 1 < H; ++y)
            for (size_t x = 1; x + 1 < W; ++x) {
                double angle = dir[y][x];
                if (angle < 0) angle += 3.14159265358979;
                double q = 0, r = 0;

                // 0°
                if (angle < 0.3927 || angle >= 2.7489) {
                    q = mag[y][x+1]; r = mag[y][x-1];
                }
                // 45°
                else if (angle < 1.1781) {
                    q = mag[y-1][x+1]; r = mag[y+1][x-1];
                }
                // 90°
                else if (angle < 1.9635) {
                    q = mag[y-1][x]; r = mag[y+1][x];
                }
                // 135°
                else {
                    q = mag[y-1][x-1]; r = mag[y+1][x+1];
                }

                if (mag[y][x] >= q && mag[y][x] >= r)
                    nms[y][x] = mag[y][x];
            }

        // Step 4: ヒステリシス閾値 (BFS)
        Image<T> result(H, W, 1);
        double lo = static_cast<double>(lowThresh);
        double hi = static_cast<double>(highThresh);

        // strong edges
        std::vector<std::vector<int>> label(H, std::vector<int>(W, 0));
        // 0=suppressed, 1=weak, 2=strong
        std::queue<std::pair<size_t, size_t>> q;
        for (size_t y = 0; y < H; ++y)
            for (size_t x = 0; x < W; ++x) {
                if (nms[y][x] >= hi) {
                    label[y][x] = 2;
                    q.push({y, x});
                } else if (nms[y][x] >= lo) {
                    label[y][x] = 1;
                }
            }

        // BFS: strong に隣接する weak を strong に昇格
        int dx8[] = {-1,0,1,-1,1,-1,0,1};
        int dy8[] = {-1,-1,-1,0,0,1,1,1};
        while (!q.empty()) {
            auto [cy, cx] = q.front(); q.pop();
            for (int d = 0; d < 8; ++d) {
                int ny = static_cast<int>(cy) + dy8[d];
                int nx = static_cast<int>(cx) + dx8[d];
                if (ny >= 0 && ny < static_cast<int>(H) &&
                    nx >= 0 && nx < static_cast<int>(W) &&
                    label[ny][nx] == 1) {
                    label[ny][nx] = 2;
                    q.push({static_cast<size_t>(ny), static_cast<size_t>(nx)});
                }
            }
        }

        for (size_t y = 0; y < H; ++y)
            for (size_t x = 0; x < W; ++x)
                result.at(y, x, 0) = (label[y][x] == 2) ? T{1} : T{0};

        return result;
    }

} // namespace calx::image

#endif // CALX_IMAGE_IMGPROC_CANNY_HPP
