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

// morphology.hpp — モルフォロジー演算 (scipy.ndimage 対応)
//
// 新規実装:
//   - 膨張 (dilation), 収縮 (erosion)
//   - オープニング, クロージング
//   - ラベリング (connected components), 領域計測

#ifndef CALX_IMAGE_MORPHOLOGY_HPP
#define CALX_IMAGE_MORPHOLOGY_HPP

#include <math/core/matrix.hpp>
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <vector>

namespace calx::image {

    /// 構造要素 (structuring element) を生成
    /// @param size   サイズ (奇数)
    /// @param cross  true: 十字型, false: 正方形
    inline Matrix<int> structuringElement(int size, bool cross = false) {
        assert(size > 0 && size % 2 == 1);
        Matrix<int> se(size, size, cross ? 0 : 1);
        if (cross) {
            int half = size / 2;
            for (int i = 0; i < size; ++i) {
                se(half, i) = 1;
                se(i, half) = 1;
            }
        }
        return se;
    }

    /// 膨張 (dilation)
    /// @param image   2 値画像 (0/非0)
    /// @param se      構造要素
    template<typename T>
    Matrix<T> dilate(const Matrix<T>& image, const Matrix<int>& se) {
        int H = static_cast<int>(image.rows());
        int W = static_cast<int>(image.cols());
        int kH = static_cast<int>(se.rows());
        int kW = static_cast<int>(se.cols());
        int hkH = kH / 2;
        int hkW = kW / 2;

        Matrix<T> result(H, W, T{0});

        for (int y = 0; y < H; ++y) {
            for (int x = 0; x < W; ++x) {
                T maxVal = T{0};
                for (int ky = 0; ky < kH; ++ky) {
                    for (int kx = 0; kx < kW; ++kx) {
                        if (se(ky, kx) == 0) continue;
                        int iy = y + ky - hkH;
                        int ix = x + kx - hkW;
                        if (iy >= 0 && iy < H && ix >= 0 && ix < W) {
                            if (image(iy, ix) > maxVal) maxVal = image(iy, ix);
                        }
                    }
                }
                result(y, x) = maxVal;
            }
        }
        return result;
    }

    /// 収縮 (erosion)
    template<typename T>
    Matrix<T> erode(const Matrix<T>& image, const Matrix<int>& se) {
        int H = static_cast<int>(image.rows());
        int W = static_cast<int>(image.cols());
        int kH = static_cast<int>(se.rows());
        int kW = static_cast<int>(se.cols());
        int hkH = kH / 2;
        int hkW = kW / 2;

        Matrix<T> result(H, W, T{0});

        for (int y = 0; y < H; ++y) {
            for (int x = 0; x < W; ++x) {
                T minVal = std::numeric_limits<T>::max();
                bool hasNeighbor = false;
                for (int ky = 0; ky < kH; ++ky) {
                    for (int kx = 0; kx < kW; ++kx) {
                        if (se(ky, kx) == 0) continue;
                        int iy = y + ky - hkH;
                        int ix = x + kx - hkW;
                        if (iy >= 0 && iy < H && ix >= 0 && ix < W) {
                            if (image(iy, ix) < minVal) minVal = image(iy, ix);
                            hasNeighbor = true;
                        } else {
                            minVal = T{0};
                            hasNeighbor = true;
                        }
                    }
                }
                result(y, x) = hasNeighbor ? minVal : T{0};
            }
        }
        return result;
    }

    /// オープニング (erosion → dilation)
    template<typename T>
    Matrix<T> opening(const Matrix<T>& image, const Matrix<int>& se) {
        return dilate(erode(image, se), se);
    }

    /// クロージング (dilation → erosion)
    template<typename T>
    Matrix<T> closing(const Matrix<T>& image, const Matrix<int>& se) {
        return erode(dilate(image, se), se);
    }

    /// 連結成分ラベリング (4 連結)
    /// @param image   2 値画像 (0 = 背景, 非0 = 前景)
    /// @return        ラベル行列 (0 = 背景, 1,2,... = 各連結成分)
    template<typename T>
    Matrix<int> labelComponents(const Matrix<T>& image) {
        int H = static_cast<int>(image.rows());
        int W = static_cast<int>(image.cols());
        Matrix<int> labels(H, W, 0);

        // Union-Find
        std::vector<int> parent;
        parent.push_back(0);  // 背景

        auto find = [&](int x) -> int {
            while (parent[x] != x) {
                parent[x] = parent[parent[x]];
                x = parent[x];
            }
            return x;
        };

        auto unite = [&](int a, int b) {
            a = find(a);
            b = find(b);
            if (a != b) {
                if (a > b) std::swap(a, b);
                parent[b] = a;
            }
        };

        int nextLabel = 1;

        // 1st pass
        for (int y = 0; y < H; ++y) {
            for (int x = 0; x < W; ++x) {
                if (image(y, x) == T{0}) continue;

                int up = (y > 0 && image(y - 1, x) != T{0}) ? labels(y - 1, x) : 0;
                int left = (x > 0 && image(y, x - 1) != T{0}) ? labels(y, x - 1) : 0;

                if (up == 0 && left == 0) {
                    labels(y, x) = nextLabel;
                    parent.push_back(nextLabel);
                    ++nextLabel;
                } else if (up != 0 && left == 0) {
                    labels(y, x) = up;
                } else if (up == 0 && left != 0) {
                    labels(y, x) = left;
                } else {
                    labels(y, x) = std::min(up, left);
                    unite(up, left);
                }
            }
        }

        // 2nd pass: ラベルを正規化
        std::vector<int> remap(nextLabel, 0);
        int finalLabel = 0;
        for (int i = 1; i < nextLabel; ++i) {
            int root = find(i);
            if (remap[root] == 0) {
                remap[root] = ++finalLabel;
            }
            remap[i] = remap[root];
        }

        for (int y = 0; y < H; ++y)
            for (int x = 0; x < W; ++x)
                if (labels(y, x) > 0)
                    labels(y, x) = remap[labels(y, x)];

        return labels;
    }

    /// ラベルごとの領域計測結果
    struct RegionProps {
        int label;
        int area;               // ピクセル数
        double centroidX;       // 重心 X
        double centroidY;       // 重心 Y
    };

    /// ラベル行列から各領域の計測値を求める
    inline std::vector<RegionProps> measureRegions(const Matrix<int>& labels) {
        int H = static_cast<int>(labels.rows());
        int W = static_cast<int>(labels.cols());

        // 最大ラベルを求める
        int maxLabel = 0;
        for (int y = 0; y < H; ++y)
            for (int x = 0; x < W; ++x)
                if (labels(y, x) > maxLabel) maxLabel = labels(y, x);

        std::vector<int> area(maxLabel + 1, 0);
        std::vector<double> sumX(maxLabel + 1, 0);
        std::vector<double> sumY(maxLabel + 1, 0);

        for (int y = 0; y < H; ++y) {
            for (int x = 0; x < W; ++x) {
                int l = labels(y, x);
                if (l > 0) {
                    ++area[l];
                    sumX[l] += x;
                    sumY[l] += y;
                }
            }
        }

        std::vector<RegionProps> props;
        for (int i = 1; i <= maxLabel; ++i) {
            if (area[i] > 0) {
                props.push_back({i, area[i],
                                 sumX[i] / area[i], sumY[i] / area[i]});
            }
        }
        return props;
    }

} // namespace calx::image

#endif // CALX_IMAGE_MORPHOLOGY_HPP
