diff --git a/src/core/algo/GeoVolumeBuilder.cpp b/src/core/algo/GeoVolumeBuilder.cpp new file mode 100644 index 0000000..2e09313 --- /dev/null +++ b/src/core/algo/GeoVolumeBuilder.cpp @@ -0,0 +1,335 @@ +// GeoVolumeBuilder:按真实 RTK 几何把多条测线插值进统一路向网格。 +// +// 编排层(io_gpr + core 采样核 + store),故与 StreamingVolumeBuilder 一样编进 +// geopro_store(不污染纯 geopro_core)。命名空间保持 geopro::core(接口契约)。 +#include "core/algo/GeoVolumeBuilder.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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(f)), + std::istreambuf_iterator()); + return geopro::io::gpr::parseIprHeader(text); +} + +// 一条线的标尺(轨迹局部米 + 通道横偏 + 头)。 +struct LineScale { + std::vector trackM; // 该线轨迹(局部米,旋转前) + std::vector chanX; // 通道横偏(米,文件顺序) + geopro::io::gpr::IprHeader header; + std::int64_t totalTraces = 0; +}; + +// 全线总道数 = min 通道(fileBytes/(samples*2))(与 assembler 口径一致)。 +std::int64_t totalTracesOf(const std::vector& iprb, int samples) { + std::int64_t minTr = std::numeric_limits::max(); + const std::int64_t per = static_cast(samples) * 2; + if (per <= 0) throw std::runtime_error("samples<=0"); + for (const auto& p : iprb) { + const std::int64_t bytes = static_cast(fs::file_size(p)); + minTr = std::min(minTr, bytes / per); + } + return minTr; +} + +} // namespace + +// ---- PCA 主轴角(可测纯函数)---- +double principalAxisAngle(const std::vector& xs, + const std::vector& 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(n); + my /= static_cast(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(n); + syy /= static_cast(n); + sxy /= static_cast(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& 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 tracks(lines.size()); + double minLat = std::numeric_limits::infinity(); + double minLon = std::numeric_limits::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 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(ordF)), + std::istreambuf_iterator()); + 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 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(1, nt / 200); + for (std::int64_t t = 0; t < nt; t += stride) { + const double frac = nt > 1 ? static_cast(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::infinity(); + double minY = std::numeric_limits::infinity(); + double maxX = -std::numeric_limits::infinity(); + double maxY = -std::numeric_limits::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(1, nt / 200); + for (std::int64_t t = 0; t < nt; t += stride) { + const double frac = nt > 1 ? static_cast(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(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(nx) * ny * nz; + + // ---- 3) 逐线逐道逐通道逐样本 → grid sum/count(重叠均值)---- + std::vector sum(cellCount, 0.0); + std::vector cnt(cellCount, 0); + + auto cellIdx = [&](int gi, int gj, int gk) -> std::size_t { + return (static_cast(gk) * ny + gj) * nx + gi; + }; + + constexpr std::int64_t kChunk = 256; // 逐线分段读道,内存有界 + double vmin = std::numeric_limits::infinity(); + double vmax = -std::numeric_limits::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(sc.chanX.size()); + + // 各通道 .iprb 区间读(与 chanX 文件顺序对齐:iprb[c] ↔ chanX[c])。 + const auto& iprbPaths = lines[i].iprb; + if (static_cast(iprbPaths.size()) != nChan) { + // 通道数与 .ord 有效通道不一致:按较小者对齐,避免越界。 + } + const int useChan = std::min(nChan, static_cast(iprbPaths.size())); + + for (std::int64_t t0 = 0; t0 < nt; t0 += kChunk) { + const std::int64_t t1 = std::min(nt, t0 + kChunk); + // 逐通道读该道段。 + std::vector 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(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(std::lround((rx - minX) / spec.cellXY)); + const int gj = static_cast(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(std::lround(depth / spec.cellZ)); + if (gk < 0 || gk >= nz) continue; + const double v = static_cast( + bs[c].data[static_cast(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(cellCount); + return r; +} + +} // namespace geopro::core diff --git a/src/core/algo/GeoVolumeBuilder.hpp b/src/core/algo/GeoVolumeBuilder.hpp new file mode 100644 index 0000000..10babbf --- /dev/null +++ b/src/core/algo/GeoVolumeBuilder.hpp @@ -0,0 +1,57 @@ +#ifndef GEOPRO_CORE_ALGO_GEOVOLUMEBUILDER_HPP +#define GEOPRO_CORE_ALGO_GEOVOLUMEBUILDER_HPP + +#include +#include +#include + +namespace geopro::core { + +// 统一路向网格规格:横向(X'沿路 / Y'横路)cell = cellXY,深度 cell = cellZ。 +struct GeoGridSpec { + double cellXY = 0.5; // 米 + double cellZ = 0.1; // 米 +}; + +// 一组二维路向点的主成分朝向(协方差最大特征向量),用于把路向对齐到网格长轴 X'。 +// 返回主轴方向角 rotRad(弧度,相对东向 +X),调用方旋转 -rotRad 使路向 → X'。 +// 点少于 2 / 退化(零方差)→ 返回 0。 +double principalAxisAngle(const std::vector& xs, + const std::vector& ys); + +// 真实数据建体的产物指标(维度/原点/旋转角/填充率)。 +struct GeoBuildResult { + int nx = 0, ny = 0, nz = 0; + double oxM = 0, oyM = 0; // 旋转后网格原点(路向坐标系,米) + double rotRad = 0; // 路向主轴角(弧度,East 起算) + std::int64_t filled = 0; // 非空 cell 数 + std::int64_t total = 0; // 总 cell 数 +}; + +// 一条线的输入文件束(各通道 .iprb 已按通道号排序 + 该线 .ord + .gps)。 +// .iprh 与每个 .iprb 同目录同名(扩展名替换为 .iprh),由读取层自行推导。 +struct GeoLineInput { + std::vector iprb; // 该线 14 通道 .iprb(按通道号升序) + std::string ord; // 该线 .ord(通道横偏) + std::string gps; // 该线 .gps(RTK 轨迹) +}; + +// 把多条测线按各自真实 RTK 几何插值进一个统一路向网格,重叠 cell 取均值, +// 量化为 int16 写 ChunkedVolumeStore + 金字塔。 +// +// 步骤: +// 1) 各线 .gps → 经纬,全体共用原点(最小 lat/lon) → 局部米; +// 2) PCA 求路向主轴角 rotRad,旋转使路向 = X'、横路 = Y';Z = 深度; +// 3) 每线逐道(里程均匀分布插值定位)逐通道(横偏垂直航向摆放)逐样本 → 网格 cell, +// 累加 sum + count; +// 4) 空 cell = blank;非空 = sum/count 均值; +// 5) 扫值域定 Quant(offset=中点),量化写 store + buildPyramid。 +// +// 内存有界:逐线流式读道(readIprbRange 区间),只持 grid 的 sum/count 累加缓冲。 +GeoBuildResult buildGeoVolume(const std::vector& lines, + const GeoGridSpec& spec, + const std::string& outStoreDir, int pyramidLevels); + +} // namespace geopro::core + +#endif // GEOPRO_CORE_ALGO_GEOVOLUMEBUILDER_HPP diff --git a/src/data/store/CMakeLists.txt b/src/data/store/CMakeLists.txt index ab139b5..9c16f22 100644 --- a/src/data/store/CMakeLists.txt +++ b/src/data/store/CMakeLists.txt @@ -6,7 +6,10 @@ find_package(Qt6 COMPONENTS Core REQUIRED) add_library(geopro_store STATIC ChunkedVolumeStore.cpp # 流式建体(B4):编排 io_gpr+core+store,沿 X 分 slab 建 level0 体(命名空间 geopro::data)。 - ${CMAKE_CURRENT_SOURCE_DIR}/../StreamingVolumeBuilder.cpp) + ${CMAKE_CURRENT_SOURCE_DIR}/../StreamingVolumeBuilder.cpp + # build-geo(G1):按真实 RTK 几何把多线插值进统一路向网格(命名空间 geopro::core, + # 文件物理位于 src/core/algo/,但编排 io_gpr+core+store 故编进 store,不污染纯 core)。 + ${CMAKE_SOURCE_DIR}/src/core/algo/GeoVolumeBuilder.cpp) # include 根 = src/,使 #include "data/store/..." 与 "core/algo/..." 可解析 # (geopro_tests 链 geopro_store 后透传)。 diff --git a/src/io/gpr/CMakeLists.txt b/src/io/gpr/CMakeLists.txt index a66381b..e670c63 100644 --- a/src/io/gpr/CMakeLists.txt +++ b/src/io/gpr/CMakeLists.txt @@ -1,5 +1,6 @@ add_library(geopro_io_gpr STATIC - IprHeader.cpp IprbReader.cpp GprGeometry.cpp GprSurveyAssembler.cpp) + IprHeader.cpp IprbReader.cpp GprGeometry.cpp GprSurveyAssembler.cpp + GpsTrack.cpp) target_include_directories(geopro_io_gpr PUBLIC ${CMAKE_SOURCE_DIR}/src) target_compile_features(geopro_io_gpr PUBLIC cxx_std_17) # GprSurveyAssembler 返回 geopro::core::GprSurvey(头文件内联,仅需 include 解析)。 diff --git a/src/io/gpr/GpsTrack.cpp b/src/io/gpr/GpsTrack.cpp new file mode 100644 index 0000000..fd00805 --- /dev/null +++ b/src/io/gpr/GpsTrack.cpp @@ -0,0 +1,109 @@ +#include "io/gpr/GpsTrack.hpp" + +#include +#include +#include +#include + +namespace geopro::io::gpr { + +namespace { +constexpr double kMetersPerDegLat = 111320.0; +constexpr double kPi = 3.14159265358979323846; + +// 严格把整列解析为 double;失败返回 false(不抛,供逐行容错跳过)。 +bool parseDouble(const std::string& s, double& out) { + try { + std::size_t used = 0; + out = std::stod(s, &used); + return used == s.size(); + } catch (...) { + return false; + } +} +} // namespace + +GpsTrack parseGps(const std::string& path) { + std::ifstream f(path); + if (!f) throw std::runtime_error("parseGps: 打不开 " + path); + + GpsTrack track; + std::string line; + while (std::getline(f, line)) { + // 去掉可能的 \r(CRLF)。 + if (!line.empty() && line.back() == '\r') line.pop_back(); + + std::istringstream tok(line); + std::vector cols; + std::string c; + while (tok >> c) cols.push_back(c); + if (cols.size() < 7) continue; // 列数不足 + + GpsPt p; + if (!parseDouble(cols[2], p.lat)) continue; // 纬 + if (!parseDouble(cols[4], p.lon)) continue; // 经 + if (!parseDouble(cols[6], p.elev)) continue; // 高 + track.pts.push_back(p); + } + return track; +} + +XY lonLatToLocalM(double lat, double lon, double lat0, double lon0) { + XY xy; + xy.x = (lon - lon0) * kMetersPerDegLat * std::cos(lat0 * kPi / 180.0); + xy.y = (lat - lat0) * kMetersPerDegLat; + return xy; +} + +PosHeading interpAlongTrack(const std::vector& trackM, double frac) { + PosHeading r; + if (trackM.empty()) return r; + if (trackM.size() == 1) { + r.pos = trackM.front(); + return r; + } + + // 各段长度 + 累积里程。 + const std::size_t n = trackM.size(); + std::vector cum(n, 0.0); + for (std::size_t i = 1; i < n; ++i) { + const double dx = trackM[i].x - trackM[i - 1].x; + const double dy = trackM[i].y - trackM[i - 1].y; + cum[i] = cum[i - 1] + std::hypot(dx, dy); + } + const double total = cum.back(); + + // 退化(零长轨迹):返回起点,航向取首段方向或 (1,0)。 + if (total <= 0.0) { + r.pos = trackM.front(); + return r; + } + + if (frac <= 0.0) frac = 0.0; + if (frac >= 1.0) frac = 1.0; + const double target = frac * total; + + // 找包含 target 里程的段 [i-1, i]。 + std::size_t seg = 1; + while (seg < n && cum[seg] < target) ++seg; + if (seg >= n) seg = n - 1; + + const double segLen = cum[seg] - cum[seg - 1]; + const double localFrac = segLen > 0.0 ? (target - cum[seg - 1]) / segLen : 0.0; + + const XY& a = trackM[seg - 1]; + const XY& b = trackM[seg]; + r.pos.x = a.x + (b.x - a.x) * localFrac; + r.pos.y = a.y + (b.y - a.y) * localFrac; + + double dx = b.x - a.x; + double dy = b.y - a.y; + const double len = std::hypot(dx, dy); + if (len > 0.0) { + r.hx = dx / len; + r.hy = dy / len; + } + return r; +} + +} // namespace geopro::io::gpr diff --git a/src/io/gpr/GpsTrack.hpp b/src/io/gpr/GpsTrack.hpp new file mode 100644 index 0000000..5996851 --- /dev/null +++ b/src/io/gpr/GpsTrack.hpp @@ -0,0 +1,46 @@ +#ifndef GEOPRO_IO_GPR_GPSTRACK_HPP +#define GEOPRO_IO_GPR_GPSTRACK_HPP + +#include +#include + +namespace geopro::io::gpr { + +// 一个 RTK 轨迹点(度,度,米)。 +struct GpsPt { + double lat = 0; // 纬度(度) + double lon = 0; // 经度(度) + double elev = 0; // 高程(米) +}; + +// 一条测线的 RTK 轨迹。 +struct GpsTrack { + std::vector pts; +}; + +// 解析 .gps 文本。每行 tab/空白分隔: +// 日期 \t 时间 \t 纬 \t N \t 经 \t E \t 高 \t M \t 卫星 +// 取 纬(col2)/经(col4)/高(col6)。列数不足或非数字行跳过。文件打不开抛 std::runtime_error。 +GpsTrack parseGps(const std::string& path); + +// 局部米坐标(等距投影,绕给定原点 lat0/lon0)。 +// x_east = (lon-lon0)*111320*cos(lat0°) +// y_north = (lat-lat0)*111320 +struct XY { + double x = 0; // 东向米 + double y = 0; // 北向米 +}; +XY lonLatToLocalM(double lat, double lon, double lat0, double lon0); + +// 沿轨迹按里程插值:frac∈[0,1] 返回该里程分数处的局部米坐标 + 该处航向(单位向量)。 +// frac<=0 → 起点;frac>=1 → 终点;空/单点轨迹航向 (1,0)。 +struct PosHeading { + XY pos; + double hx = 1; // 航向单位向量 x + double hy = 0; // 航向单位向量 y +}; +PosHeading interpAlongTrack(const std::vector& trackM, double frac); + +} // namespace geopro::io::gpr + +#endif // GEOPRO_IO_GPR_GPSTRACK_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e50fdf4..66a2a60 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -62,6 +62,9 @@ target_sources(geopro_tests PRIVATE data/store/test_pyramid.cpp) target_sources(geopro_tests PRIVATE data/store/test_streaming_write.cpp) # data 层:StreamingVolumeBuilder(流式建体 B4;与非流式 buildGprVolume+write 逐 brick+meta 对拍)。 target_sources(geopro_tests PRIVATE data/test_streaming_builder.cpp) +# core/algo:GeoVolumeBuilder(G1 build-geo:PCA 路向旋转 + 多线统一网格重叠均值; +# 编排 io_gpr+core+store,符号编在 geopro_store,故归此节)。 +target_sources(geopro_tests PRIVATE core/test_geo_volume_builder.cpp) target_link_libraries(geopro_tests PRIVATE geopro_store) # net 层:RSA 加密器。测试需直接用 OpenSSL 生成/解密密钥,故显式 find_package @@ -207,6 +210,8 @@ target_sources(geopro_tests PRIVATE io/gpr/test_iprb_reader.cpp) target_sources(geopro_tests PRIVATE io/gpr/test_gpr_geometry.cpp) # GprSurveyAssembler:若干通道 .iprb + .ord -> GprSurvey(samples 校验/traces 对齐/Y 升序重排)。 target_sources(geopro_tests PRIVATE io/gpr/test_gpr_survey_assembler.cpp) +# GpsTrack:.gps 解析 + 经纬→局部米 + 沿轨迹里程插值/航向(G1 build-geo 基础,纯 C++17)。 +target_sources(geopro_tests PRIVATE io/gpr/test_gps_track.cpp) target_link_libraries(geopro_tests PRIVATE geopro_io_gpr) add_subdirectory(spike) # spike S3: banded contour 渲染验证 diff --git a/tests/core/test_geo_volume_builder.cpp b/tests/core/test_geo_volume_builder.cpp new file mode 100644 index 0000000..1c3343a --- /dev/null +++ b/tests/core/test_geo_volume_builder.cpp @@ -0,0 +1,154 @@ +#include "core/algo/GeoVolumeBuilder.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "data/store/ChunkedVolumeStore.hpp" + +namespace fs = std::filesystem; +using namespace geopro::core; + +// ---- principalAxisAngle 纯函数 ---- +TEST(GeoVolumeBuilder, PcaFindsRoadDirection) { + // 点沿 +Y(北向)排布 → 主轴角约 ±90°(±pi/2)。 + std::vector xs, ys; + for (int i = 0; i < 20; ++i) { + xs.push_back(0.01 * i); // 微小横向噪声 + ys.push_back(static_cast(i)); + } + const double ang = principalAxisAngle(xs, ys); + EXPECT_NEAR(std::abs(ang), 3.14159265358979323846 / 2.0, 0.05); +} + +TEST(GeoVolumeBuilder, PcaAlongEastIsZero) { + std::vector xs, ys; + for (int i = 0; i < 20; ++i) { + xs.push_back(static_cast(i)); + ys.push_back(0.01 * i); + } + EXPECT_NEAR(principalAxisAngle(xs, ys), 0.0, 0.05); +} + +TEST(GeoVolumeBuilder, PcaDegenerateReturnsZero) { + EXPECT_DOUBLE_EQ(principalAxisAngle({}, {}), 0.0); + EXPECT_DOUBLE_EQ(principalAxisAngle({1.0}, {1.0}), 0.0); +} + +// ---- 合成小线建体 ---- +namespace { + +// 写一通道 .iprb + .iprh:值恒为 fixedVal(便于校验重叠均值)。 +void writeChannel(const fs::path& iprb, int samples, int traces, + std::int16_t fixedVal) { + fs::path iprh = fs::path(iprb).replace_extension(".iprh"); + std::ofstream h(iprh); + h << "SAMPLES: " << samples << "\n"; + h << "LAST TRACE: " << (traces - 1) << "\n"; + h << "CHANNELS: 2\n"; + h << "TIMEWINDOW: 100.0\n"; + h << "SOIL VELOCITY: 100.0\n"; // m/µs → 1e8 m/s + h << "DISTANCE INTERVAL: 0.05\n"; + h.close(); + + std::ofstream b(iprb, std::ios::binary); + for (int t = 0; t < traces; ++t) + for (int s = 0; s < samples; ++s) + b.write(reinterpret_cast(&fixedVal), sizeof(fixedVal)); +} + +// 写一条南北直线的 .gps(lat 从 lat0 递增到 lat1,lon 固定)。 +void writeGps(const fs::path& path, double lat0, double lat1, double lon, + int pts) { + std::ofstream f(path); + for (int i = 0; i < pts; ++i) { + const double frac = pts > 1 ? static_cast(i) / (pts - 1) : 0.0; + const double lat = lat0 + (lat1 - lat0) * frac; + f << "2023-06-03\t00:00:00:000\t" << std::fixed + << std::setprecision(10) << lat << "\tN\t" << lon << "\tE\t9.0\tM\t4\n"; + } +} + +// 写两通道 .ord(横偏 -0.5 / +0.5,末列=1 有效)。 +void writeOrd(const fs::path& path) { + std::ofstream f(path); + f << "0 -0.500000 -1.5 1\n"; + f << "1 0.500000 -1.5 1\n"; +} + +GeoLineInput makeLine(const fs::path& dir, const std::string& tag, double lat0, + double lat1, double lon, int traces, std::int16_t val) { + const int samples = 8; + writeChannel(dir / (tag + "_A01.iprb"), samples, traces, val); + writeChannel(dir / (tag + "_A02.iprb"), samples, traces, val); + writeOrd(dir / (tag + ".ord")); + writeGps(dir / (tag + ".gps"), lat0, lat1, lon, traces); + GeoLineInput in; + in.iprb = {(dir / (tag + "_A01.iprb")).string(), + (dir / (tag + "_A02.iprb")).string()}; + in.ord = (dir / (tag + ".ord")).string(); + in.gps = (dir / (tag + ".gps")).string(); + return in; +} + +} // namespace + +TEST(GeoVolumeBuilder, BuildsSyntheticLinesOverlapAveraged) { + const fs::path tmp = fs::temp_directory_path() / "geopro_geovol_test"; + std::error_code ec; + fs::remove_all(tmp, ec); + fs::create_directories(tmp); + + // 两条同位置南北线(lat 30.200→30.201,~111m),值不同 → 重叠 cell 取均值。 + const double lat0 = 30.200, lat1 = 30.201, lon = 120.244; + std::vector lines = { + makeLine(tmp, "synA_001", lat0, lat1, lon, /*traces=*/40, /*val=*/100), + makeLine(tmp, "synB_002", lat0, lat1, lon, /*traces=*/40, /*val=*/300), + }; + + const std::string store = (tmp / "store").string(); + GeoGridSpec spec{/*cellXY=*/0.5, /*cellZ=*/0.1}; + GeoBuildResult r = buildGeoVolume(lines, spec, store, /*pyramidLevels=*/1); + + // 维度合理:约 111m 长 → nx 在数百量级;横路窄(~1m 阵列) → ny 小;nz>1。 + EXPECT_GT(r.nx, 50); + EXPECT_GE(r.ny, 1); + EXPECT_GT(r.nz, 1); + EXPECT_GT(r.filled, 0); + EXPECT_LE(r.filled, r.total); + + // store 可读,维度一致。 + geopro::data::ChunkedVolumeStore s(store); + EXPECT_EQ(s.meta().nx, r.nx); + EXPECT_EQ(s.meta().ny, r.ny); + EXPECT_EQ(s.meta().nz, r.nz); + EXPECT_EQ(s.levels(), 2); // level0 + 1 + + // 重叠均值:两线值 100/300,命中同 cell → 均值 200。扫所有块找非 blank 体素, + // 其物理值应接近 200(量化误差内)。 + const auto& m = s.meta(); + bool foundNonBlank = false; + double sampleVal = 0.0; + for (int bz = 0; bz < s.bricksZ() && !foundNonBlank; ++bz) + for (int by = 0; by < s.bricksY() && !foundNonBlank; ++by) + for (int bx = 0; bx < s.bricksX() && !foundNonBlank; ++bx) { + auto vox = s.readBrick(bx, by, bz); + for (std::int16_t q : vox) { + if (q != geopro::core::ScalarVolumeI16::kBlank) { + sampleVal = m.quant.toPhys(q); + foundNonBlank = true; + break; + } + } + } + ASSERT_TRUE(foundNonBlank); + EXPECT_NEAR(sampleVal, 200.0, 2.0); + + fs::remove_all(tmp, ec); +} diff --git a/tests/io/gpr/test_gps_track.cpp b/tests/io/gpr/test_gps_track.cpp new file mode 100644 index 0000000..410417d --- /dev/null +++ b/tests/io/gpr/test_gps_track.cpp @@ -0,0 +1,92 @@ +#include "io/gpr/GpsTrack.hpp" + +#include + +#include +#include +#include +#include + +using namespace geopro::io::gpr; + +// 经纬→米:纬 1° ≈ 111320 m;经按 cos(lat0) 缩。 +TEST(GpsTrack, LonLatToLocalMatchesPhysics) { + const double lat0 = 30.204, lon0 = 120.244; + // 原点处归零。 + XY o = lonLatToLocalM(lat0, lon0, lat0, lon0); + EXPECT_NEAR(o.x, 0.0, 1e-9); + EXPECT_NEAR(o.y, 0.0, 1e-9); + + // 北移 0.001° 纬 → y ≈ 111.32 m,x≈0。 + XY north = lonLatToLocalM(lat0 + 0.001, lon0, lat0, lon0); + EXPECT_NEAR(north.y, 111.32, 0.01); + EXPECT_NEAR(north.x, 0.0, 1e-6); + + // 东移 0.001° 经 → x ≈ 111.32*cos(30.204°) ≈ 96.2 m。 + XY east = lonLatToLocalM(lat0, lon0 + 0.001, lat0, lon0); + const double expX = 111.32 * std::cos(lat0 * 3.14159265358979323846 / 180.0); + EXPECT_NEAR(east.x, expX, 0.01); + EXPECT_NEAR(east.y, 0.0, 1e-6); +} + +// 直线轨迹:frac=0/0.5/1 → 位置/航向正确。 +TEST(GpsTrack, InterpStraightLine) { + std::vector tr = {{0, 0}, {10, 0}, {20, 0}}; // 沿 +X 直线,长 20 + auto p0 = interpAlongTrack(tr, 0.0); + EXPECT_NEAR(p0.pos.x, 0.0, 1e-9); + EXPECT_NEAR(p0.hx, 1.0, 1e-9); + EXPECT_NEAR(p0.hy, 0.0, 1e-9); + + auto pm = interpAlongTrack(tr, 0.5); + EXPECT_NEAR(pm.pos.x, 10.0, 1e-9); + EXPECT_NEAR(pm.pos.y, 0.0, 1e-9); + + auto p1 = interpAlongTrack(tr, 1.0); + EXPECT_NEAR(p1.pos.x, 20.0, 1e-9); + EXPECT_NEAR(p1.hx, 1.0, 1e-9); +} + +// 折线(L 形):里程一半落在拐角后,航向应是第二段方向。 +TEST(GpsTrack, InterpPolyline) { + std::vector tr = {{0, 0}, {10, 0}, {10, 10}}; // 先 +X 10,再 +Y 10,总长 20 + auto pm = interpAlongTrack(tr, 0.5); // 里程 10 = 正好拐角 + EXPECT_NEAR(pm.pos.x, 10.0, 1e-9); + EXPECT_NEAR(pm.pos.y, 0.0, 1e-9); + + auto p75 = interpAlongTrack(tr, 0.75); // 里程 15 → 第二段中点 (10,5) + EXPECT_NEAR(p75.pos.x, 10.0, 1e-9); + EXPECT_NEAR(p75.pos.y, 5.0, 1e-9); + EXPECT_NEAR(p75.hx, 0.0, 1e-9); // 第二段沿 +Y + EXPECT_NEAR(p75.hy, 1.0, 1e-9); +} + +// 通道横向定位:航向 (hx,hy),垂直 = (-hy,hx);通道位置 = trace 位置 + o_c×垂直。 +TEST(GpsTrack, ChannelLateralPlacement) { + auto ph = interpAlongTrack({{0, 0}, {0, 10}}, 0.5); // 沿 +Y,航向 (0,1) + EXPECT_NEAR(ph.hx, 0.0, 1e-9); + EXPECT_NEAR(ph.hy, 1.0, 1e-9); + // 垂直(航向) = (-hy, hx) = (-1, 0)。偏移 +0.5 → 通道在 (-0.5, 5)。 + const double perpX = -ph.hy, perpY = ph.hx; + const double oc = 0.5; + const double cx = ph.pos.x + oc * perpX; + const double cy = ph.pos.y + oc * perpY; + EXPECT_NEAR(cx, -0.5, 1e-9); + EXPECT_NEAR(cy, 5.0, 1e-9); +} + +// 解析真实格式的 .gps(tab 分隔,含 N/E/M 标记列)。 +TEST(GpsTrack, ParsesGpsFile) { + auto tmp = std::filesystem::temp_directory_path() / "geopro_gps_test.gps"; + { + std::ofstream f(tmp); + f << "2023-06-03\t14:42:23:000\t30.21402519\tN\t120.24466077\tE\t9.390\tM\t4\r\n"; + f << "2023-06-03\t14:42:23:203\t30.21402388\tN\t120.24466074\tE\t9.386\tM\t4\r\n"; + f << "garbage line should be skipped\n"; + } + GpsTrack t = parseGps(tmp.string()); + std::filesystem::remove(tmp); + ASSERT_EQ(t.pts.size(), 2u); + EXPECT_NEAR(t.pts[0].lat, 30.21402519, 1e-8); + EXPECT_NEAR(t.pts[0].lon, 120.24466077, 1e-8); + EXPECT_NEAR(t.pts[0].elev, 9.390, 1e-6); +} diff --git a/tools/gpr_poc/main.cpp b/tools/gpr_poc/main.cpp index 3c83d1e..41d9da6 100644 --- a/tools/gpr_poc/main.cpp +++ b/tools/gpr_poc/main.cpp @@ -26,6 +26,7 @@ #include "Probe.hpp" +#include "core/algo/GeoVolumeBuilder.hpp" #include "core/algo/GprVolumeBuilder.hpp" #include "core/algo/IInterpolator.hpp" #include "core/model/GprSurvey.hpp" @@ -681,6 +682,117 @@ int cmdBuildStream(int argc, char** argv) { return 0; } +// ============================================================================ +// build-geo:按真实 RTK 几何把多线插值进统一路向网格(Task G1) +// ============================================================================ +// +// 消除 build-stream 顺序拼接的退化扁带:各线 .gps RTK 轨迹 → 经纬 → 局部米 → +// PCA 路向旋转 → 道按里程均匀分布定位、14 通道横偏垂直航向摆放 → 全部插进统一 +// 路向网格(≈4472×43×81)重叠取均值 → 量化写 ChunkedVolumeStore + 金字塔。 +int cmdBuildGeo(int argc, char** argv) { + const Args a = parseArgs(argc, argv, 2); + if (a.positional.empty()) { + std::cerr << "用法: gpr_poc build-geo [--cellXY 0.5] [--cellZ 0.1] " + "[--out ] [--levels 4] [--maxLines N]\n"; + return 2; + } + const std::string dir = a.positional[0]; + const double cellXY = std::stod(a.get("cellXY", "0.5")); + const double cellZ = std::stod(a.get("cellZ", "0.1")); + const int levels = std::stoi(a.get("levels", "4")); + const int maxLines = std::stoi(a.get("maxLines", "0")); + const std::string out = + a.get("out", (fs::temp_directory_path() / "gpr_store_geo").string()); + + std::cout << "[build-geo] dir=" << dir << " cellXY=" << cellXY + << " cellZ=" << cellZ << " levels=" << levels << " out=" << out + << "\n"; + + // 发现线号 → 各线 iprb/ord(复用 discoverLine)+ .gps(按线号匹配)。 + std::vector lineNos = discoverLines(dir); + if (maxLines > 0 && static_cast(lineNos.size()) > maxLines) + lineNos.resize(maxLines); + std::cout << "[build-geo] 发现测线数=" << lineNos.size() << "\n"; + if (lineNos.empty()) { + std::cerr << "[build-geo] 错误: 未发现任何 .ord 测线\n"; + return 1; + } + + std::vector lines; + for (const std::string& ln : lineNos) { + const LineFiles lf = discoverLine(dir, ln); + if (lf.iprb.empty() || lf.ord.empty()) { + std::cerr << "[build-geo] 警告: 线 " << ln << " 缺 iprb/ord,跳过\n"; + continue; + } + // .gps:匹配 "*_.gps"。 + std::string gps; + for (const auto& e : fs::directory_iterator(dir)) { + if (!e.is_regular_file()) continue; + if (e.path().extension().string() != ".gps") continue; + const std::string stem = e.path().stem().string(); + const std::size_t us = stem.find_last_of('_'); + if (us != std::string::npos && stem.substr(us + 1) == ln) { + gps = e.path().string(); + break; + } + } + if (gps.empty()) { + std::cerr << "[build-geo] 警告: 线 " << ln << " 缺 .gps,跳过\n"; + continue; + } + geopro::core::GeoLineInput in; + in.iprb = lf.iprb; + in.ord = lf.ord; + in.gps = gps; + lines.push_back(std::move(in)); + } + if (lines.empty()) { + std::cerr << "[build-geo] 错误: 无可用测线\n"; + return 1; + } + std::cout << "[build-geo] 可用测线数=" << lines.size() << "\n"; + + fs::create_directories(out); + Stopwatch sw; + const geopro::core::GeoGridSpec spec{cellXY, cellZ}; + const geopro::core::GeoBuildResult r = + geopro::core::buildGeoVolume(lines, spec, out, levels); + const double buildMs = sw.elapsedMs(); + + const double fillRate = + r.total > 0 ? static_cast(r.filled) / r.total : 0.0; + const double rotDeg = r.rotRad * 180.0 / 3.14159265358979323846; + const std::int64_t dataBytes = storeDataBytes(out); + const double peak = Probe::peakMemMB(); + + std::cout << "\n=== build-geo 指标(真实 RTK 几何统一路向体)===\n"; + std::cout << "测线数 : " << lines.size() << "\n"; + std::cout << "体维度 : " << r.nx << " x " << r.ny << " x " << r.nz + << "\n"; + std::cout << "总 cell : " << r.total << "\n"; + std::cout << "非空 cell : " << r.filled << "\n"; + std::cout << "填充率 : " << (fillRate * 100.0) << " %\n"; + std::cout << "路向旋转角 : " << rotDeg << " ° (" << r.rotRad << " rad)\n"; + std::cout << "网格原点(m) : (" << r.oxM << ", " << r.oyM << ")\n"; + std::cout << "data.bin(B) : " << dataBytes << " (" + << dataBytes / (1024.0 * 1024.0) << " MB)\n"; + std::cout << "建体总耗时(ms) : " << buildMs << "\n"; + std::cout << "峰值内存(MB) : " << peak << "\n"; + + writeMetricLine( + "build-geo,lines=" + std::to_string(lines.size()) + + ",cellXY=" + std::to_string(cellXY) + ",cellZ=" + std::to_string(cellZ) + + ",nx=" + std::to_string(r.nx) + ",ny=" + std::to_string(r.ny) + + ",nz=" + std::to_string(r.nz) + ",total=" + std::to_string(r.total) + + ",filled=" + std::to_string(r.filled) + + ",fillRate=" + std::to_string(fillRate) + + ",rotDeg=" + std::to_string(rotDeg) + + ",buildMs=" + std::to_string(buildMs) + + ",peakMB=" + std::to_string(peak)); + return 0; +} + int cmdLoad(int argc, char** argv) { const Args a = parseArgs(argc, argv, 2); if (a.positional.empty()) { @@ -3822,6 +3934,8 @@ void usage() { " gpr_poc build-stream [--cellXY 0.05] [--cellZ 0.05] " "[--out ] [--levels 3] [--sliceXBricks 8] " "[--maxLines N]\n" + " gpr_poc build-geo [--cellXY 0.5] [--cellZ 0.1] " + "[--out ] [--levels 4] [--maxLines N]\n" " gpr_poc load \n" " gpr_poc selftest\n" " gpr_poc offscreen-smoke\n" @@ -3855,6 +3969,7 @@ int main(int argc, char** argv) { try { if (cmd == "build") return cmdBuild(argc, argv); if (cmd == "build-stream") return cmdBuildStream(argc, argv); + if (cmd == "build-geo") return cmdBuildGeo(argc, argv); if (cmd == "load") return cmdLoad(argc, argv); if (cmd == "selftest") return cmdSelftest(); if (cmd == "offscreen-smoke") return cmdOffscreenSmoke();