#include "VoxelFromScatters.hpp"

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>

#include "actors/VoxelActor.hpp"
#include "algo/IdwInterpolator.hpp"

namespace geopro::render {

namespace {

// 体素网格各维上限：防止退化输入(如坐标异常)撑出超大体拖垮内存/渲染。
constexpr int kMaxDim = 400;

// 钳制维度到 [1, kMaxDim]。
int clampDim(double ext, double cell)
{
    int n = static_cast<int>(ext / cell) + 1;
    if (n < 1) n = 1;
    if (n > kMaxDim) n = kMaxDim;
    return n;
}

}  // namespace

VoxelResult buildVoxelFromScatters(const std::vector<geopro::core::ScatterField>& profiles,
                                   const geopro::core::ColorScale& cs,
                                   const geopro::core::CrsTransform& crs,
                                   const geopro::core::GeoLocalFrame& frame,
                                   double cellXY,
                                   double cellZ,
                                   double power,
                                   double maxDist,
                                   double zDisplayScale)
{
    // 1) 配准所有点到世界局部米 + 深度，组装 IDW 输入点集。
    geopro::core::PointSet pts;
    for (const auto& s : profiles) {
        const std::size_t n = s.v.size();
        if (s.projX.size() < n || s.projY.size() < n || s.y.size() < n) continue;  // 字段不齐跳过
        for (std::size_t i = 0; i < n; ++i) {
            const auto ll = crs.forward(s.projX[i], s.projY[i]);  // (lon, lat)，always_xy
            const auto local = frame.toLocal(ll.y, ll.x);         // (x East, y North) 米
            pts.x.push_back(local.x);
            pts.y.push_back(local.y);
            pts.z.push_back(-s.y[i]);  // 深度向下：z 取负(与帘面一致)
            pts.v.push_back(s.v[i]);
        }
    }
    if (pts.v.empty()) return VoxelResult{};

    // 2) 由点集包络定 GridSpec(角点对齐，与 IDW/vtkImageData 一致)。
    double minx = pts.x[0], maxx = pts.x[0];
    double miny = pts.y[0], maxy = pts.y[0];
    double minz = pts.z[0], maxz = pts.z[0];
    for (std::size_t i = 1; i < pts.v.size(); ++i) {
        minx = std::min(minx, pts.x[i]); maxx = std::max(maxx, pts.x[i]);
        miny = std::min(miny, pts.y[i]); maxy = std::max(maxy, pts.y[i]);
        minz = std::min(minz, pts.z[i]); maxz = std::max(maxz, pts.z[i]);
    }

    geopro::core::GridSpec spec{};
    spec.ox = minx; spec.oy = miny; spec.oz = minz;
    spec.dx = cellXY; spec.dy = cellXY; spec.dz = cellZ;
    spec.nx = clampDim(maxx - minx, cellXY);
    spec.ny = clampDim(maxy - miny, cellXY);
    spec.nz = clampDim(maxz - minz, cellZ);
    spec.power = power;
    spec.maxDist = maxDist;

    // 3) IDW → ScalarVolume(maxDist 外 NaN 留空)。
    const geopro::core::IdwInterpolator idw;
    const geopro::core::ScalarVolume vol = idw.interpolate(pts, spec);

    // 4) 色阶范围：优先 colorBar 真实分段值，否则数据实测。
    double vmin, vmax;
    const std::vector<double> stops = cs.stopValues();
    if (stops.size() >= 2) {
        vmin = stops.front(); vmax = stops.back();
    } else {
        vmin = std::numeric_limits<double>::infinity();
        vmax = -std::numeric_limits<double>::infinity();
        for (double v : vol.data()) {
            if (std::isnan(v)) continue;
            vmin = std::min(vmin, v); vmax = std::max(vmax, v);
        }
        if (!(vmin < vmax)) { vmin = 0.0; vmax = 1.0; }
    }

    // 5) 体绘制 + 暴露 image(供切片)。纵向夸张烤进 image 的 z 原点/间距(IDW 已用真实 cellZ 采样)，
    //    使体绘制/切片/帘面纵向一致(切片穿过体素而非在旁)。
    VoxelResult out;
    out.volume = buildVoxel(vol, cs, spec.ox, spec.oy, spec.oz * zDisplayScale, spec.dx, spec.dy,
                            spec.dz * zDisplayScale, vmin, vmax, out.image);
    return out;
}

}  // namespace geopro::render