72 lines
2.6 KiB
C++
72 lines
2.6 KiB
C++
#include "algo/VolumeBuilder.hpp"
|
||
|
||
#include <algorithm>
|
||
#include <cmath>
|
||
#include <cstddef>
|
||
#include <limits>
|
||
#include <stdexcept>
|
||
|
||
#include "algo/IdwInterpolator.hpp"
|
||
|
||
namespace geopro::core {
|
||
|
||
namespace {
|
||
// 某轴：优先用 cell 间距；若包络 ext 过大致格数超 kMaxVolumeDim，则**放大间距**使 maxDim 格跨满 ext
|
||
// （分辨率降低，但**不截断**——否则跨 TM 多剖面相距 > maxDim×cell 时，远端剖面落网格外、丢失）。
|
||
void fitAxis(double ext, double cell, double& outCell, int& outN) {
|
||
if (!(ext > 0.0) || !(cell > 0.0)) { outCell = (cell > 0.0 ? cell : 1.0); outN = 1; return; }
|
||
int n = static_cast<int>(ext / cell) + 1;
|
||
if (n <= kMaxVolumeDim) {
|
||
outCell = cell;
|
||
outN = (n < 1) ? 1 : n;
|
||
return;
|
||
}
|
||
outN = kMaxVolumeDim;
|
||
outCell = ext / static_cast<double>(kMaxVolumeDim - 1); // maxDim 格覆盖全 ext
|
||
}
|
||
} // namespace
|
||
|
||
BuiltVolume buildVolume(const PointSet& pts, double cellXY, double cellZ,
|
||
double power, double maxDist) {
|
||
if (pts.v.empty()) {
|
||
throw std::invalid_argument("buildVolume: empty point set");
|
||
}
|
||
|
||
// 1) 点集包络。
|
||
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]);
|
||
}
|
||
|
||
// 2) GridSpec（角点对齐 = 原点取包络最小角）。间距优先用 cell；包络过大时放大间距以覆盖全程
|
||
// （fitAxis），避免跨 TM 多剖面相距过远时远端被截断。
|
||
GridSpec spec{};
|
||
spec.ox = minx; spec.oy = miny; spec.oz = minz;
|
||
fitAxis(maxx - minx, cellXY, spec.dx, spec.nx);
|
||
fitAxis(maxy - miny, cellXY, spec.dy, spec.ny);
|
||
fitAxis(maxz - minz, cellZ, spec.dz, spec.nz);
|
||
spec.power = power;
|
||
spec.maxDist = maxDist;
|
||
|
||
// 3) IDW（maxDist 外 NaN 留空）。
|
||
const IdwInterpolator idw;
|
||
ScalarVolume vol = idw.interpolate(pts, spec);
|
||
|
||
// 4) 数据实测值域（仅有限值）。无有限值 → 退化 {0,1}。
|
||
double vmin = std::numeric_limits<double>::infinity();
|
||
double 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; }
|
||
|
||
return BuiltVolume{std::move(vol), spec, vmin, vmax};
|
||
}
|
||
|
||
} // namespace geopro::core
|