336 lines
12 KiB
C++
336 lines
12 KiB
C++
// GeoVolumeBuilder：按真实 RTK 几何把多条测线插值进统一路向网格。
|
||
//
|
||
// 编排层（io_gpr + core 采样核 + store），故与 StreamingVolumeBuilder 一样编进
|
||
// geopro_store（不污染纯 geopro_core）。命名空间保持 geopro::core（接口契约）。
|
||
#include "core/algo/GeoVolumeBuilder.hpp"
|
||
|
||
#include <algorithm>
|
||
#include <array>
|
||
#include <cmath>
|
||
#include <cstdint>
|
||
#include <filesystem>
|
||
#include <fstream>
|
||
#include <iterator>
|
||
#include <limits>
|
||
#include <stdexcept>
|
||
#include <vector>
|
||
|
||
#include "core/algo/GprVolumeBuilder.hpp" // BuiltI16
|
||
#include "core/model/ScalarVolumeI16.hpp" // ScalarVolumeI16::kBlank, Quant
|
||
#include "data/store/ChunkedVolumeStore.hpp"
|
||
#include "io/gpr/GprGeometry.hpp" // parseChannelXOffsets, depthOfSample
|
||
#include "io/gpr/GpsTrack.hpp"
|
||
#include "io/gpr/IprHeader.hpp"
|
||
#include "io/gpr/IprbReader.hpp"
|
||
|
||
namespace fs = std::filesystem;
|
||
|
||
namespace geopro::core {
|
||
|
||
namespace {
|
||
|
||
constexpr double kPi = 3.14159265358979323846;
|
||
|
||
// 读 .iprh 文本 → 解析头（与 .iprb 同名）。
|
||
geopro::io::gpr::IprHeader readHeaderFor(const std::string& iprbPath) {
|
||
fs::path h = fs::path(iprbPath).replace_extension(".iprh");
|
||
std::ifstream f(h);
|
||
if (!f) throw std::runtime_error("GeoVolumeBuilder: 打不开 iprh " + h.string());
|
||
std::string text((std::istreambuf_iterator<char>(f)),
|
||
std::istreambuf_iterator<char>());
|
||
return geopro::io::gpr::parseIprHeader(text);
|
||
}
|
||
|
||
// 一条线的标尺（轨迹局部米 + 通道横偏 + 头）。
|
||
struct LineScale {
|
||
std::vector<geopro::io::gpr::XY> trackM; // 该线轨迹（局部米，旋转前）
|
||
std::vector<double> chanX; // 通道横偏（米，文件顺序）
|
||
geopro::io::gpr::IprHeader header;
|
||
std::int64_t totalTraces = 0;
|
||
};
|
||
|
||
// 全线总道数 = min 通道(fileBytes/(samples*2))（与 assembler 口径一致）。
|
||
std::int64_t totalTracesOf(const std::vector<std::string>& iprb, int samples) {
|
||
std::int64_t minTr = std::numeric_limits<std::int64_t>::max();
|
||
const std::int64_t per = static_cast<std::int64_t>(samples) * 2;
|
||
if (per <= 0) throw std::runtime_error("samples<=0");
|
||
for (const auto& p : iprb) {
|
||
const std::int64_t bytes = static_cast<std::int64_t>(fs::file_size(p));
|
||
minTr = std::min(minTr, bytes / per);
|
||
}
|
||
return minTr;
|
||
}
|
||
|
||
} // namespace
|
||
|
||
// ---- PCA 主轴角（可测纯函数）----
|
||
double principalAxisAngle(const std::vector<double>& xs,
|
||
const std::vector<double>& ys) {
|
||
const std::size_t n = std::min(xs.size(), ys.size());
|
||
if (n < 2) return 0.0;
|
||
|
||
double mx = 0, my = 0;
|
||
for (std::size_t i = 0; i < n; ++i) {
|
||
mx += xs[i];
|
||
my += ys[i];
|
||
}
|
||
mx /= static_cast<double>(n);
|
||
my /= static_cast<double>(n);
|
||
|
||
double sxx = 0, syy = 0, sxy = 0;
|
||
for (std::size_t i = 0; i < n; ++i) {
|
||
const double dx = xs[i] - mx, dy = ys[i] - my;
|
||
sxx += dx * dx;
|
||
syy += dy * dy;
|
||
sxy += dx * dy;
|
||
}
|
||
sxx /= static_cast<double>(n);
|
||
syy /= static_cast<double>(n);
|
||
sxy /= static_cast<double>(n);
|
||
|
||
// 2×2 对称协方差矩阵最大特征向量方向 = 0.5*atan2(2*sxy, sxx-syy)。
|
||
if (std::abs(sxy) < 1e-15 && std::abs(sxx - syy) < 1e-15) return 0.0;
|
||
return 0.5 * std::atan2(2.0 * sxy, sxx - syy);
|
||
}
|
||
|
||
GeoBuildResult buildGeoVolume(const std::vector<GeoLineInput>& lines,
|
||
const GeoGridSpec& spec,
|
||
const std::string& outStoreDir,
|
||
int pyramidLevels) {
|
||
using geopro::io::gpr::interpAlongTrack;
|
||
using geopro::io::gpr::lonLatToLocalM;
|
||
using geopro::io::gpr::PosHeading;
|
||
using geopro::io::gpr::XY;
|
||
|
||
if (lines.empty()) throw std::runtime_error("buildGeoVolume: 无测线");
|
||
if (spec.cellXY <= 0 || spec.cellZ <= 0)
|
||
throw std::runtime_error("buildGeoVolume: cell 必须 > 0");
|
||
|
||
// ---- 1) 各线 .gps → 局部米（共用原点 = 全体最小 lat/lon）----
|
||
std::vector<geopro::io::gpr::GpsTrack> tracks(lines.size());
|
||
double minLat = std::numeric_limits<double>::infinity();
|
||
double minLon = std::numeric_limits<double>::infinity();
|
||
for (std::size_t i = 0; i < lines.size(); ++i) {
|
||
tracks[i] = geopro::io::gpr::parseGps(lines[i].gps);
|
||
for (const auto& p : tracks[i].pts) {
|
||
minLat = std::min(minLat, p.lat);
|
||
minLon = std::min(minLon, p.lon);
|
||
}
|
||
}
|
||
if (!std::isfinite(minLat) || !std::isfinite(minLon))
|
||
throw std::runtime_error("buildGeoVolume: 轨迹为空");
|
||
|
||
// 各线标尺：轨迹局部米 + 通道横偏 + 头 + 总道数。
|
||
std::vector<LineScale> scales(lines.size());
|
||
int samples = 0;
|
||
geopro::io::gpr::IprHeader hdr0{};
|
||
for (std::size_t i = 0; i < lines.size(); ++i) {
|
||
LineScale& sc = scales[i];
|
||
sc.trackM.reserve(tracks[i].pts.size());
|
||
for (const auto& p : tracks[i].pts)
|
||
sc.trackM.push_back(lonLatToLocalM(p.lat, p.lon, minLat, minLon));
|
||
|
||
std::ifstream ordF(lines[i].ord);
|
||
if (!ordF) throw std::runtime_error("buildGeoVolume: 打不开 ord " + lines[i].ord);
|
||
std::string ordText((std::istreambuf_iterator<char>(ordF)),
|
||
std::istreambuf_iterator<char>());
|
||
sc.chanX = geopro::io::gpr::parseChannelXOffsets(ordText);
|
||
|
||
if (lines[i].iprb.empty())
|
||
throw std::runtime_error("buildGeoVolume: 线无 iprb");
|
||
sc.header = readHeaderFor(lines[i].iprb.front());
|
||
sc.totalTraces = totalTracesOf(lines[i].iprb, sc.header.samples);
|
||
if (i == 0) {
|
||
samples = sc.header.samples;
|
||
hdr0 = sc.header;
|
||
}
|
||
}
|
||
if (samples <= 0) throw std::runtime_error("buildGeoVolume: samples<=0");
|
||
|
||
// ---- 2) PCA 求路向主轴角，旋转使路向 = X'、横路 = Y' ----
|
||
// 收集全体通道平面点（trace 位置 + 通道横偏）求主轴。
|
||
std::vector<double> px, py;
|
||
for (std::size_t i = 0; i < scales.size(); ++i) {
|
||
const LineScale& sc = scales[i];
|
||
const std::int64_t nt = sc.totalTraces;
|
||
if (nt < 1 || sc.trackM.size() < 2) continue;
|
||
// 稀疏采样道（每线最多 ~200 个采样点）求 PCA，省内存。
|
||
const std::int64_t stride = std::max<std::int64_t>(1, nt / 200);
|
||
for (std::int64_t t = 0; t < nt; t += stride) {
|
||
const double frac = nt > 1 ? static_cast<double>(t) / (nt - 1) : 0.0;
|
||
PosHeading ph = interpAlongTrack(sc.trackM, frac);
|
||
px.push_back(ph.pos.x);
|
||
py.push_back(ph.pos.y);
|
||
}
|
||
}
|
||
const double rotRad = principalAxisAngle(px, py);
|
||
const double cosR = std::cos(-rotRad), sinR = std::sin(-rotRad);
|
||
auto rotate = [&](double x, double y, double& rx, double& ry) {
|
||
rx = x * cosR - y * sinR;
|
||
ry = x * sinR + y * cosR;
|
||
};
|
||
|
||
// ---- 求旋转后包围盒（含通道横偏）→ 网格维度 + 原点 ----
|
||
double minX = std::numeric_limits<double>::infinity();
|
||
double minY = std::numeric_limits<double>::infinity();
|
||
double maxX = -std::numeric_limits<double>::infinity();
|
||
double maxY = -std::numeric_limits<double>::infinity();
|
||
for (std::size_t i = 0; i < scales.size(); ++i) {
|
||
const LineScale& sc = scales[i];
|
||
const std::int64_t nt = sc.totalTraces;
|
||
if (nt < 1 || sc.trackM.size() < 2 || sc.chanX.empty()) continue;
|
||
const std::int64_t stride = std::max<std::int64_t>(1, nt / 200);
|
||
for (std::int64_t t = 0; t < nt; t += stride) {
|
||
const double frac = nt > 1 ? static_cast<double>(t) / (nt - 1) : 0.0;
|
||
PosHeading ph = interpAlongTrack(sc.trackM, frac);
|
||
const double perpX = -ph.hy, perpY = ph.hx; // 垂直航向（左法向）
|
||
for (double oc : {sc.chanX.front(), sc.chanX.back()}) {
|
||
const double cx = ph.pos.x + oc * perpX;
|
||
const double cy = ph.pos.y + oc * perpY;
|
||
double rx, ry;
|
||
rotate(cx, cy, rx, ry);
|
||
minX = std::min(minX, rx);
|
||
maxX = std::max(maxX, rx);
|
||
minY = std::min(minY, ry);
|
||
maxY = std::max(maxY, ry);
|
||
}
|
||
}
|
||
}
|
||
if (!std::isfinite(minX)) throw std::runtime_error("buildGeoVolume: 无有效几何");
|
||
|
||
// 深度范围：用首线头取最深样本深度。
|
||
const double maxDepth = geopro::io::gpr::depthOfSample(samples - 1, hdr0);
|
||
|
||
auto cells = [](double range, double cell) {
|
||
if (cell <= 0.0 || range <= 0.0) return 1;
|
||
return static_cast<int>(std::ceil(range / cell)) + 1;
|
||
};
|
||
const int nx = cells(maxX - minX, spec.cellXY); // 沿路（长轴 X'）
|
||
const int ny = cells(maxY - minY, spec.cellXY); // 横路 Y'
|
||
const int nz = cells(maxDepth, spec.cellZ); // 深度 Z
|
||
|
||
const std::size_t cellCount =
|
||
static_cast<std::size_t>(nx) * ny * nz;
|
||
|
||
// ---- 3) 逐线逐道逐通道逐样本 → grid sum/count（重叠均值）----
|
||
std::vector<double> sum(cellCount, 0.0);
|
||
std::vector<std::uint16_t> cnt(cellCount, 0);
|
||
|
||
auto cellIdx = [&](int gi, int gj, int gk) -> std::size_t {
|
||
return (static_cast<std::size_t>(gk) * ny + gj) * nx + gi;
|
||
};
|
||
|
||
constexpr std::int64_t kChunk = 256; // 逐线分段读道，内存有界
|
||
double vmin = std::numeric_limits<double>::infinity();
|
||
double vmax = -std::numeric_limits<double>::infinity();
|
||
|
||
for (std::size_t i = 0; i < scales.size(); ++i) {
|
||
const LineScale& sc = scales[i];
|
||
const std::int64_t nt = sc.totalTraces;
|
||
if (nt < 1 || sc.trackM.size() < 2 || sc.chanX.empty()) continue;
|
||
const int nChan = static_cast<int>(sc.chanX.size());
|
||
|
||
// 各通道 .iprb 区间读（与 chanX 文件顺序对齐：iprb[c] ↔ chanX[c]）。
|
||
const auto& iprbPaths = lines[i].iprb;
|
||
if (static_cast<int>(iprbPaths.size()) != nChan) {
|
||
// 通道数与 .ord 有效通道不一致：按较小者对齐，避免越界。
|
||
}
|
||
const int useChan = std::min<int>(nChan, static_cast<int>(iprbPaths.size()));
|
||
|
||
for (std::int64_t t0 = 0; t0 < nt; t0 += kChunk) {
|
||
const std::int64_t t1 = std::min(nt, t0 + kChunk);
|
||
// 逐通道读该道段。
|
||
std::vector<geopro::io::gpr::BScan> bs(useChan);
|
||
for (int c = 0; c < useChan; ++c)
|
||
bs[c] = geopro::io::gpr::readIprbRange(iprbPaths[c], sc.header, t0, t1);
|
||
|
||
for (std::int64_t t = t0; t < t1; ++t) {
|
||
const double frac = nt > 1 ? static_cast<double>(t) / (nt - 1) : 0.0;
|
||
PosHeading ph = interpAlongTrack(sc.trackM, frac);
|
||
const double perpX = -ph.hy, perpY = ph.hx;
|
||
const std::int64_t lt = t - t0; // BScan 内局部道索引
|
||
|
||
for (int c = 0; c < useChan; ++c) {
|
||
const double oc = sc.chanX[c];
|
||
const double cx = ph.pos.x + oc * perpX;
|
||
const double cy = ph.pos.y + oc * perpY;
|
||
double rx, ry;
|
||
rotate(cx, cy, rx, ry);
|
||
const int gi = static_cast<int>(std::lround((rx - minX) / spec.cellXY));
|
||
const int gj = static_cast<int>(std::lround((ry - minY) / spec.cellXY));
|
||
if (gi < 0 || gi >= nx || gj < 0 || gj >= ny) continue;
|
||
|
||
for (int s = 0; s < samples; ++s) {
|
||
const double depth = geopro::io::gpr::depthOfSample(s, sc.header);
|
||
const int gk = static_cast<int>(std::lround(depth / spec.cellZ));
|
||
if (gk < 0 || gk >= nz) continue;
|
||
const double v = static_cast<double>(
|
||
bs[c].data[static_cast<std::size_t>(lt) * samples + s]);
|
||
const std::size_t idx = cellIdx(gi, gj, gk);
|
||
sum[idx] += v;
|
||
++cnt[idx];
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// ---- 4) 均值 + 求值域 ----
|
||
std::int64_t filled = 0;
|
||
for (std::size_t idx = 0; idx < cellCount; ++idx) {
|
||
if (cnt[idx] == 0) continue;
|
||
const double mean = sum[idx] / cnt[idx];
|
||
sum[idx] = mean; // 复用 sum 存均值
|
||
vmin = std::min(vmin, mean);
|
||
vmax = std::max(vmax, mean);
|
||
++filled;
|
||
}
|
||
if (!(vmin <= vmax)) {
|
||
vmin = 0.0;
|
||
vmax = 0.0;
|
||
}
|
||
|
||
// ---- 5) 量化 + 写 store + 金字塔 ----
|
||
Quant quant;
|
||
quant.scale = (vmax > vmin) ? (vmax - vmin) / 64000.0 : 1.0;
|
||
quant.offset = 0.5 * (vmin + vmax);
|
||
|
||
ScalarVolumeI16 vol(nx, ny, nz);
|
||
for (int gk = 0; gk < nz; ++gk)
|
||
for (int gj = 0; gj < ny; ++gj)
|
||
for (int gi = 0; gi < nx; ++gi) {
|
||
const std::size_t idx = cellIdx(gi, gj, gk);
|
||
vol.at(gi, gj, gk) = (cnt[idx] == 0)
|
||
? ScalarVolumeI16::kBlank
|
||
: quant.toQ(sum[idx]);
|
||
}
|
||
|
||
BuiltI16 built;
|
||
built.vol = std::move(vol);
|
||
built.quant = quant;
|
||
built.origin = {minX, minY, 0.0};
|
||
built.spacing = {spec.cellXY, spec.cellXY, spec.cellZ};
|
||
built.vminPhys = vmin;
|
||
built.vmaxPhys = vmax;
|
||
|
||
fs::create_directories(outStoreDir);
|
||
geopro::data::ChunkedVolumeStore::write(outStoreDir, built);
|
||
if (pyramidLevels > 0) {
|
||
geopro::data::ChunkedVolumeStore store(outStoreDir);
|
||
store.buildPyramid(pyramidLevels);
|
||
}
|
||
|
||
GeoBuildResult r;
|
||
r.nx = nx;
|
||
r.ny = ny;
|
||
r.nz = nz;
|
||
r.oxM = minX;
|
||
r.oyM = minY;
|
||
r.rotRad = rotRad;
|
||
r.filled = filled;
|
||
r.total = static_cast<std::int64_t>(cellCount);
|
||
return r;
|
||
}
|
||
|
||
} // namespace geopro::core
|