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

// template_matching.hpp — テンプレートマッチング
//
// SQDIFF:  差分二乗和 (最小が最良一致)
// CCORR:   相互相関
// CCOEFF:  相関係数 (正規化なし)

#ifndef CALX_IMAGE_IMGPROC_TEMPLATE_MATCHING_HPP
#define CALX_IMAGE_IMGPROC_TEMPLATE_MATCHING_HPP

#include <image/Image.hpp>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>

namespace calx::image {

    enum class MatchMethod { SQDIFF, CCORR, CCOEFF };

    /// テンプレートマッチング
    /// @param img     入力画像 (1ch)
    /// @param templ   テンプレート画像 (1ch)
    /// @param method  マッチング手法
    /// @return        結果マップ (H - tH + 1) x (W - tW + 1)
    template<typename T>
    Matrix<T> matchTemplate(const Image<T>& img, const Image<T>& templ,
                            MatchMethod method = MatchMethod::CCOEFF) {
        assert(img.numChannels() >= 1 && templ.numChannels() >= 1);

        size_t H = img.height(), W = img.width();
        size_t tH = templ.height(), tW = templ.width();
        assert(H >= tH && W >= tW);

        size_t rH = H - tH + 1;
        size_t rW = W - tW + 1;
        Matrix<T> result(rH, rW, T{0});

        const auto& imgCh = img.channel(0);
        const auto& tplCh = templ.channel(0);

        // テンプレートの平均値 (CCOEFF で使用)
        T tplMean = T{0};
        if (method == MatchMethod::CCOEFF) {
            for (size_t ty = 0; ty < tH; ++ty)
                for (size_t tx = 0; tx < tW; ++tx)
                    tplMean += tplCh(ty, tx);
            tplMean /= static_cast<T>(tH * tW);
        }

        for (size_t y = 0; y < rH; ++y) {
            for (size_t x = 0; x < rW; ++x) {
                T val = T{0};

                if (method == MatchMethod::SQDIFF) {
                    for (size_t ty = 0; ty < tH; ++ty)
                        for (size_t tx = 0; tx < tW; ++tx) {
                            T diff = imgCh(y + ty, x + tx) - tplCh(ty, tx);
                            val += diff * diff;
                        }
                } else if (method == MatchMethod::CCORR) {
                    for (size_t ty = 0; ty < tH; ++ty)
                        for (size_t tx = 0; tx < tW; ++tx)
                            val += imgCh(y + ty, x + tx) * tplCh(ty, tx);
                } else {
                    // CCOEFF: Σ (I(y+ty,x+tx) - I_mean) * (T(ty,tx) - T_mean)
                    T imgMean = T{0};
                    for (size_t ty = 0; ty < tH; ++ty)
                        for (size_t tx = 0; tx < tW; ++tx)
                            imgMean += imgCh(y + ty, x + tx);
                    imgMean /= static_cast<T>(tH * tW);

                    for (size_t ty = 0; ty < tH; ++ty)
                        for (size_t tx = 0; tx < tW; ++tx)
                            val += (imgCh(y + ty, x + tx) - imgMean)
                                 * (tplCh(ty, tx) - tplMean);
                }

                result(y, x) = val;
            }
        }
        return result;
    }

    /// 結果マップから最小/最大位置を取得
    /// @param result   マッチング結果マップ
    /// @param findMax  true: 最大位置, false: 最小位置
    /// @return         (y, x) 位置
    template<typename T>
    std::pair<size_t, size_t> minMaxLoc(const Matrix<T>& result, bool findMax = true) {
        assert(result.rows() > 0 && result.cols() > 0);

        size_t bestY = 0, bestX = 0;
        T bestVal = result(0, 0);

        for (size_t y = 0; y < result.rows(); ++y) {
            for (size_t x = 0; x < result.cols(); ++x) {
                T v = result(y, x);
                if (findMax ? (v > bestVal) : (v < bestVal)) {
                    bestVal = v;
                    bestY = y;
                    bestX = x;
                }
            }
        }
        return {bestY, bestX};
    }

} // namespace calx::image

#endif // CALX_IMAGE_IMGPROC_TEMPLATE_MATCHING_HPP
