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

/**
 * @file seamless_clone.hpp
 * @brief Seamless Cloning / Poisson Blending (OCV-26)
 *
 * 前景画像の勾配を維持しつつ背景画像になじませる。
 * Poisson 方程式を Jacobi 反復で解く。
 */

#ifndef CALX_IMAGE_IMGPROC_SEAMLESS_CLONE_HPP
#define CALX_IMAGE_IMGPROC_SEAMLESS_CLONE_HPP

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

namespace calx::image {

/// Seamless Clone (Poisson blending)
/// src: 前景画像, dst: 背景画像, mask: 合成領域 (非0), cx/cy: dst 上の貼付位置
template<typename T>
Image<T> seamlessClone(const Image<T>& src, const Image<T>& dst,
                        const Image<T>& mask, int cx, int cy,
                        int iterations = 500) {
    int sh = static_cast<int>(src.height()), sw = static_cast<int>(src.width());
    int dh = static_cast<int>(dst.height()), dw = static_cast<int>(dst.width());
    int ch = std::min(src.numChannels(), dst.numChannels());
    int ox = cx - sw / 2, oy = cy - sh / 2;  // オフセット

    Image<T> result = dst;

    // マスク領域の判定
    auto inMask = [&](int sy, int sx) -> bool {
        if (sy < 0 || sy >= sh || sx < 0 || sx >= sw) return false;
        return mask(sy, sx) != T(0);
    };

    // src の勾配 (ラプラシアン) を計算
    // ∇²f = f(y-1,x) + f(y+1,x) + f(x-1,y) + f(x+1,y) - 4f(x,y)
    auto srcLap = [&](int sy, int sx, int c) -> double {
        double v = 0;
        if (sy > 0)    v += double(src.at(sy-1, sx, c));
        if (sy < sh-1) v += double(src.at(sy+1, sx, c));
        if (sx > 0)    v += double(src.at(sy, sx-1, c));
        if (sx < sw-1) v += double(src.at(sy, sx+1, c));
        v -= 4.0 * double(src.at(sy, sx, c));
        return v;
    };

    // Jacobi 反復で Poisson 方程式を解く
    for (int iter = 0; iter < iterations; ++iter) {
        for (int sy = 1; sy < sh - 1; ++sy) {
            for (int sx = 1; sx < sw - 1; ++sx) {
                if (!inMask(sy, sx)) continue;
                int dy = oy + sy, dx = ox + sx;
                if (dy < 0 || dy >= dh || dx < 0 || dx >= dw) continue;

                for (int c = 0; c < ch; ++c) {
                    double lap = srcLap(sy, sx, c);
                    // 隣接ピクセルの現在値
                    double sum = 0;
                    int cnt = 0;
                    for (auto [ndy, ndx] : std::vector<std::pair<int,int>>{{-1,0},{1,0},{0,-1},{0,1}}) {
                        int ny = dy + ndy, nx = dx + ndx;
                        if (ny < 0 || ny >= dh || nx < 0 || nx >= dw) continue;
                        sum += double(result.at(ny, nx, c));
                        ++cnt;
                    }
                    if (cnt > 0) {
                        double v = (sum + lap) / cnt;
                        result.at(dy, dx, c) = T(std::clamp(v, 0.0, 1.0));
                    }
                }
            }
        }
    }
    return result;
}

} // namespace calx::image

#endif
