102 lines
3.9 KiB
C++
102 lines
3.9 KiB
C++
#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
|