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

/**
 * @file hdr.hpp
 * @brief HDR 合成 + トーンマッピング (OCV-25)
 *
 * 複数露出画像からの HDR 合成 (Debevec) + トーンマッピング (Reinhard, Drago)
 */

#ifndef CALX_IMAGE_IMGPROC_HDR_HPP
#define CALX_IMAGE_IMGPROC_HDR_HPP

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

namespace calx::image {

/// 複数露出画像を HDR に合成 (簡易版: 重み付き平均)
template<typename T>
Image<T> mergeExposures(const std::vector<Image<T>>& images,
                         const std::vector<double>& exposureTimes) {
    if (images.empty()) return {};
    int h = static_cast<int>(images[0].height());
    int w = static_cast<int>(images[0].width());
    int ch = images[0].numChannels();
    int n = static_cast<int>(images.size());

    Image<T> hdr(h, w, ch, T(0));

    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x)
            for (int c = 0; c < ch; ++c) {
                double sumW = 0, sumVal = 0;
                for (int i = 0; i < n; ++i) {
                    double v = double(images[i].at(y, x, c));
                    // 重み: 中間値に近いほど信頼性が高い (hat function)
                    double w = 1.0 - std::abs(2.0 * v - 1.0);
                    w = std::max(w, 0.01);
                    // log ドメインで合成
                    double logE = std::log(v + 1e-6) - std::log(exposureTimes[i]);
                    sumVal += w * logE;
                    sumW += w;
                }
                hdr.at(y, x, c) = T(std::exp(sumVal / sumW));
            }
    return hdr;
}

/// Reinhard トーンマッピング
template<typename T>
Image<T> tonemapReinhard(const Image<T>& hdr, T gamma = T(1), T intensity = T(0),
                          T lightAdapt = T(1), T colorAdapt = T(0)) {
    int h = static_cast<int>(hdr.height());
    int w = static_cast<int>(hdr.width());
    int ch = hdr.numChannels();

    // 対数平均輝度
    double sumLog = 0;
    int count = 0;
    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x) {
            double lum = 0;
            for (int c = 0; c < ch; ++c) lum += double(hdr.at(y, x, c));
            lum /= ch;
            sumLog += std::log(lum + 1e-6);
            ++count;
        }
    double logAvg = std::exp(sumLog / count);
    double key = 0.18 * std::pow(2.0, double(intensity));

    Image<T> result(h, w, ch);
    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x)
            for (int c = 0; c < ch; ++c) {
                double v = double(hdr.at(y, x, c)) * key / logAvg;
                v = v / (1.0 + v);  // Reinhard operator
                v = std::pow(v, 1.0 / double(gamma));
                result.at(y, x, c) = T(std::clamp(v, 0.0, 1.0));
            }
    return result;
}

/// Drago トーンマッピング
template<typename T>
Image<T> tonemapDrago(const Image<T>& hdr, T gamma = T(1), T bias = T(0.85)) {
    int h = static_cast<int>(hdr.height());
    int w = static_cast<int>(hdr.width());
    int ch = hdr.numChannels();

    double maxLum = 0;
    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x) {
            double lum = 0;
            for (int c = 0; c < ch; ++c) lum += double(hdr.at(y, x, c));
            maxLum = std::max(maxLum, lum / ch);
        }

    Image<T> result(h, w, ch);
    double b = double(bias);
    double logBase = std::log(2.0 + 8.0 * std::pow(maxLum / (maxLum + 1e-6), std::log(b) / std::log(0.5)));

    for (int y = 0; y < h; ++y)
        for (int x = 0; x < w; ++x)
            for (int c = 0; c < ch; ++c) {
                double v = double(hdr.at(y, x, c));
                double mapped = std::log(v + 1.0) / logBase;
                mapped = std::pow(mapped, 1.0 / double(gamma));
                result.at(y, x, c) = T(std::clamp(mapped, 0.0, 1.0));
            }
    return result;
}

} // namespace calx::image

#endif
