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

/**
 * @file threshold.hpp
 * @brief 閾値処理 (OCV-10: threshold / adaptiveThreshold)
 *
 * 大域閾値 (固定, Otsu, Triangle) + 適応的閾値 (ガウシアン / 平均)
 */

#ifndef CALX_IMAGE_IMGPROC_THRESHOLD_HPP
#define CALX_IMAGE_IMGPROC_THRESHOLD_HPP

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

namespace calx::image {

enum class ThresholdType { Binary, BinaryInv, Truncate, ToZero, ToZeroInv };
enum class AdaptiveMethod { Mean, Gaussian };

/// 大域閾値処理
template<typename T>
Image<T> threshold(const Image<T>& img, T thresh, T maxVal,
                    ThresholdType type = ThresholdType::Binary) {
    Image<T> result(img.height(), img.width());
    for (std::size_t y = 0; y < img.height(); ++y)
        for (std::size_t x = 0; x < img.width(); ++x) {
            T v = img(y, x);
            switch (type) {
            case ThresholdType::Binary:     result(y,x) = (v > thresh) ? maxVal : T(0); break;
            case ThresholdType::BinaryInv:  result(y,x) = (v > thresh) ? T(0) : maxVal; break;
            case ThresholdType::Truncate:   result(y,x) = (v > thresh) ? thresh : v; break;
            case ThresholdType::ToZero:     result(y,x) = (v > thresh) ? v : T(0); break;
            case ThresholdType::ToZeroInv:  result(y,x) = (v > thresh) ? T(0) : v; break;
            }
        }
    return result;
}

/// Otsu の自動閾値決定 (8bit 画像向け、0-1 float にも対応)
template<typename T>
T otsuThreshold(const Image<T>& img) {
    constexpr int nBins = 256;
    std::vector<double> hist(nBins, 0);
    double total = 0;

    for (std::size_t y = 0; y < img.height(); ++y)
        for (std::size_t x = 0; x < img.width(); ++x) {
            int bin = std::clamp(int(double(img(y,x)) * (nBins - 1)), 0, nBins - 1);
            hist[bin]++;
            total++;
        }

    double sumAll = 0;
    for (int i = 0; i < nBins; ++i) sumAll += i * hist[i];

    double sumB = 0, wB = 0, maxVar = 0;
    int bestT = 0;

    for (int t = 0; t < nBins; ++t) {
        wB += hist[t];
        if (wB == 0) continue;
        double wF = total - wB;
        if (wF == 0) break;
        sumB += t * hist[t];
        double mB = sumB / wB;
        double mF = (sumAll - sumB) / wF;
        double var = wB * wF * (mB - mF) * (mB - mF);
        if (var > maxVar) { maxVar = var; bestT = t; }
    }
    return T(double(bestT) / (nBins - 1));
}

/// 適応的閾値処理
template<typename T>
Image<T> adaptiveThreshold(const Image<T>& img, T maxVal, AdaptiveMethod method,
                            ThresholdType type, int blockSize, T C_offset) {
    int h = static_cast<int>(img.height());
    int w = static_cast<int>(img.width());
    int half = blockSize / 2;

    Image<T> result(h, w);

    // 積分画像 (Mean 用)
    std::vector<std::vector<double>> integral;
    if (method == AdaptiveMethod::Mean) {
        integral.resize(h + 1, std::vector<double>(w + 1, 0));
        for (int y = 0; y < h; ++y)
            for (int x = 0; x < w; ++x)
                integral[y+1][x+1] = double(img(y,x))
                    + integral[y][x+1] + integral[y+1][x] - integral[y][x];
    }

    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x) {
            double localMean;
            if (method == AdaptiveMethod::Mean) {
                int y0 = std::max(0, y - half), y1 = std::min(h - 1, y + half);
                int x0 = std::max(0, x - half), x1 = std::min(w - 1, x + half);
                double area = (y1-y0+1) * (x1-x0+1);
                double sum = integral[y1+1][x1+1] - integral[y0][x1+1]
                           - integral[y1+1][x0] + integral[y0][x0];
                localMean = sum / area;
            } else {
                // Gaussian 重み付き平均
                double sum = 0, wSum = 0;
                double sigma = blockSize / 6.0;
                for (int dy = -half; dy <= half; ++dy)
                    for (int dx = -half; dx <= half; ++dx) {
                        int ny = y + dy, nx = x + dx;
                        if (ny < 0 || ny >= h || nx < 0 || nx >= w) continue;
                        double w = std::exp(-(dx*dx + dy*dy) / (2*sigma*sigma));
                        sum += w * double(img(ny, nx));
                        wSum += w;
                    }
                localMean = sum / wSum;
            }

            T thresh = T(localMean - double(C_offset));
            T v = img(y, x);
            switch (type) {
            case ThresholdType::Binary:    result(y,x) = (v > thresh) ? maxVal : T(0); break;
            case ThresholdType::BinaryInv: result(y,x) = (v > thresh) ? T(0) : maxVal; break;
            default: result(y,x) = (v > thresh) ? maxVal : T(0); break;
            }
        }
    return result;
}

} // namespace calx::image

#endif
