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

/**
 * @file moments.hpp
 * @brief 画像モーメント (OCV-9: moments / HuMoments)
 *
 * 空間モーメント, 中心モーメント, 正規化モーメント, Hu 不変モーメント
 */

#ifndef CALX_IMAGE_IMGPROC_MOMENTS_HPP
#define CALX_IMAGE_IMGPROC_MOMENTS_HPP

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

namespace calx::image {

/// 画像モーメント
struct ImageMoments {
    // 空間モーメント m_{pq} = ΣΣ x^p y^q I(x,y)
    double m00=0, m10=0, m01=0, m20=0, m11=0, m02=0, m30=0, m21=0, m12=0, m03=0;

    // 中心モーメント μ_{pq} = ΣΣ (x-x̄)^p (y-ȳ)^q I(x,y)
    double mu20=0, mu11=0, mu02=0, mu30=0, mu21=0, mu12=0, mu03=0;

    // 正規化中心モーメント η_{pq} = μ_{pq} / m00^((p+q)/2+1)
    double nu20=0, nu11=0, nu02=0, nu30=0, nu21=0, nu12=0, nu03=0;

    // 重心
    double xBar() const { return (m00 > 0) ? m10 / m00 : 0; }
    double yBar() const { return (m00 > 0) ? m01 / m00 : 0; }
};

/// 画像モーメントを計算
template<typename T>
ImageMoments moments(const Image<T>& img) {
    ImageMoments m;
    int h = static_cast<int>(img.height());
    int w = static_cast<int>(img.width());

    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x) {
            double I = static_cast<double>(img(y, x));
            double xx = x, yy = y;
            m.m00 += I;
            m.m10 += xx * I;       m.m01 += yy * I;
            m.m20 += xx*xx * I;    m.m11 += xx*yy * I;    m.m02 += yy*yy * I;
            m.m30 += xx*xx*xx * I; m.m21 += xx*xx*yy * I;
            m.m12 += xx*yy*yy * I; m.m03 += yy*yy*yy * I;
        }

    // 中心モーメント
    double cx = m.xBar(), cy = m.yBar();
    m.mu20 = m.m20 - cx * m.m10;
    m.mu11 = m.m11 - cx * m.m01;
    m.mu02 = m.m02 - cy * m.m01;
    m.mu30 = m.m30 - 3*cx*m.m20 + 2*cx*cx*m.m10;
    m.mu21 = m.m21 - 2*cx*m.m11 - cy*m.m20 + 2*cx*cx*m.m01;
    m.mu12 = m.m12 - 2*cy*m.m11 - cx*m.m02 + 2*cy*cy*m.m10;
    m.mu03 = m.m03 - 3*cy*m.m02 + 2*cy*cy*m.m01;

    // 正規化
    if (m.m00 > 0) {
        auto norm = [&](double mu, int p, int q) {
            return mu / std::pow(m.m00, double(p + q) / 2.0 + 1.0);
        };
        m.nu20 = norm(m.mu20, 2, 0); m.nu11 = norm(m.mu11, 1, 1);
        m.nu02 = norm(m.mu02, 0, 2);
        m.nu30 = norm(m.mu30, 3, 0); m.nu21 = norm(m.mu21, 2, 1);
        m.nu12 = norm(m.mu12, 1, 2); m.nu03 = norm(m.mu03, 0, 3);
    }
    return m;
}

/// Hu の 7 つの不変モーメント (回転・スケール・平行移動に不変)
inline std::array<double, 7> huMoments(const ImageMoments& m) {
    double n20=m.nu20, n02=m.nu02, n11=m.nu11;
    double n30=m.nu30, n21=m.nu21, n12=m.nu12, n03=m.nu03;
    std::array<double, 7> hu;
    hu[0] = n20 + n02;
    hu[1] = (n20 - n02)*(n20 - n02) + 4*n11*n11;
    hu[2] = (n30 - 3*n12)*(n30 - 3*n12) + (3*n21 - n03)*(3*n21 - n03);
    hu[3] = (n30 + n12)*(n30 + n12) + (n21 + n03)*(n21 + n03);
    hu[4] = (n30-3*n12)*(n30+n12)*((n30+n12)*(n30+n12)-3*(n21+n03)*(n21+n03))
           +(3*n21-n03)*(n21+n03)*(3*(n30+n12)*(n30+n12)-(n21+n03)*(n21+n03));
    hu[5] = (n20-n02)*((n30+n12)*(n30+n12)-(n21+n03)*(n21+n03))
           +4*n11*(n30+n12)*(n21+n03);
    hu[6] = (3*n21-n03)*(n30+n12)*((n30+n12)*(n30+n12)-3*(n21+n03)*(n21+n03))
           -(n30-3*n12)*(n21+n03)*(3*(n30+n12)*(n30+n12)-(n21+n03)*(n21+n03));
    return hu;
}

} // namespace calx::image

#endif
