From f713c1536632ebc21c96d9466dd1fa3c235acf01 Mon Sep 17 00:00:00 2001 From: gaozheng Date: Thu, 25 Jun 2026 10:41:02 +0800 Subject: [PATCH] =?UTF-8?q?feat(gpr3dv):=20=E7=A7=BB=E6=A4=8D=E7=B2=BE?= =?UTF-8?q?=E7=A1=AE=E5=9D=90=E6=A0=87/=E8=BD=A8=E8=BF=B9/=E4=B8=96?= =?UTF-8?q?=E7=95=8C=E7=BD=91=E6=A0=BC(CGCS2000)+=E6=B5=8B=E7=BB=98?= =?UTF-8?q?=E7=BA=A7=E9=80=90=E7=BA=BF=E4=B8=96=E7=95=8C=E5=AF=B9=E9=BD=90?= =?UTF-8?q?=E5=BB=BA=E4=BD=93?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 复制 CoordinateTransform/TrajectoryCalculator/CScanGridder/PosParser(逐字节一致)进 external/gpr3dviewer;新增 Gpr3dvSurveyVolumeBridge 按 CGCS2000+逐道GPS轨迹建世界对齐体; gpr_poc 加 build-survey-all/view-survey-all(各体自带世界origin,精确就位跟GPS弯)。 --- external/gpr3dviewer/CMakeLists.txt | 9 + external/gpr3dviewer/CScanGridder.cpp | 336 +++++++++++ external/gpr3dviewer/CScanGridder.h | 48 ++ external/gpr3dviewer/CoordinateTransform.cpp | 219 +++++++ external/gpr3dviewer/CoordinateTransform.h | 48 ++ external/gpr3dviewer/PosParser.cpp | 73 +++ external/gpr3dviewer/PosParser.h | 17 + external/gpr3dviewer/TrajectoryCalculator.cpp | 242 ++++++++ external/gpr3dviewer/TrajectoryCalculator.h | 49 ++ src/io/gpr/CMakeLists.txt | 12 +- src/io/gpr/Gpr3dvSurveyVolumeBridge.cpp | 270 +++++++++ src/io/gpr/Gpr3dvSurveyVolumeBridge.hpp | 71 +++ tools/gpr_poc/main.cpp | 534 ++++++++++++++++++ 13 files changed, 1925 insertions(+), 3 deletions(-) create mode 100644 external/gpr3dviewer/CScanGridder.cpp create mode 100644 external/gpr3dviewer/CScanGridder.h create mode 100644 external/gpr3dviewer/CoordinateTransform.cpp create mode 100644 external/gpr3dviewer/CoordinateTransform.h create mode 100644 external/gpr3dviewer/PosParser.cpp create mode 100644 external/gpr3dviewer/PosParser.h create mode 100644 external/gpr3dviewer/TrajectoryCalculator.cpp create mode 100644 external/gpr3dviewer/TrajectoryCalculator.h create mode 100644 src/io/gpr/Gpr3dvSurveyVolumeBridge.cpp create mode 100644 src/io/gpr/Gpr3dvSurveyVolumeBridge.hpp diff --git a/external/gpr3dviewer/CMakeLists.txt b/external/gpr3dviewer/CMakeLists.txt index 140e85d..4066343 100644 --- a/external/gpr3dviewer/CMakeLists.txt +++ b/external/gpr3dviewer/CMakeLists.txt @@ -23,6 +23,15 @@ add_library(geopro_gpr3dv STATIC Rd3Parser.cpp RadarProcessor.cpp PerformanceLogger.cpp + # P8 测绘级精确坐标：原样拷自 3DGPRViewer（算法零改动，diff 逐字节一致）。 + # CoordinateTransform —— CGCS2000 高斯-克吕格 3°带正反算（零 Qt）。 + # TrajectoryCalculator —— RTK GPS + 天线几何 → 每通道每道 CGCS 世界坐标。 + # CScanGridder —— 世界坐标逐深度 IDW 网格化 C-scan（堆成世界对齐体）。 + # PosParser —— .pos / center.ccc GPS 解析。 + CoordinateTransform.cpp + TrajectoryCalculator.cpp + CScanGridder.cpp + PosParser.cpp third_party/kissfft/kiss_fft.c third_party/kissfft/kiss_fftr.c ) diff --git a/external/gpr3dviewer/CScanGridder.cpp b/external/gpr3dviewer/CScanGridder.cpp new file mode 100644 index 0000000..cd4cd8d --- /dev/null +++ b/external/gpr3dviewer/CScanGridder.cpp @@ -0,0 +1,336 @@ +#include "CScanGridder.h" + +#include +#include +#include +#include +#include +#include + +namespace { + +struct CScanSamplePoint { + double x = 0.0; + double y = 0.0; + float amplitude = 0.0f; + int line = -1; + int channel = -1; + int trace = -1; +}; + +struct BinKey { + int x = 0; + int y = 0; + + bool operator==(const BinKey &other) const + { + return x == other.x && y == other.y; + } +}; + +uint qHash(const BinKey &key, uint seed = 0) +{ + return ::qHash(key.x, seed) ^ (::qHash(key.y, seed + 0x9e3779b9U) << 1); +} + +static double distanceSquared(double ax, double ay, double bx, double by) +{ + const double dx = ax - bx; + const double dy = ay - by; + return dx * dx + dy * dy; +} + +static int gridIndex(int x, int y, int width) +{ + return y * width + x; +} + +static QVector gaussianKernel(double sigma) +{ + sigma = std::max(0.1, sigma); + const int radius = std::max(1, static_cast(std::ceil(sigma * 3.0))); + QVector kernel(radius * 2 + 1); + double sum = 0.0; + for (int i = -radius; i <= radius; ++i) { + const double v = std::exp(-(i * i) / (2.0 * sigma * sigma)); + kernel[i + radius] = static_cast(v); + sum += v; + } + if (sum > 0.0) { + for (float &v : kernel) v = static_cast(v / sum); + } + return kernel; +} + +static void smoothMasked(CScanGridResult &grid, double sigma) +{ + if (!grid.valid || grid.width <= 0 || grid.height <= 0) return; + + const QVector kernel = gaussianKernel(sigma); + const int radius = kernel.size() / 2; + const int count = grid.width * grid.height; + + QVector tmpValues(count, 0.0f); + QVector tmpWeights(count, 0.0f); + + for (int y = 0; y < grid.height; ++y) { + for (int x = 0; x < grid.width; ++x) { + double sum = 0.0; + double weight = 0.0; + for (int k = -radius; k <= radius; ++k) { + const int sx = x + k; + if (sx < 0 || sx >= grid.width) continue; + const int idx = gridIndex(sx, y, grid.width); + if (!grid.validMask[idx]) continue; + const double w = kernel[k + radius]; + sum += grid.values[idx] * w; + weight += w; + } + const int outIdx = gridIndex(x, y, grid.width); + if (weight > 1e-6) { + tmpValues[outIdx] = static_cast(sum / weight); + tmpWeights[outIdx] = static_cast(weight); + } + } + } + + QVector outValues(count, 0.0f); + QVector outMask(count, 0); + + for (int y = 0; y < grid.height; ++y) { + for (int x = 0; x < grid.width; ++x) { + double sum = 0.0; + double weight = 0.0; + for (int k = -radius; k <= radius; ++k) { + const int sy = y + k; + if (sy < 0 || sy >= grid.height) continue; + const int idx = gridIndex(x, sy, grid.width); + if (tmpWeights[idx] <= 1e-6f) continue; + const double w = kernel[k + radius] * tmpWeights[idx]; + sum += tmpValues[idx] * w; + weight += w; + } + const int outIdx = gridIndex(x, y, grid.width); + if (weight > 1e-4) { + outValues[outIdx] = static_cast(sum / weight); + outMask[outIdx] = 1; + } + } + } + + grid.values = std::move(outValues); + grid.validMask = std::move(outMask); +} + +static float percentile(QVector values, double percent) +{ + if (values.isEmpty()) return 0.0f; + std::sort(values.begin(), values.end()); + percent = std::clamp(percent, 0.0, 100.0); + const double pos = percent / 100.0 * (values.size() - 1); + const int lo = static_cast(std::floor(pos)); + const int hi = static_cast(std::ceil(pos)); + if (lo == hi) return values[lo]; + const double t = pos - lo; + return static_cast(values[lo] * (1.0 - t) + values[hi] * t); +} + +static void updateDisplayRange(CScanGridResult &grid, double lowPercent, double highPercent) +{ + QVector validValues; + validValues.reserve(grid.values.size()); + for (int i = 0; i < grid.values.size(); ++i) { + if (i < grid.validMask.size() && grid.validMask[i]) { + validValues.append(grid.values[i]); + } + } + + if (validValues.isEmpty()) { + grid.displayMin = 0.0f; + grid.displayMax = 1.0f; + return; + } + + grid.displayMin = percentile(validValues, lowPercent); + grid.displayMax = percentile(validValues, highPercent); + if (std::abs(grid.displayMax - grid.displayMin) < 1e-6f) { + grid.displayMax = grid.displayMin + 1.0f; + } +} + +} // namespace + +CScanGridResult CScanGridder::build(const QVector &surveyLines, CScanGridOptions options) +{ + CScanGridResult result; + + options.cellSizeM = std::max(0.005, options.cellSizeM); + options.searchRadiusM = std::max(options.cellSizeM, options.searchRadiusM); + options.idwPower = std::max(0.1, options.idwPower); + options.smoothingSigmaCells = std::clamp(options.smoothingSigmaCells, 0.1, 5.0); + if (options.clipLowPercent >= options.clipHighPercent) { + options.clipLowPercent = 2.0; + options.clipHighPercent = 98.0; + } + options.maxGridWidth = std::max(64, options.maxGridWidth); + options.maxGridHeight = std::max(64, options.maxGridHeight); + options.maxGridCells = std::max(4096, options.maxGridCells); + + QVector samples; + double minX = std::numeric_limits::max(); + double minY = std::numeric_limits::max(); + double maxX = std::numeric_limits::lowest(); + double maxY = std::numeric_limits::lowest(); + + for (int li = 0; li < surveyLines.size(); ++li) { + const SurveyLine &line = surveyLines[li]; + if (!line.hasValidTrajectories()) continue; + + const GPRDataModel &data = line.data; + + const int channelCount = std::min(static_cast(line.channelTrajectories.size()), data.getChannelCount()); + const int traceCount = data.getTraceCountPerChannel(); + const int sampleCount = data.getSampleCount(); + if (channelCount <= 0 || traceCount <= 0 || sampleCount <= 0) continue; + + const int sampleIndex = qBound(0, options.sampleIndex, sampleCount - 1); + for (int c = 0; c < channelCount; ++c) { + const QVector &traj = line.channelTrajectories[c]; + const int n = std::min(traceCount, static_cast(traj.size())); + for (int t = 0; t < n; ++t) { + const RadarTrace *trace = data.getTrace(c, t); + if (!trace || sampleIndex >= trace->amplitudes.size()) continue; + + const QVector3D &pos = traj[t]; + const float amp = static_cast(trace->amplitudes[sampleIndex]); + // Grid X uses CGCS easting (Y), grid Y uses CGCS northing (X), + // so the generated C-scan image aligns with map screen axes. + samples.append({pos.y(), pos.x(), amp, li, c, t}); + minX = std::min(minX, static_cast(pos.y())); + minY = std::min(minY, static_cast(pos.x())); + maxX = std::max(maxX, static_cast(pos.y())); + maxY = std::max(maxY, static_cast(pos.x())); + } + } + } + + if (samples.isEmpty() || maxX <= minX || maxY <= minY) { + return result; + } + + result.dataMinX = minX; + result.dataMinY = minY; + result.dataMaxX = maxX; + result.dataMaxY = maxY; + + minX -= options.searchRadiusM; + minY -= options.searchRadiusM; + maxX += options.searchRadiusM; + maxY += options.searchRadiusM; + + double cellSize = options.cellSizeM; + int width = static_cast(std::ceil((maxX - minX) / cellSize)) + 1; + int height = static_cast(std::ceil((maxY - minY) / cellSize)) + 1; + while ((width > options.maxGridWidth || height > options.maxGridHeight + || static_cast(width) * height > options.maxGridCells) && cellSize < 10.0) { + cellSize *= 1.5; + width = static_cast(std::ceil((maxX - minX) / cellSize)) + 1; + height = static_cast(std::ceil((maxY - minY) / cellSize)) + 1; + } + + if (width <= 0 || height <= 0 || width > options.maxGridWidth || height > options.maxGridHeight + || static_cast(width) * height > options.maxGridCells) { + return result; + } + + result.valid = true; + result.originX = minX; + result.originY = minY; + result.cellSizeM = cellSize; + result.width = width; + result.height = height; + + const int cellCount = width * height; + result.values.fill(0.0f, cellCount); + result.validMask.fill(0, cellCount); + result.nearestLine.fill(-1, cellCount); + result.nearestChannel.fill(-1, cellCount); + result.nearestTrace.fill(-1, cellCount); + + const double binSize = options.searchRadiusM; + QHash> bins; + bins.reserve(samples.size()); + for (int i = 0; i < samples.size(); ++i) { + const CScanSamplePoint &p = samples[i]; + const BinKey key{static_cast(std::floor((p.x - minX) / binSize)), + static_cast(std::floor((p.y - minY) / binSize))}; + bins[key].append(i); + } + + const double radiusSq = options.searchRadiusM * options.searchRadiusM; + const int searchBins = std::max(1, static_cast(std::ceil(options.searchRadiusM / binSize))); + + for (int y = 0; y < height; ++y) { + const double cy = minY + y * cellSize; + const int by = static_cast(std::floor((cy - minY) / binSize)); + for (int x = 0; x < width; ++x) { + const double cx = minX + x * cellSize; + const int bx = static_cast(std::floor((cx - minX) / binSize)); + + double sum = 0.0; + double weightSum = 0.0; + double bestDistSq = std::numeric_limits::max(); + int bestIndex = -1; + bool exact = false; + float exactValue = 0.0f; + + for (int oy = -searchBins; oy <= searchBins; ++oy) { + for (int ox = -searchBins; ox <= searchBins; ++ox) { + const BinKey key{bx + ox, by + oy}; + auto it = bins.constFind(key); + if (it == bins.constEnd()) continue; + + for (int sampleIdx : it.value()) { + const CScanSamplePoint &p = samples[sampleIdx]; + const double d2 = distanceSquared(cx, cy, p.x, p.y); + if (d2 > radiusSq) continue; + if (d2 < bestDistSq) { + bestDistSq = d2; + bestIndex = sampleIdx; + } + if (d2 < 1e-8) { + exact = true; + exactValue = p.amplitude; + bestIndex = sampleIdx; + break; + } + const double w = 1.0 / std::pow(std::sqrt(d2), options.idwPower); + sum += p.amplitude * w; + weightSum += w; + } + if (exact) break; + } + if (exact) break; + } + + const int idx = gridIndex(x, y, width); + if (exact || weightSum > 0.0) { + result.values[idx] = exact ? exactValue : static_cast(sum / weightSum); + result.validMask[idx] = 1; + if (bestIndex >= 0) { + const CScanSamplePoint &p = samples[bestIndex]; + result.nearestLine[idx] = p.line; + result.nearestChannel[idx] = p.channel; + result.nearestTrace[idx] = p.trace; + } + } + } + } + + if (options.smoothingEnabled) { + smoothMasked(result, options.smoothingSigmaCells); + } + updateDisplayRange(result, options.clipLowPercent, options.clipHighPercent); + + return result; +} diff --git a/external/gpr3dviewer/CScanGridder.h b/external/gpr3dviewer/CScanGridder.h new file mode 100644 index 0000000..aada1e6 --- /dev/null +++ b/external/gpr3dviewer/CScanGridder.h @@ -0,0 +1,48 @@ +#ifndef CSCANGRIDDER_H +#define CSCANGRIDDER_H + +#include "GPRDataModel.h" + +#include + +struct CScanGridOptions { + int sampleIndex = 0; + double cellSizeM = 0.05; + double searchRadiusM = 0.50; + double idwPower = 2.0; + bool smoothingEnabled = true; + double smoothingSigmaCells = 0.9; + int maxGridWidth = 2048; + int maxGridHeight = 2048; + int maxGridCells = 250000; + double clipLowPercent = 2.0; + double clipHighPercent = 98.0; +}; + +struct CScanGridResult { + bool valid = false; + double originX = 0.0; + double originY = 0.0; + double dataMinX = 0.0; + double dataMinY = 0.0; + double dataMaxX = 0.0; + double dataMaxY = 0.0; + double cellSizeM = 0.05; + int width = 0; + int height = 0; + QVector values; + QVector validMask; + QVector nearestLine; + QVector nearestChannel; + QVector nearestTrace; + float displayMin = 0.0f; + float displayMax = 1.0f; +}; + +class CScanGridder +{ +public: + static CScanGridResult build(const QVector &surveyLines, CScanGridOptions options); +}; + +#endif // CSCANGRIDDER_H diff --git a/external/gpr3dviewer/CoordinateTransform.cpp b/external/gpr3dviewer/CoordinateTransform.cpp new file mode 100644 index 0000000..be8141e --- /dev/null +++ b/external/gpr3dviewer/CoordinateTransform.cpp @@ -0,0 +1,219 @@ +#include "CoordinateTransform.h" +#include + +// 静态辅助：预计算的子午线弧长系数 +static inline double computeA0(double e2) { + return 1.0 + 3.0/4.0 * e2 + 45.0/64.0 * e2*e2 + 175.0/256.0 * e2*e2*e2 + + 11025.0/16384.0 * e2*e2*e2*e2; +} +static inline double computeA2(double e2) { + return 3.0/4.0 * e2 + 15.0/16.0 * e2*e2 + 525.0/512.0 * e2*e2*e2 + + 2205.0/2048.0 * e2*e2*e2*e2; +} +static inline double computeA4(double e2) { + return 15.0/64.0 * e2*e2 + 105.0/256.0 * e2*e2*e2 + + 2205.0/4096.0 * e2*e2*e2*e2; +} +static inline double computeA6(double e2) { + return 35.0/512.0 * e2*e2*e2 + 315.0/2048.0 * e2*e2*e2*e2; +} +static inline double computeA8(double e2) { + return 315.0/16384.0 * e2*e2*e2*e2; +} + +int CoordinateTransform::detectZoneFromCgcsY(double cgcsY) +{ + // 中国范围内 CGCS2000 坐标 Y 通常为正值且包含带号 + // 格式：zone * 1,000,000 + 500,000 + easting + double absY = std::abs(cgcsY); + if (absY > 1e6) { + int zone = static_cast(std::floor(absY / 1e6)); + if (zone >= 1 && zone <= 60) + return zone; + } + // 如果无法推断，返回 0（表示需要手动指定） + return 0; +} + +double CoordinateTransform::centralMeridianFromZone(int zone) +{ + // 3度带：中央经线 = 带号 * 3° + return zone * 3.0; +} + +double CoordinateTransform::meridianArc(double latRad) +{ + const double e2 = E2; + const double a = A; + + const double A0 = computeA0(e2); + const double A2 = computeA2(e2); + const double A4 = computeA4(e2); + const double A6 = computeA6(e2); + const double A8 = computeA8(e2); + + double s2 = std::sin(2.0 * latRad); + double s4 = std::sin(4.0 * latRad); + double s6 = std::sin(6.0 * latRad); + double s8 = std::sin(8.0 * latRad); + + return a * (1.0 - e2) * (A0 * latRad - A2/2.0 * s2 + A4/4.0 * s4 + - A6/6.0 * s6 + A8/8.0 * s8); +} + +double CoordinateTransform::meridianArcInverse(double x) +{ + const double e2 = E2; + const double a = A; + + const double A0 = computeA0(e2); + const double A2 = computeA2(e2); + const double A4 = computeA4(e2); + const double A6 = computeA6(e2); + const double A8 = computeA8(e2); + + const double c0 = a * (1.0 - e2) * A0; + + double bf = x / c0; + for (int i = 0; i < 8; ++i) { + double s2 = std::sin(2.0 * bf); + double s4 = std::sin(4.0 * bf); + double s6 = std::sin(6.0 * bf); + double s8 = std::sin(8.0 * bf); + + double fx = a * (1.0 - e2) * (A0 * bf - A2/2.0 * s2 + A4/4.0 * s4 + - A6/6.0 * s6 + A8/8.0 * s8) - x; + double fpx = a * (1.0 - e2) * (A0 - A2 * std::cos(2.0 * bf) + + A4 * std::cos(4.0 * bf) + - A6 * std::cos(6.0 * bf) + + A8 * std::cos(8.0 * bf)); + + double delta = fx / fpx; + bf -= delta; + if (std::abs(delta) < 1e-12) + break; + } + return bf; +} + +void CoordinateTransform::cgcs2000ToWgs84(double cgcsX, double cgcsY, int zone, + double &latRad, double &lonRad) +{ + const double a = A; + const double e2 = E2; + const double e12 = E12; + + // 提取东偏移（去掉带号和 500km 常数） + double easting = cgcsY; + if (zone > 0) { + easting = cgcsY - zone * 1e6 - FALSE_EASTING; + } else { + easting = cgcsY - FALSE_EASTING; + } + double northing = cgcsX; + double l0 = centralMeridianFromZone(zone) * DEG2RAD; + + // 底点纬度 + double bf = meridianArcInverse(northing); + + double sinBf = std::sin(bf); + double cosBf = std::cos(bf); + double tanBf = std::tan(bf); + + double Nf = a / std::sqrt(1.0 - e2 * sinBf * sinBf); + double etaf2 = e12 * cosBf * cosBf; + double tf = tanBf; + double Vf2 = 1.0 + etaf2; + + double y = easting; + double y2 = y * y; + double y4 = y2 * y2; + double y6 = y4 * y2; + + double Nf2 = Nf * Nf; + double Nf4 = Nf2 * Nf2; + double Nf6 = Nf4 * Nf2; + + // 纬度修正 + double dB = tf / Vf2 * ( + y2 / (2.0 * Nf2) * (1.0 + etaf2) + - y4 / (24.0 * Nf4) * (5.0 + 3.0 * tf*tf + 6.0 * etaf2 + - 6.0 * tf*tf * etaf2 + - 3.0 * etaf2*etaf2 + - 9.0 * tf*tf * etaf2*etaf2) + + y6 / (720.0 * Nf6) * (61.0 + 90.0 * tf*tf + 45.0 * tf*tf*tf*tf) + ); + + latRad = bf - dB; + + // 经差 + double dl = y / (Nf * cosBf) * ( + 1.0 + - y2 / (6.0 * Nf2) * (1.0 + 2.0 * tf*tf + etaf2) + + y4 / (120.0 * Nf4) * (5.0 + 28.0 * tf*tf + 24.0 * tf*tf*tf*tf + + 6.0 * etaf2 + 8.0 * tf*tf * etaf2) + - y6 / (5040.0 * Nf6) * (61.0 + 662.0 * tf*tf + 1320.0 * tf*tf*tf*tf + + 720.0 * tf*tf*tf*tf*tf*tf) + ); + + lonRad = l0 + dl; +} + +void CoordinateTransform::wgs84ToCgcs2000(double latRad, double lonRad, int zone, + double &cgcsX, double &cgcsY) +{ + const double a = A; + const double e2 = E2; + const double e12 = E12; + + double l0 = centralMeridianFromZone(zone) * DEG2RAD; + double l = lonRad - l0; + + double sinB = std::sin(latRad); + double cosB = std::cos(latRad); + double tanB = std::tan(latRad); + + double N = a / std::sqrt(1.0 - e2 * sinB * sinB); + double t2 = tanB * tanB; + double t4 = t2 * t2; + double eta2 = e12 * cosB * cosB; + double eta4 = eta2 * eta2; + + double l2 = l * l; + double l4 = l2 * l2; + double l6 = l4 * l2; + + // 子午线弧长 + double X = meridianArc(latRad); + + // 北坐标 + cgcsX = X + N / 2.0 * sinB * cosB * l2 * ( + 1.0 + l2 / 12.0 * (5.0 - t2 + 9.0 * eta2 + 4.0 * eta4) + + l4 / 360.0 * (61.0 - 58.0 * t2 + t4) + ); + + // 东坐标（自然值，不含带号和 500km） + double yNatural = N * cosB * l * ( + 1.0 + l2 / 6.0 * (1.0 - t2 + eta2) + + l4 / 120.0 * (5.0 - 18.0 * t2 + t4 + 14.0 * eta2 - 58.0 * t2 * eta2) + ); + + // 加上带号和 500km 常数 + if (zone > 0) { + cgcsY = zone * 1e6 + FALSE_EASTING + yNatural; + } else { + cgcsY = FALSE_EASTING + yNatural; + } +} + +void CoordinateTransform::wgs84ToWebMercator(double latRad, double lonRad, double &mx, double &my) +{ + mx = EARTH_RADIUS * lonRad; + my = EARTH_RADIUS * std::log(std::tan(PI / 4.0 + latRad / 2.0)); +} + +void CoordinateTransform::webMercatorToWgs84(double mx, double my, double &latRad, double &lonRad) +{ + lonRad = mx / EARTH_RADIUS; + latRad = 2.0 * std::atan(std::exp(my / EARTH_RADIUS)) - PI / 2.0; +} diff --git a/external/gpr3dviewer/CoordinateTransform.h b/external/gpr3dviewer/CoordinateTransform.h new file mode 100644 index 0000000..cfca21b --- /dev/null +++ b/external/gpr3dviewer/CoordinateTransform.h @@ -0,0 +1,48 @@ +#ifndef COORDINATETRANSFORM_H +#define COORDINATETRANSFORM_H + +#include + +class CoordinateTransform { +public: + // CGCS2000 椭球参数 + static constexpr double A = 6378137.0; // 长半轴 (m) + static constexpr double F = 1.0 / 298.257222101; // 扁率 + static constexpr double E2 = 2.0 * F - F * F; // 第一偏心率平方 + static constexpr double E12 = E2 / (1.0 - E2); // 第二偏心率平方 + static constexpr double FALSE_EASTING = 500000.0; // 东偏移 (m) + + static constexpr double EARTH_RADIUS = 6378137.0; // Web Mercator 地球半径 + static constexpr double PI = 3.14159265358979323846; + static constexpr double DEG2RAD = PI / 180.0; + static constexpr double RAD2DEG = 180.0 / PI; + + // 从 CGCS2000 东向坐标(Y)自动推断 3° 带带号 + // 假设 Y 格式为：zone*1e6 + 500000 + easting，或纯自然值 + static int detectZoneFromCgcsY(double cgcsY); + + // 带号 → 中央经线（度） + static double centralMeridianFromZone(int zone); + + // CGCS2000 平面坐标 → WGS84 经纬度（弧度） + // cgcsX: 北向坐标(m), cgcsY: 东向坐标(m, 可含带号) + static void cgcs2000ToWgs84(double cgcsX, double cgcsY, int zone, + double &latRad, double &lonRad); + + // WGS84 经纬度（弧度） → CGCS2000 平面坐标 + static void wgs84ToCgcs2000(double latRad, double lonRad, int zone, + double &cgcsX, double &cgcsY); + + // WGS84 经纬度（弧度） ↔ Web Mercator 平面坐标（米） + static void wgs84ToWebMercator(double latRad, double lonRad, double &mx, double &my); + static void webMercatorToWgs84(double mx, double my, double &latRad, double &lonRad); + +private: + // 子午线弧长（从赤道到纬度 B 的弧长） + static double meridianArc(double latRad); + + // 子午线弧长反解：已知弧长 X，求底点纬度 Bf + static double meridianArcInverse(double x); +}; + +#endif // COORDINATETRANSFORM_H diff --git a/external/gpr3dviewer/PosParser.cpp b/external/gpr3dviewer/PosParser.cpp new file mode 100644 index 0000000..0538978 --- /dev/null +++ b/external/gpr3dviewer/PosParser.cpp @@ -0,0 +1,73 @@ +#include "PosParser.h" +#include +#include +#include +#include + +bool PosParser::loadPosFile(const QString &posFilePath, QVector &outPositions) { + outPositions.clear(); + + QFile file(posFilePath); + if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) { + qDebug() << "PosParser: Cannot open pos file:" << posFilePath; + return false; + } + + QTextStream in(&file); + QRegularExpression reLineComment("^\\s*%"); + + while (!in.atEnd()) { + QString line = in.readLine().trimmed(); + if (line.isEmpty()) continue; + if (reLineComment.match(line).hasMatch()) continue; // 跳过注释行 + + QStringList parts = line.split(QRegularExpression("\\s+"), Qt::SkipEmptyParts); + if (parts.size() < 4) continue; + + bool okX = false, okY = false, okZ = false; + // 列顺序: 道号(忽略) X Y Z + double x = parts[1].toDouble(&okX); + double y = parts[2].toDouble(&okY); + double z = parts[3].toDouble(&okZ); + + if (okX && okY && okZ) { + outPositions.append(QVector3D(static_cast(x), static_cast(y), static_cast(z))); + } + } + + file.close(); + qDebug() << "PosParser: Loaded" << outPositions.size() << "GPS points from" << posFilePath; + return !outPositions.isEmpty(); +} + +bool PosParser::loadCenterCcc(const QString &cccFilePath, double &outLat, double &outLon) { + QFile file(cccFilePath); + if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) { + qDebug() << "PosParser: Cannot open center.ccc:" << cccFilePath; + return false; + } + + QTextStream in(&file); + QRegularExpression reLineComment("^\\s*%"); + + while (!in.atEnd()) { + QString line = in.readLine().trimmed(); + if (line.isEmpty()) continue; + if (reLineComment.match(line).hasMatch()) continue; + + QStringList parts = line.split(QRegularExpression("\\s+"), Qt::SkipEmptyParts); + if (parts.size() < 2) continue; + + bool okLat = false, okLon = false; + outLat = parts[0].toDouble(&okLat); + outLon = parts[1].toDouble(&okLon); + + if (okLat && okLon) { + file.close(); + return true; + } + } + + file.close(); + return false; +} diff --git a/external/gpr3dviewer/PosParser.h b/external/gpr3dviewer/PosParser.h new file mode 100644 index 0000000..cd1353c --- /dev/null +++ b/external/gpr3dviewer/PosParser.h @@ -0,0 +1,17 @@ +#ifndef POSPARSER_H +#define POSPARSER_H + +#include +#include +#include + +class PosParser { +public: + // 加载 .pos 文件（4列：道号 X Y Z） + static bool loadPosFile(const QString &posFilePath, QVector &outPositions); + + // 加载 center.ccc 项目中心坐标（2列：纬度 经度） + static bool loadCenterCcc(const QString &cccFilePath, double &outLat, double &outLon); +}; + +#endif // POSPARSER_H diff --git a/external/gpr3dviewer/TrajectoryCalculator.cpp b/external/gpr3dviewer/TrajectoryCalculator.cpp new file mode 100644 index 0000000..2516454 --- /dev/null +++ b/external/gpr3dviewer/TrajectoryCalculator.cpp @@ -0,0 +1,242 @@ +#include "TrajectoryCalculator.h" +#include +#include + +namespace { + +double planarDistance(const QVector3D &a, const QVector3D &b) +{ + const double dx = static_cast(a.x()) - b.x(); + const double dy = static_cast(a.y()) - b.y(); + return std::sqrt(dx * dx + dy * dy); +} + +QVector3D lerpPoint(const QVector3D &a, const QVector3D &b, double t) +{ + return a * static_cast(1.0 - t) + b * static_cast(t); +} + +QVector3D catmullRom(const QVector3D &p0, const QVector3D &p1, const QVector3D &p2, const QVector3D &p3, double t) +{ + const float tf = static_cast(t); + const float t2 = tf * tf; + const float t3 = t2 * tf; + return 0.5f * ((2.0f * p1) + (-p0 + p2) * tf + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * t2 + + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * t3); +} + +} // namespace + +double TrajectoryCalculator::headingAt(int traceIdx, const QVector &gpsPositions) +{ + const int n = gpsPositions.size(); + if (n <= 1) return 0.0; + + if (traceIdx < n - 1) { + double deltaX = gpsPositions[traceIdx + 1].x() - gpsPositions[traceIdx].x(); // 北向变化 + double deltaY = gpsPositions[traceIdx + 1].y() - gpsPositions[traceIdx].y(); // 东向变化 + return std::atan2(deltaY, deltaX); + } else { + // 最后一点复用前一点方向 + double deltaX = gpsPositions[traceIdx].x() - gpsPositions[traceIdx - 1].x(); + double deltaY = gpsPositions[traceIdx].y() - gpsPositions[traceIdx - 1].y(); + return std::atan2(deltaY, deltaX); + } +} + +bool TrajectoryCalculator::computeTrajectories(SurveyLine &line, const SurveyGeometry &geom) +{ + const QVector &gps = line.gpsPositions; + const int nRtk = gps.size(); + if (nRtk <= 0) return false; + + const int chCount = geom.channelCount; + if (chCount <= 0) return false; + + // 通道相对x偏移优先使用头文件CH_X_OFFSETS推导出的逐通道数组；旧工程无数组时按首通道对称补齐。 + QVector chXRel = geom.channelXRel; + if (chXRel.size() != chCount) { + chXRel.resize(chCount); + const double lastXRel = -geom.ch1XRel; + for (int c = 0; c < chCount; ++c) { + chXRel[c] = (chCount > 1) + ? geom.ch1XRel + (lastXRel - geom.ch1XRel) * c / (chCount - 1.0) + : geom.ch1XRel; + } + } + + // 预分配 channelTrajectories + line.channelTrajectories.resize(chCount); + for (int c = 0; c < chCount; ++c) { + line.channelTrajectories[c].resize(nRtk); + } + + // 遍历每个 RTK 点 + for (int i = 0; i < nRtk; ++i) { + double thetaY = headingAt(i, gps); + double thetaX = thetaY + M_PI / 2.0; + + // 旋转矩阵 R = [cos(theta_x), cos(theta_y); sin(theta_x), sin(theta_y)] + // 作用：将设备相对坐标系（x=右, y=前）转换为绝对坐标系（X=北, Y=东） + double r11 = std::cos(thetaX); + double r12 = std::cos(thetaY); + double r21 = std::sin(thetaX); + double r22 = std::sin(thetaY); + + // RTK 相对天线中心的偏移向量（设备坐标系） + double rtkLocalX = geom.rtkOffsetX; // 横向偏移 + double rtkLocalY = geom.rtkOffsetY; // 纵向偏移（前方为正） + + // 转换为绝对坐标系偏移 + double rtkAbsOffsetX = r11 * rtkLocalX + r12 * rtkLocalY; + double rtkAbsOffsetY = r21 * rtkLocalX + r22 * rtkLocalY; + + // 天线中心绝对坐标 = RTK - 偏移量 + double antX = gps[i].x() - rtkAbsOffsetX; + double antY = gps[i].y() - rtkAbsOffsetY; + double antZ = gps[i].z(); + + // 遍历每个通道 + for (int c = 0; c < chCount; ++c) { + // 通道相对坐标（设备坐标系） + double chLocalX = chXRel[c]; + double chLocalY = 0.0; // MATLAB 中 ch_y_rel = zeros(...) + + // 转换为绝对坐标系偏移 + double chAbsOffsetX = r11 * chLocalX + r12 * chLocalY; + double chAbsOffsetY = r21 * chLocalX + r22 * chLocalY; + + // 通道绝对坐标 + double chX = antX + chAbsOffsetX; + double chY = antY + chAbsOffsetY; + double chZ = antZ; + + line.channelTrajectories[c][i] = QVector3D(static_cast(chX), + static_cast(chY), + static_cast(chZ)); + } + } + + // 同步更新 GPRDataModel 中 traces 的 position + // 数据存储顺序：trace0_ch0, trace0_ch1, ... trace0_chN, trace1_ch0, ... + // 即 globalIdx = traceInChannel * chCount + channel + const int tracesPerChannel = line.data.getTracesPerChannel(); + for (int c = 0; c < chCount; ++c) { + for (int t = 0; t < tracesPerChannel && t < nRtk; ++t) { + int globalIdx = t * chCount + c; + if (globalIdx >= 0 && globalIdx < line.data.traces.size()) { + line.data.traces[globalIdx].position = line.channelTrajectories[c][t]; + } + } + } + + return true; +} + +QVector TrajectoryCalculator::resampleTrajectoryLinear(const QVector &input, int targetCount) +{ + QVector output; + if (targetCount <= 0 || input.isEmpty()) return output; + if (input.size() == 1 || targetCount == 1) { + output.resize(targetCount); + std::fill(output.begin(), output.end(), input.first()); + return output; + } + + output.reserve(targetCount); + const double scale = static_cast(input.size() - 1) / qMax(1, targetCount - 1); + for (int i = 0; i < targetCount; ++i) { + const double src = i * scale; + const int i0 = qBound(0, static_cast(std::floor(src)), input.size() - 1); + const int i1 = qBound(0, i0 + 1, input.size() - 1); + output.append(lerpPoint(input[i0], input[i1], src - i0)); + } + return output; +} + +QVector TrajectoryCalculator::resampleTrajectorySpline(const QVector &input, int targetCount) +{ + if (input.size() < 4) return resampleTrajectoryLinear(input, targetCount); + QVector output; + if (targetCount <= 0) return output; + if (targetCount == 1) return QVector{input.first()}; + + output.reserve(targetCount); + const double scale = static_cast(input.size() - 1) / qMax(1, targetCount - 1); + for (int i = 0; i < targetCount; ++i) { + const double src = i * scale; + const int i1 = qBound(0, static_cast(std::floor(src)), input.size() - 1); + const int i2 = qBound(0, i1 + 1, input.size() - 1); + const int i0 = qBound(0, i1 - 1, input.size() - 1); + const int i3 = qBound(0, i2 + 1, input.size() - 1); + output.append(catmullRom(input[i0], input[i1], input[i2], input[i3], src - i1)); + } + return output; +} + +TrajectoryFilterResult TrajectoryCalculator::filterAndInterpolateTrajectory(const QVector &input, + const TrajectoryFilterOptions &options, + int targetCount) +{ + TrajectoryFilterResult result; + const int n = input.size(); + result.keepMask = QVector(n, true); + result.distanceOutlierMask = QVector(n, false); + result.speedOutlierMask = QVector(n, false); + result.angleOutlierMask = QVector(n, false); + if (n <= 0 || targetCount <= 0) return result; + + for (int i = 1; i < n; ++i) { + const double dist = planarDistance(input[i - 1], input[i]); + if (options.distanceFilterEnabled && dist > options.maxSegmentDistanceM) { + result.distanceOutlierMask[i] = true; + result.keepMask[i] = false; + } + if (options.speedFilterEnabled && options.traceIntervalSec > 1e-9 && dist / options.traceIntervalSec > options.maxSpeedMps) { + result.speedOutlierMask[i] = true; + result.keepMask[i] = false; + } + } + + if (options.angleFilterEnabled) { + const double thresholdRad = options.maxTurnAngleDeg * M_PI / 180.0; + for (int i = 1; i + 1 < n; ++i) { + const double d0 = planarDistance(input[i - 1], input[i]); + const double d1 = planarDistance(input[i], input[i + 1]); + if (d0 < 1e-6 || d1 < 1e-6) continue; + const double ax = input[i].x() - input[i - 1].x(); + const double ay = input[i].y() - input[i - 1].y(); + const double bx = input[i + 1].x() - input[i].x(); + const double by = input[i + 1].y() - input[i].y(); + const double cosv = std::clamp((ax * bx + ay * by) / (d0 * d1), -1.0, 1.0); + const double turn = std::acos(cosv); + if (turn > thresholdRad) { + result.angleOutlierMask[i] = true; + result.keepMask[i] = false; + } + } + } + + if (options.preserveEndpoints && n > 0) { + result.keepMask[0] = true; + result.keepMask[n - 1] = true; + } + + QVector kept; + kept.reserve(n); + for (int i = 0; i < n; ++i) { + if (result.keepMask.value(i, true)) kept.append(input[i]); + } + if (kept.size() < 2) { + kept = input; + result.warnings.append(QStringLiteral("过滤后有效轨迹点过少，已使用原始轨迹插值。")); + } + + if (options.interpolationMode == TrajectoryFilterOptions::InterpolationMode::Spline) { + if (kept.size() < 4) result.warnings.append(QStringLiteral("样条插值点数不足，已回退到线性插值。")); + result.outputPositions = resampleTrajectorySpline(kept, targetCount); + } else { + result.outputPositions = resampleTrajectoryLinear(kept, targetCount); + } + return result; +} diff --git a/external/gpr3dviewer/TrajectoryCalculator.h b/external/gpr3dviewer/TrajectoryCalculator.h new file mode 100644 index 0000000..07aeedc --- /dev/null +++ b/external/gpr3dviewer/TrajectoryCalculator.h @@ -0,0 +1,49 @@ +#ifndef TRAJECTORYCALCULATOR_H +#define TRAJECTORYCALCULATOR_H + +#include "GPRDataModel.h" +#include "SurveyGeometry.h" +#include +#include + +struct TrajectoryFilterOptions { + bool distanceFilterEnabled = true; + double maxSegmentDistanceM = 1.0; + bool speedFilterEnabled = false; + double maxSpeedMps = 60.0; + double traceIntervalSec = 0.05; + bool angleFilterEnabled = true; + double maxTurnAngleDeg = 60.0; + enum class InterpolationMode { Linear, Spline }; + InterpolationMode interpolationMode = InterpolationMode::Linear; + bool preserveEndpoints = true; + bool preserveManualEdits = true; +}; + +struct TrajectoryFilterResult { + QVector outputPositions; + QVector keepMask; + QVector distanceOutlierMask; + QVector speedOutlierMask; + QVector angleOutlierMask; + QStringList warnings; +}; + +class TrajectoryCalculator { +public: + // 计算第 traceIdx 个 RTK 点的行进方向角（弧度） + // 返回 theta_y = atan2(dY, dX)，其中 X=北向, Y=东向 + static double headingAt(int traceIdx, const QVector &gpsPositions); + + // 根据 RTK 轨迹和天线几何参数，计算所有通道的绝对坐标 + // 结果写入 line.channelTrajectories 和 line.data.traces[].position + static bool computeTrajectories(SurveyLine &line, const SurveyGeometry &geom); + + static QVector resampleTrajectoryLinear(const QVector &input, int targetCount); + static QVector resampleTrajectorySpline(const QVector &input, int targetCount); + static TrajectoryFilterResult filterAndInterpolateTrajectory(const QVector &input, + const TrajectoryFilterOptions &options, + int targetCount); +}; + +#endif // TRAJECTORYCALCULATOR_H diff --git a/src/io/gpr/CMakeLists.txt b/src/io/gpr/CMakeLists.txt index 1d6bcd2..4afd503 100644 --- a/src/io/gpr/CMakeLists.txt +++ b/src/io/gpr/CMakeLists.txt @@ -9,11 +9,17 @@ target_link_libraries(geopro_io_gpr PUBLIC geopro_core) # 纯 C++17 解析层，零 Qt/VTK。顶层全局开启了 AUTOMOC/UIC/RCC，这里显式关闭保持纯净。 set_target_properties(geopro_io_gpr PROPERTIES AUTOMOC OFF AUTOUIC OFF AUTORCC OFF) -# gpr3dv 桥接层(P2)：把 vendored gpr3dv(P1) 处理后立方体 → geopro 量化体(BuiltI16)。 +# gpr3dv 桥接层(P2/P8)：把 vendored gpr3dv 处理链 → geopro 量化体(BuiltI16)。 # 独立 target(不污染上面的纯 C++17 解析层)：链 geopro_gpr3dv(Qt)+geopro_core。 -add_library(geopro_gpr3dv_bridge STATIC Gpr3dvVolumeBridge.cpp) +# Gpr3dvVolumeBridge —— P2 线局部坐标体(X=道/Y=通道/Z=样本，origin≈0)。 +# Gpr3dvSurveyVolumeBridge —— P8 测绘级 CGCS2000 世界对齐体(逐道跟 GPS 轨迹)； +# 额外用 GpsTrack(geopro_io_gpr) 解析 .gps。 +add_library(geopro_gpr3dv_bridge STATIC + Gpr3dvVolumeBridge.cpp + Gpr3dvSurveyVolumeBridge.cpp) target_include_directories(geopro_gpr3dv_bridge PUBLIC ${CMAKE_SOURCE_DIR}/src) target_compile_features(geopro_gpr3dv_bridge PUBLIC cxx_std_17) -target_link_libraries(geopro_gpr3dv_bridge PUBLIC geopro_core geopro_gpr3dv) +target_link_libraries(geopro_gpr3dv_bridge PUBLIC + geopro_core geopro_gpr3dv geopro_io_gpr) set_target_properties(geopro_gpr3dv_bridge PROPERTIES AUTOMOC OFF AUTOUIC OFF AUTORCC OFF) diff --git a/src/io/gpr/Gpr3dvSurveyVolumeBridge.cpp b/src/io/gpr/Gpr3dvSurveyVolumeBridge.cpp new file mode 100644 index 0000000..62d9d2b --- /dev/null +++ b/src/io/gpr/Gpr3dvSurveyVolumeBridge.cpp @@ -0,0 +1,270 @@ +#include "io/gpr/Gpr3dvSurveyVolumeBridge.hpp" + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "CScanGridder.h" +#include "CoordinateTransform.h" +#include "GPRDataModel.h" +#include "IprhParser.h" +#include "RadarProcessor.h" +#include "SurveyGeometry.h" +#include "TrajectoryCalculator.h" + +#include "core/model/ScalarVolumeI16.hpp" +#include "io/gpr/GpsTrack.hpp" + +namespace geopro::io::gpr { + +namespace { + +// 全体 traces 平均绝对幅值——「处理是否生效」标量证据(与 P2 桥接同口径)。 +double meanAbsAmplitude(const GPRDataModel& model) { + long double sum = 0.0L; + long long count = 0; + for (const RadarTrace& tr : model.traces) { + for (short v : tr.amplitudes) { + sum += static_cast(v < 0 ? -v : v); + ++count; + } + } + return count > 0 ? static_cast(sum / count) : 0.0; +} + +// 深度采样间距(米)：(timeWindow/samples) ns × 波速(m/ns) / 2(往返)。 +// 与 GPRDataModel::calculateDepth 同式。取不到 → 1.0(单位间距)。 +double depthSpacingZ(const GPRDataModel::Header& h) { + if (h.samplesPerTrace <= 0 || h.timeWindowNs <= 0.0) return 1.0; + const double timePerSample = h.timeWindowNs / h.samplesPerTrace; // ns + const double vel = h.waveVelocity > 0.0 ? h.waveVelocity : 0.1; // m/ns + const double dz = timePerSample * vel / 2.0; + return dz > 1e-12 ? dz : 1.0; +} + +double nowMs(std::chrono::steady_clock::time_point t0) { + const auto t1 = std::chrono::steady_clock::now(); + return std::chrono::duration(t1 - t0).count(); +} + +constexpr double kDeg2Rad = 3.14159265358979323846 / 180.0; + +} // namespace + +geopro::core::BuiltI16 buildLineVolumeSurvey(const std::string& lineDir, + const std::string& linePrefix, + const std::string& gpsPath, + SurveyBridgeMetrics* metricsOut, + int coarse, double cellSizeM, + double searchRadiusM) { + const int stride = coarse > 1 ? coarse : 1; // 深度下采样步长(≥1) + const QString dir = QString::fromLocal8Bit(lineDir.c_str()); + const QString base = QString::fromLocal8Bit(linePrefix.c_str()); + + // 1) 走 P1 链：load → buildVolumeData(处理前) → runPipeline → 再 buildVolumeData。 + const auto tLoad = std::chrono::steady_clock::now(); + GPRDataModel model; + if (!IprhParser::loadImpulseMultiChannel(dir, base, model)) { + throw std::runtime_error("loadImpulseMultiChannel 失败: " + lineDir + " / " + + linePrefix); + } + model.buildVolumeData(); + const double meanBefore = meanAbsAmplitude(model); + const double loadMs = nowMs(tLoad); + + const auto tPipe = std::chrono::steady_clock::now(); + RadarProcessor::ProcPipeline pipeline; + pipeline.setDefaultFlow(); + GPRDataModel processed = RadarProcessor::runPipeline(model, pipeline); + // 零时校正会改 samplesPerTrace，且 runPipeline 不重建 volumeData → 必须再建一次。 + processed.buildVolumeData(); + const double meanAfter = meanAbsAmplitude(processed); + const double pipelineMs = nowMs(tPipe); + + // 2) 几何 + GPS→CGCS2000 世界坐标(逐 RTK 点)。 + const auto tTraj = std::chrono::steady_clock::now(); + // 值拷贝(非引用)：下面 processed 会被 std::move 进 SurveyLine，之后仍要用 header 算 dz， + // 引用会变悬垂(use-after-move)。Header 是小结构(几个 QVector/QString)，拷贝代价可忽略。 + const GPRDataModel::Header h = processed.header; + const int channels = h.numberOfChannels > 0 ? h.numberOfChannels : 1; + + // SurveyGeometry：逐通道天线横偏(归一到天线中心)，来自头 CH_X/Y_OFFSETS。 + SurveyGeometry geom = + SurveyGeometry::fromHeaderOffsets(channels, h.chXOffsets, h.chYOffsets); + + // .gps 经纬(度) → CGCS2000 平面(北X, 东Y) 米。带号由首点经度自动定(3°带)。 + const geopro::io::gpr::GpsTrack track = geopro::io::gpr::parseGps(gpsPath); + if (track.pts.size() < 2) { + throw std::runtime_error("GPS 轨迹点 <2(无法定位/航向): " + gpsPath); + } + const double lon0Deg = track.pts.front().lon; + // 3°带号 = round(中央经线/3)；中央经线最近 3° 倍数。 + const int zone = static_cast(std::lround(lon0Deg / 3.0)); + const double centralMeridianDeg = CoordinateTransform::centralMeridianFromZone(zone); + geom.cgcsZone = zone; + geom.centralMeridianDeg = centralMeridianDeg; + + QVector gpsCgcs; + gpsCgcs.reserve(static_cast(track.pts.size())); + for (const auto& p : track.pts) { + double cx = 0.0, cy = 0.0; // cx=北(northing), cy=东(easting，含带号) + CoordinateTransform::wgs84ToCgcs2000(p.lat * kDeg2Rad, p.lon * kDeg2Rad, + zone, cx, cy); + // SurveyLine.gpsPositions 约定 x=北, y=东, z=高(与 TrajectoryCalculator 一致)。 + gpsCgcs.append(QVector3D(static_cast(cx), static_cast(cy), + static_cast(p.elev))); + } + + SurveyLine line; + line.data = std::move(processed); + line.geometry = geom; + line.gpsPositions = gpsCgcs; + if (!TrajectoryCalculator::computeTrajectories(line, geom)) { + throw std::runtime_error("computeTrajectories 失败: " + linePrefix); + } + if (!line.hasValidTrajectories()) { + throw std::runtime_error("无有效通道轨迹: " + linePrefix); + } + const double trajMs = nowMs(tTraj); + + const int samples = line.data.getSampleCount(); + if (samples <= 0) { + throw std::runtime_error("处理后样本数为 0: " + linePrefix); + } + + // 3) 逐深度 IDW 世界网格(CScanGridder)→ 堆成世界轴对齐体。 + // 全部深度共享同一水平网格(dims 仅由轨迹 AABB + cellSize 定，与深度无关)， + // 故用深度 0 的网格几何作体的世界框，逐深度填一层。 + const auto tGrid = std::chrono::steady_clock::now(); + const QVector lines{line}; + + CScanGridOptions baseOpts; + if (cellSizeM > 0.0) baseOpts.cellSizeM = cellSizeM; + if (searchRadiusM > 0.0) baseOpts.searchRadiusM = searchRadiusM; + // coarse：水平格距 ×coarse 省空间(searchRadius 同步放大保证邻域覆盖)。 + if (stride > 1) { + baseOpts.cellSizeM *= stride; + baseOpts.searchRadiusM = + std::max(baseOpts.searchRadiusM, baseOpts.cellSizeM * 2.0); + } + // 关闭逐层百分位裁剪对体值无副作用(displayMin/Max 仅显示用，不改 values)； + // 平滑保留(与原版 C-scan 一致)。 + + // 输出深度采样数(coarse 沿深度下采样)：nzOut = ceil(samples/stride)。 + const int nzOut = (samples + stride - 1) / stride; + + // 先建深度 0 网格锁定世界框。 + CScanGridOptions opt0 = baseOpts; + opt0.sampleIndex = 0; + const CScanGridResult g0 = CScanGridder::build(lines, opt0); + if (!g0.valid || g0.width <= 0 || g0.height <= 0) { + throw std::runtime_error("CScanGridder 深度0 网格无效(轨迹/数据为空?): " + + linePrefix); + } + const int nx = g0.width; // 东向 easting 格数 + const int ny = g0.height; // 北向 northing 格数 + const int nz = nzOut; // 深度 + const double cellSize = g0.cellSizeM; + const double originX = g0.originX; // CGCS easting(含带号) 米 + const double originY = g0.originY; // CGCS northing 米 + + // 量化区间：先扫所有输出深度网格的有效值(命中 validMask)定 [vmin,vmax]。 + // (CScanGridder values 为处理后幅值的 IDW 插值结果，连续浮点。) + double vmin = std::numeric_limits::infinity(); + double vmax = -std::numeric_limits::infinity(); + // 缓存每层网格(避免二次计算)：网格层数 = nzOut，每层 nx*ny float，明星路约 + // 数百×数千×数百 → 单线峰值可控(逐线建，建完即释放)。 + std::vector layers; + layers.reserve(static_cast(nzOut)); + std::int64_t filledCells = 0; + for (int zo = 0; zo < nzOut; ++zo) { + const int s = zo * stride; + CScanGridOptions opt = baseOpts; + opt.sampleIndex = std::min(s, samples - 1); + CScanGridResult g = + (zo == 0) ? g0 : CScanGridder::build(lines, opt); + // 守护：所有深度网格 dims 必须与深度0一致(理论恒等)；不一致则该层置空跳过。 + if (!g.valid || g.width != nx || g.height != ny) { + g.valid = false; + } else { + for (int i = 0; i < g.values.size() && i < g.validMask.size(); ++i) { + if (!g.validMask[i]) continue; + ++filledCells; + const double v = static_cast(g.values[i]); + if (v < vmin) vmin = v; + if (v > vmax) vmax = v; + } + } + layers.push_back(std::move(g)); + } + if (!(vmin <= vmax)) { // 无任何有效值(理论不至)。 + vmin = 0.0; + vmax = 0.0; + } + + geopro::core::Quant quant; + quant.scale = (vmax > vmin) ? (vmax - vmin) / 64000.0 : 1.0; + quant.offset = 0.5 * (vmin + vmax); // 中点 → 防溢出 + + // 4) 填体：逐层逐 cell。空 cell(validMask=0)→ kBlank(渲染透明)。 + // 轴 X=东(easting), Y=北(northing), Z=深度。 + geopro::core::BuiltI16 built; + built.vol = geopro::core::ScalarVolumeI16(nx, ny, nz); + for (int zo = 0; zo < nzOut; ++zo) { + const CScanGridResult& g = layers[static_cast(zo)]; + for (int y = 0; y < ny; ++y) { + for (int x = 0; x < nx; ++x) { + const int gi = y * nx + x; + std::int16_t q = geopro::core::ScalarVolumeI16::kBlank; + if (g.valid && gi < g.values.size() && gi < g.validMask.size() && + g.validMask[gi]) { + q = quant.toQ(static_cast(g.values[gi])); + } + built.vol.at(x, y, zo) = q; + } + } + } + const double gridMs = nowMs(tGrid); + + // 5) 几何：origin=CGCS2000 世界米(东, 北, 0)，spacing=(cell, cell, dz×stride)。 + const double dz = depthSpacingZ(h) * stride; + built.quant = quant; + built.origin = {originX, originY, 0.0}; + built.spacing = {cellSize, cellSize, dz}; + built.vminPhys = vmin; + built.vmaxPhys = vmax; + + if (metricsOut) { + metricsOut->nx = nx; + metricsOut->ny = ny; + metricsOut->nz = nz; + metricsOut->rtkPoints = static_cast(track.pts.size()); + metricsOut->channels = channels; + metricsOut->cgcsZone = zone; + metricsOut->centralMeridianDeg = centralMeridianDeg; + metricsOut->meanAbsBefore = meanBefore; + metricsOut->meanAbsAfter = meanAfter; + metricsOut->vminPhys = vmin; + metricsOut->vmaxPhys = vmax; + metricsOut->cellSizeM = cellSize; + metricsOut->dz = dz; + metricsOut->originX = originX; + metricsOut->originY = originY; + metricsOut->filledCells = filledCells; + metricsOut->totalCells = static_cast(nx) * ny * nz; + metricsOut->loadMs = loadMs; + metricsOut->pipelineMs = pipelineMs; + metricsOut->trajMs = trajMs; + metricsOut->gridMs = gridMs; + } + + return built; +} + +} // namespace geopro::io::gpr diff --git a/src/io/gpr/Gpr3dvSurveyVolumeBridge.hpp b/src/io/gpr/Gpr3dvSurveyVolumeBridge.hpp new file mode 100644 index 0000000..c893d5f --- /dev/null +++ b/src/io/gpr/Gpr3dvSurveyVolumeBridge.hpp @@ -0,0 +1,71 @@ +#ifndef GEOPRO_IO_GPR_GPR3DVSURVEYVOLUMEBRIDGE_HPP +#define GEOPRO_IO_GPR_GPR3DVSURVEYVOLUMEBRIDGE_HPP + +#include + +#include "core/algo/GprVolumeBuilder.hpp" // geopro::core::BuiltI16 + +namespace geopro::io::gpr { + +// 测绘级(survey-grade)精确坐标逐线建体桥接(Task P8)。 +// +// 与 Gpr3dvVolumeBridge(线局部坐标 X=道/Y=通道/Z=样本，origin≈0)不同： +// 本桥接用 vendored 3DGPRViewer 的精确坐标/轨迹/网格代码(算法零改动)， +// 让单条测线严格按 CGCS2000 大地坐标、逐通道逐道跟 GPS 轨迹建成【世界轴对齐体】， +// 体的 origin/spacing 为真实 CGCS2000 米，多线共享同一参考系 → view-all 直接按 origin +// 摆放即精确就位(跟路的弯，非起点+航向刚体近似)。 +// +// 链路(全程走 P1 原版 + P8 拷入模块，算法零改动)： +// IprhParser::loadImpulseMultiChannel(dir, base, model) // 多通道交错 + 头几何 +// → model.buildVolumeData() // 处理前立方体(before 统计) +// → RadarProcessor::runPipeline(默认链) // 信号处理(11 步选配) +// → processed.buildVolumeData() // 处理后立方体 +// geopro::io::gpr::parseGps(.gps) → 逐 RTK 点经纬(度) +// → CoordinateTransform::wgs84ToCgcs2000 // CGCS2000 平面(北X,东Y) 米 +// → SurveyLine.gpsPositions(QVector3D: x=北, y=东, z=高) +// SurveyGeometry::fromHeaderOffsets(头 CH_X/Y_OFFSETS) // 逐通道天线横偏 +// TrajectoryCalculator::computeTrajectories(line, geom) // 每通道每道 CGCS 世界坐标(跟弯) +// 逐深度 sampleIndex：CScanGridder::build({line}, opts) // 世界坐标 IDW 网格 C-scan +// → 沿深度堆成世界轴对齐规则体(空 cell=kBlank) +// 量化 int16(中点 offset) → BuiltI16，origin=CGCS2000 世界米、spacing=(cell,cell,dz)。 +// +// 轴映射(geopro 体)：X=东向(CGCS easting)、Y=北向(CGCS northing)、Z=深度。 +// (与 CScanGridder 同：gridX=pos.y()东、gridY=pos.x()北，图像与地图屏幕轴对齐。) +struct SurveyBridgeMetrics { + int nx = 0, ny = 0, nz = 0; // 体维度(东 × 北 × 深度) + int rtkPoints = 0; // .gps RTK 点数 + int channels = 0; // 通道数 + int cgcsZone = 0; // CGCS2000 3°带号 + double centralMeridianDeg = 0.0; // 中央经线(度) + double meanAbsBefore = 0.0; // 处理前 traces 平均绝对幅值 + double meanAbsAfter = 0.0; // 处理后 traces 平均绝对幅值 + double vminPhys = 0.0, vmaxPhys = 0.0; + double cellSizeM = 0.0; // 水平格距(米) + double dz = 0.0; // 深度采样间距(米) + double originX = 0.0, originY = 0.0; // CGCS2000 世界 origin(东, 北) 米 + std::int64_t filledCells = 0; // 非空体素数(validMask 命中累计) + std::int64_t totalCells = 0; // 总体素数 + double loadMs = 0.0; // 读 + 建立方体 + double pipelineMs = 0.0; // runPipeline + 再建立方体 + double trajMs = 0.0; // GPS→CGCS + 轨迹 + double gridMs = 0.0; // 逐深度 IDW 网格 + 量化填体 +}; + +// 走 P8 测绘级链建单线世界对齐体。 +// lineDir/linePrefix 同 build-line(如 "D:/Downloads/明星路", "明星路_001")。 +// gpsPath 为该线 .gps 绝对路径(经纬度，geopro parseGps 解析)。 +// coarse(≥1)：省空间档。水平格距 ×coarse、深度每 coarse 个采样取 1(保形)；coarse≤1 全分辨率。 +// cellSizeM/searchRadiusM:基础水平格距/搜索半径(米)，0=用 CScanGridder 默认(0.05/0.50)。 +// metricsOut 非空时回填维度/世界 origin/量化/耗时(供 CLI 报告，不编造)。 +// 失败(加载失败/GPS 无效/体为空) → 抛 std::runtime_error。 +geopro::core::BuiltI16 buildLineVolumeSurvey(const std::string& lineDir, + const std::string& linePrefix, + const std::string& gpsPath, + SurveyBridgeMetrics* metricsOut, + int coarse = 1, + double cellSizeM = 0.0, + double searchRadiusM = 0.0); + +} // namespace geopro::io::gpr + +#endif // GEOPRO_IO_GPR_GPR3DVSURVEYVOLUMEBRIDGE_HPP diff --git a/tools/gpr_poc/main.cpp b/tools/gpr_poc/main.cpp index b61c5db..1532143 100644 --- a/tools/gpr_poc/main.cpp +++ b/tools/gpr_poc/main.cpp @@ -34,6 +34,7 @@ #include "core/model/ColorScale.hpp" #include "core/model/ScalarVolumeI16.hpp" #include "data/store/ChunkedVolumeStore.hpp" +#include "io/gpr/Gpr3dvSurveyVolumeBridge.hpp" #include "io/gpr/Gpr3dvVolumeBridge.hpp" #include "io/gpr/GprSurveyAssembler.hpp" #include "io/gpr/GpsTrack.hpp" @@ -1060,6 +1061,291 @@ int cmdBuildAll(int argc, char** argv) { return 0; } +// ============================================================================ +// build-survey-line / build-survey-all：测绘级精确坐标逐线世界对齐体(Task P8) +// ============================================================================ +// +// 与 build-line(线局部坐标 X=道/Y=通道/Z=样本，origin≈0)不同：本路径走 P8 测绘级桥接 +// (Gpr3dvSurveyVolumeBridge)——用 vendored 3DGPRViewer 的精确坐标/轨迹/网格代码 +// (CoordinateTransform/TrajectoryCalculator/CScanGridder，算法零改动)，让单线严格按 +// CGCS2000 大地坐标、逐通道逐道跟 GPS 轨迹建成【世界轴对齐体】(跟路的弯)。体 meta.origin +// 为真实 CGCS2000 世界米，多线共享同一参考系 → view-survey-all 直接按 origin 摆放即精确就位。 + +// 按线前缀(如 "明星路_001")在 dir 下找匹配的 .gps(尾段 _NNN 相同)。空串=未找到。 +std::string findGpsForPrefix(const std::string& dir, const std::string& prefix) { + const std::size_t us = prefix.find_last_of('_'); + if (us == std::string::npos) return ""; + const std::string num = prefix.substr(us + 1); + std::error_code ec; + for (const auto& e : fs::directory_iterator(dir, ec)) { + 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 s2 = stem.find_last_of('_'); + if (s2 != std::string::npos && stem.substr(s2 + 1) == num) + return e.path().string(); + } + return ""; +} + +// 单线测绘级建体核心：P8 桥接 → 世界对齐量化体 → 落盘 + 金字塔。 +// 异常(加载失败/GPS 无效/体空/退化)由调用方捕获，不在此中断批量。 +LineBuildResult buildOneSurveyLine(const std::string& lineDir, + const std::string& linePrefix, + const std::string& gpsPath, + const std::string& out, int levels, + int coarse, double cellSizeM, + double searchRadiusM) { + LineBuildResult r; + r.prefix = linePrefix; + + std::cout << "[build-survey-line] lineDir=" << lineDir + << " linePrefix=" << linePrefix << " gps=" << gpsPath + << " levels=" << levels << " coarse=" << coarse + << " cellSize=" << cellSizeM << " out=" << out << "\n"; + + Stopwatch swBridge; + geopro::io::gpr::SurveyBridgeMetrics bm; + geopro::core::BuiltI16 built = geopro::io::gpr::buildLineVolumeSurvey( + lineDir, linePrefix, gpsPath, &bm, coarse, cellSizeM, searchRadiusM); + const double bridgeMs = swBridge.elapsedMs(); + + const std::int64_t nx = built.vol.nx(), ny = built.vol.ny(), + nz = built.vol.nz(); + r.nx = nx; + r.ny = ny; + r.nz = nz; + if (nx < 2 || ny < 2 || nz < 2) { + r.ok = false; + r.reason = "世界体维度退化(东×北×深=" + std::to_string(nx) + "x" + + std::to_string(ny) + "x" + std::to_string(nz) + + ")，无法建可看体，跳过"; + std::cerr << "[build-survey-line] " << linePrefix << " " << r.reason << "\n"; + return r; + } + const std::int64_t voxels = nx * ny * nz; + const std::int64_t rawBytes = voxels * 2; + const double fillRate = + bm.totalCells > 0 + ? static_cast(bm.filledCells) / bm.totalCells + : 0.0; + + std::cout << "[build-survey-line] 处理前后平均绝对幅值: " << bm.meanAbsBefore + << " → " << bm.meanAbsAfter << "\n"; + std::cout << "[build-survey-line] RTK点=" << bm.rtkPoints + << " 通道=" << bm.channels << " CGCS带号=" << bm.cgcsZone + << " 中央经线=" << bm.centralMeridianDeg << "°\n"; + std::cout << "[build-survey-line] 世界体维度(东×北×深) = " << nx << " x " << ny + << " x " << nz << " 非空体素=" << bm.filledCells << " (" + << (fillRate * 100.0) << "%)\n"; + std::cout << "[build-survey-line] CGCS2000 世界 origin=(" << std::fixed + << bm.originX << ", " << bm.originY << ") spacing cell=" << bm.cellSizeM + << " dz=" << bm.dz << " (m)\n"; + + fs::create_directories(out); + Stopwatch swWrite; + geopro::data::ChunkedVolumeStore::write(out, built); + const double writeMs = swWrite.elapsedMs(); + + Stopwatch swPyr; + { + geopro::data::ChunkedVolumeStore store(out); + store.buildPyramid(levels); + } + const double pyrMs = swPyr.elapsedMs(); + + const std::int64_t dataBytes = storeDataBytes(out); + r.dataBytes = dataBytes; + const double ratio = + dataBytes > 0 ? static_cast(rawBytes) / dataBytes : 0.0; + const double peak = Probe::peakMemMB(); + + std::cout << "\n=== build-survey-line 指标(测绘级 CGCS2000 世界对齐体)===\n"; + std::cout << "桥接耗时(ms) : " << bridgeMs << " (读 " << bm.loadMs + << " + 处理 " << bm.pipelineMs << " + GPS轨迹 " << bm.trajMs + << " + 网格量化 " << bm.gridMs << ")\n"; + std::cout << "落盘耗时(ms) : " << writeMs << "\n"; + std::cout << "金字塔耗时(ms) : " << pyrMs << "\n"; + std::cout << "世界体维度 : " << nx << " x " << ny << " x " << nz << "\n"; + std::cout << "体素数 : " << voxels << " 非空 " << bm.filledCells + << " (" << (fillRate * 100.0) << "%)\n"; + std::cout << "处理后值域 : [" << bm.vminPhys << ", " << bm.vmaxPhys + << "] 量化 scale=" << built.quant.scale + << " offset=" << built.quant.offset << "\n"; + std::cout << "CGCS2000 origin: (" << std::fixed << bm.originX << ", " + << bm.originY << ") (东, 北) 米\n"; + std::cout << "世界 spacing : cell=" << bm.cellSizeM << " dz=" << bm.dz + << " (m)\n"; + std::cout << "data.bin(B) : " << dataBytes << " (" + << dataBytes / (1024.0 * 1024.0) << " MB) 压缩比 " << ratio + << " x\n"; + std::cout << "峰值内存(MB) : " << peak << "\n"; + + writeMetricLine( + "build-survey-line,prefix=" + linePrefix + + ",coarse=" + std::to_string(coarse) + ",nx=" + std::to_string(nx) + + ",ny=" + std::to_string(ny) + ",nz=" + std::to_string(nz) + + ",voxels=" + std::to_string(voxels) + + ",filled=" + std::to_string(bm.filledCells) + + ",zone=" + std::to_string(bm.cgcsZone) + + ",originX=" + std::to_string(bm.originX) + + ",originY=" + std::to_string(bm.originY) + + ",cell=" + std::to_string(bm.cellSizeM) + ",dz=" + std::to_string(bm.dz) + + ",vmin=" + std::to_string(bm.vminPhys) + + ",vmax=" + std::to_string(bm.vmaxPhys) + + ",dataB=" + std::to_string(dataBytes) + + ",bridgeMs=" + std::to_string(bridgeMs) + + ",peakMB=" + std::to_string(peak)); + + r.ok = true; + return r; +} + +int cmdBuildSurveyLine(int argc, char** argv) { + const Args a = parseArgs(argc, argv, 2); + if (a.positional.size() < 2) { + std::cerr << "用法: gpr_poc build-survey-line " + "--out [--levels 3] [--coarse F] [--cell 0.05] " + "[--radius 0.5] [--gps ]\n" + "例: gpr_poc build-survey-line \"D:/Downloads/明星路\" " + "明星路_001 --out tmp/survey001 --levels 3 --coarse 4\n"; + return 2; + } + const std::string lineDir = a.positional[0]; + const std::string linePrefix = a.positional[1]; + const int levels = std::stoi(a.get("levels", "3")); + const int coarse = std::stoi(a.get("coarse", "1")); + const double cellSizeM = std::stod(a.get("cell", "0")); + const double radiusM = std::stod(a.get("radius", "0")); + std::string gps = a.get("gps", ""); + if (gps.empty()) gps = findGpsForPrefix(lineDir, linePrefix); + const std::string out = a.get( + "out", (fs::temp_directory_path() / "gpr_store_survey_line").string()); + + if (gps.empty()) { + std::cerr << "[build-survey-line] 失败: 未找到 " << linePrefix + << " 的 .gps(可用 --gps 指定)\n"; + return 1; + } + try { + const LineBuildResult r = buildOneSurveyLine( + lineDir, linePrefix, gps, out, levels, coarse, cellSizeM, radiusM); + if (!r.ok) { + std::cerr << "[build-survey-line] 跳过: " << r.reason << "\n"; + return 1; + } + return 0; + } catch (const std::exception& e) { + std::cerr << "[build-survey-line] 失败(" << linePrefix << "): " << e.what() + << "\n"; + return 1; + } +} + +// build-survey-all：发现所有测线 → 逐条走 P8 测绘级路径建世界对齐体到 baseDir//。 +// 各线 .gps 自动按尾段 _NNN 匹配；缺 .gps 的线跳过报因。磁盘守护同 build-all。 +int cmdBuildSurveyAll(int argc, char** argv) { + const Args a = parseArgs(argc, argv, 2); + if (a.positional.empty() || !a.kv.count("outBase")) { + std::cerr << "用法: gpr_poc build-survey-all --outBase " + "[--levels 3] [--coarse F] [--cell 0.05] [--radius 0.5] " + "[--minFreeGB 3]\n" + "例: gpr_poc build-survey-all \"D:/Downloads/明星路\" " + "--outBase tmp/survey_all --levels 3 --coarse 4\n"; + return 2; + } + const std::string lineDir = a.positional[0]; + const std::string outBase = a.get("outBase", ""); + const int levels = std::stoi(a.get("levels", "3")); + const int coarse = std::stoi(a.get("coarse", "1")); + const double cellSizeM = std::stod(a.get("cell", "0")); + const double radiusM = std::stod(a.get("radius", "0")); + const double minFreeGB = std::stod(a.get("minFreeGB", "3")); + + // 发现所有测线前缀(扫 *_A.iprh，截 _A 得前缀)，与 build-all 同口径。 + std::set prefixSet; + std::error_code ec; + for (const auto& e : fs::directory_iterator(lineDir, ec)) { + if (!e.is_regular_file()) continue; + if (e.path().extension().string() != ".iprh") continue; + const std::string name = e.path().filename().string(); + const std::size_t pos = name.rfind("_A"); + if (pos == std::string::npos) continue; + std::size_t d = pos + 2; + while (d < name.size() && std::isdigit(static_cast(name[d]))) + ++d; + if (d == pos + 2) continue; + prefixSet.insert(name.substr(0, pos)); + } + if (prefixSet.empty()) { + std::cerr << "[build-survey-all] 未在 " << lineDir + << " 发现任何测线(*_A.iprh)\n"; + return 1; + } + std::vector prefixes(prefixSet.begin(), prefixSet.end()); + std::cout << "[build-survey-all] 发现 " << prefixes.size() + << " 条测线，outBase=" << outBase << " levels=" << levels + << " coarse=" << coarse << " minFreeGB=" << minFreeGB << "\n"; + + fs::create_directories(outBase, ec); + + std::vector results; + bool stoppedByDisk = false; + for (const std::string& prefix : prefixes) { + const double freeGB = freeSpaceGB(outBase); + std::cout << "\n[build-survey-all] --- " << prefix << " --- 剩余磁盘 " + << freeGB << " GB\n"; + if (freeGB >= 0.0 && freeGB < minFreeGB) { + std::cerr << "[build-survey-all] 磁盘守护触发: 剩余 " << freeGB << " GB < " + << minFreeGB << " GB，停止，未建 " << prefix << " 及其后。\n"; + stoppedByDisk = true; + break; + } + const std::string out = (fs::path(outBase) / prefix).string(); + LineBuildResult r; + r.prefix = prefix; + const std::string gps = findGpsForPrefix(lineDir, prefix); + if (gps.empty()) { + r.ok = false; + r.reason = "缺 .gps，跳过"; + std::cerr << "[build-survey-all] " << prefix << " 缺 .gps，跳过\n"; + results.push_back(r); + continue; + } + try { + r = buildOneSurveyLine(lineDir, prefix, gps, out, levels, coarse, + cellSizeM, radiusM); + } catch (const std::exception& e) { + r.ok = false; + r.reason = std::string("异常: ") + e.what(); + std::cerr << "[build-survey-all] " << prefix << " 失败: " << e.what() + << "(跳过，继续)\n"; + } + results.push_back(r); + } + + std::cout << "\n=== build-survey-all 汇总(测绘级 CGCS2000 世界对齐体) ===\n"; + int okCount = 0; + std::int64_t totalBytes = 0; + for (const auto& r : results) { + if (r.ok) { + ++okCount; + totalBytes += r.dataBytes; + std::cout << " [OK] " << r.prefix << " 维度(东×北×深)=" << r.nx << "x" + << r.ny << "x" << r.nz << " data=" + << r.dataBytes / (1024.0 * 1024.0) << " MB\n"; + } else { + std::cout << " [跳过] " << r.prefix << " 原因=" << r.reason << "\n"; + } + } + std::cout << "成功 " << okCount << "/" << prefixes.size() + << " 条，合计 data=" << totalBytes / (1024.0 * 1024.0 * 1024.0) + << " GB，剩余磁盘 " << freeSpaceGB(outBase) << " GB\n"; + if (stoppedByDisk) + std::cout << "注意: 因磁盘守护提前停止，部分测线未建。\n"; + return 0; +} + int cmdLoad(int argc, char** argv) { const Args a = parseArgs(argc, argv, 2); if (a.positional.empty()) { @@ -4013,6 +4299,243 @@ int cmdViewAll(int argc, char** argv) { return 0; } +// ============================================================================ +// view-survey-all：测绘级世界对齐体直接按 CGCS2000 origin 摆放渲染(Task P8) +// ============================================================================ +// +// 与 view-all(线局部体 + .gps 起点+航向刚体近似摆放)不同：本命令的体已是测绘级世界对齐 +// 体(build-survey-all 产出，meta.origin=CGCS2000 世界米，轴 X=东/Y=北/Z=深)。故【无需 .gps、 +// 无需航向旋转】——直接按各体 meta.origin 平移摆进同一世界框即精确就位(跟路的弯)。 +// 取公共参考原点(全体 origin 最小东/北)平移到近零局部框，保 VTK 数值稳定且相对位置严格不变。 + +// 由测绘级世界对齐体 + 公共参考原点 + Z 夸张 → vtkVolume(平移到世界位，无旋转)。 +vtkSmartPointer makeSurveyPlacedVolume( + vtkImageData* img, const geopro::data::StoreMeta& m, double refX, + double refY, double exagg, double& vminOut, double& vmaxOut) { + const GalleryVariant& v = kViewDefaultVariant; + + double vmin = m.vminPhys, vmax = m.vmaxPhys; + const ScalarPercentiles pc = + sampleScalarPercentiles(img, m.quant, 0.02, 0.98); + if (pc.samples > 0) { + vmin = pc.lo; + vmax = pc.hi; + } + vminOut = vmin; + vmaxOut = vmax; + const geopro::core::ColorScale cs = pickColor(v.color, vmin, vmax); + + GradStats gs; + if (v.useGradientOpacity) gs = sampleGradientMagnitude(img); + vtkSmartPointer prop = makeVariantProperty( + v, m.quant, cs, vmin, vmax, v.maxOpacity, + v.useGradientOpacity ? &gs : nullptr); + + // img 的 origin 是 meta.origin(CGCS 世界米，东向 ~4e7 量级)。直接把 image origin 减去公共 + // 参考原点落到近零局部框——避免 VTK 内部 float32 在 ~4e7 绝对坐标上丢精度(±1m)；相对位置严格不变。 + double iorg[3]; + img->GetOrigin(iorg); + img->SetOrigin(iorg[0] - refX, iorg[1] - refY, iorg[2]); + + vtkNew mapper; + mapper->SetInputData(img); + mapper->SetRequestedRenderMode(vtkSmartVolumeMapper::GPURenderMode); + mapper->SetAutoAdjustSampleDistances(0); + mapper->SetInteractiveAdjustSampleDistances(0); + + auto volume = vtkSmartPointer::New(); + volume->SetMapper(mapper); + volume->SetProperty(prop); + + // 体已落近零世界框；仅 Z 夸张(局部深度)。 + auto xf = vtkSmartPointer::New(); + xf->PostMultiply(); + xf->Scale(1.0, 1.0, exagg); + volume->SetUserTransform(xf); + return volume; +} + +int cmdViewSurveyAll(int argc, char** argv) { + const Args a = parseArgs(argc, argv, 2); + if (a.positional.empty()) { + std::cerr << "用法: gpr_poc view-survey-all [--preview] " + "[--exagg 8] [--level 1] [--shotDir ]\n" + "例: gpr_poc view-survey-all tmp/survey_all --preview " + "--exagg 8\n"; + return 2; + } + const std::string storesDir = a.positional[0]; + const double exagg = std::stod(a.get("exagg", "8")); + const int level = std::stoi(a.get("level", "1")); + const bool preview = a.kv.count("preview") > 0; + const std::string shotDir = a.get("shotDir", storesDir); + + std::cout << "[view-survey-all] storesDir=" << storesDir << " exagg=" << exagg + << " level=" << level + << (preview ? " [PREVIEW 离屏俯视+斜视出图]" : " [真窗口可交互]") + << "\n"; + + std::cout << "[view-survey-all] 离屏闸门复检...\n"; + if (cmdOffscreenSmoke() != 0) { + std::cout << "[view-survey-all] 闸门失败，中止。\n"; + return 1; + } + + // 1) 发现体目录(含 meta.json)，按名排序。 + std::vector storeNames; + for (const auto& e : fs::directory_iterator(storesDir)) { + if (!e.is_directory()) continue; + if (!fs::exists(e.path() / "meta.json")) continue; + storeNames.push_back(e.path().filename().string()); + } + std::sort(storeNames.begin(), storeNames.end()); + std::cout << "[view-survey-all] 发现体目录数=" << storeNames.size() << "\n"; + if (storeNames.empty()) { + std::cerr << "[view-survey-all] 错误: 未发现任何含 meta.json 的体目录\n"; + return 1; + } + + // 2) 逐体读 meta，先定公共参考原点(全体 origin 最小东/北)。 + struct SurveyPlaced { + std::string name; + vtkSmartPointer img; + geopro::data::StoreMeta meta; + }; + std::vector placed; + double refX = std::numeric_limits::infinity(); + double refY = std::numeric_limits::infinity(); + for (const std::string& nm : storeNames) { + const std::string storePath = (fs::path(storesDir) / nm).string(); + geopro::data::ChunkedVolumeStore store(storePath); + const geopro::data::StoreMeta m = store.meta(); + refX = std::min(refX, m.origin[0]); + refY = std::min(refY, m.origin[1]); + const int lv = std::max(0, std::min(level, store.levels() - 1)); + vtkSmartPointer img = buildLevelImage(store, lv, m); + std::cout << "[view-survey-all] " << nm << " level=" << lv << " 维度=" + << img->GetDimensions()[0] << "x" << img->GetDimensions()[1] << "x" + << img->GetDimensions()[2] << " CGCS origin=(" << std::fixed + << m.origin[0] << ", " << m.origin[1] << ")\n"; + placed.push_back({nm, img, m}); + } + if (placed.empty()) { + std::cerr << "[view-survey-all] 错误: 无可用体\n"; + return 1; + } + std::cout << "[view-survey-all] 公共参考原点(最小东/北)=(" << std::fixed << refX + << ", " << refY << ") (共 " << placed.size() << " 体)\n"; + + // 3) 同一 renderer 加全部世界对齐体(按 origin-ref 摆放，无旋转)。 + const int winW = 1400, winH = 900; + auto rw = preview ? makeOffscreenWindow(winW, winH) + : vtkSmartPointer::New(); + if (!preview) rw->SetSize(winW, winH); + vtkNew ren; + ren->SetBackground(kViewDefaultVariant.bg[0], kViewDefaultVariant.bg[1], + kViewDefaultVariant.bg[2]); + rw->AddRenderer(ren); + + auto capWin = vtkSmartPointer::New(); + vtkOutputWindow::SetInstance(capWin); + + int added = 0; + for (const SurveyPlaced& sp : placed) { + double vmin = 0, vmax = 0; + vtkSmartPointer vol = makeSurveyPlacedVolume( + sp.img.Get(), sp.meta, refX, refY, exagg, vmin, vmax); + ren->AddVolume(vol); + ++added; + } + std::cout << "[view-survey-all] 已加入场景体数=" << added << "\n"; + + ren->ResetCamera(); + rw->Render(); + if (capWin->textureError()) { + std::cerr << "[view-survey-all] 警告: 检测到 3D 纹理维度错误(某体超 GL 上限)，" + "可增大 --level 取更粗层。\n"; + } + + if (preview) { + fs::create_directories(shotDir); + { // 俯视(top)：看 20 条按真实 CGCS 位置铺成测区(含路的弯)。 + ren->ResetCamera(); + vtkCamera* cam = ren->GetActiveCamera(); + double fp[3]; + cam->GetFocalPoint(fp); + const double* b = ren->ComputeVisiblePropBounds(); + const double span = std::max({b[1] - b[0], b[3] - b[2], b[5] - b[4]}); + cam->SetPosition(fp[0], fp[1], fp[2] + span * 2.0); + cam->SetFocalPoint(fp[0], fp[1], fp[2]); + cam->SetViewUp(0, 1, 0); + ren->ResetCameraClippingRange(); + rw->Render(); + const std::string p = + (fs::path(shotDir) / "view-survey-all-top.png").string(); + savePng(rw.Get(), p); + std::cout << "[view-survey-all] 俯视图存: " << p << "\n"; + } + { // 斜视(oblique)。 + ren->ResetCamera(); + vtkCamera* cam = ren->GetActiveCamera(); + cam->Elevation(35.0); + cam->Azimuth(30.0); + ren->ResetCameraClippingRange(); + rw->Render(); + const std::string p = + (fs::path(shotDir) / "view-survey-all-oblique.png").string(); + savePng(rw.Get(), p); + std::cout << "[view-survey-all] 斜视图存: " << p << "\n"; + } + + rw->Render(); + Stopwatch sw; + const int frames = 60; + for (int f = 0; f < frames; ++f) rw->Render(); + const double ms = sw.elapsedMs(); + const double fps = ms > 0 ? frames * 1000.0 / ms : 0.0; + const vtkIdType nb = countNonBlackPixels(rw.Get(), winW, winH); + + std::cout << "\n=== view-survey-all --preview 测区全貌(测绘级世界对齐体)===\n"; + std::cout << "参与体数 : " << placed.size() << "\n"; + std::cout << "level : " << level << "\n"; + std::cout << "exagg(Z) : " << exagg << "\n"; + std::cout << "公共参考原点 : (" << std::fixed << refX << ", " << refY + << ") CGCS2000 米\n"; + std::cout << "非黑像素 : " << nb << " / " << (winW * winH) << "\n"; + std::cout << "fps(" << frames << "帧连渲) : " << fps << "\n"; + std::cout << "俯视图 : " + << (fs::path(shotDir) / "view-survey-all-top.png").string() + << "\n"; + std::cout << "斜视图 : " + << (fs::path(shotDir) / "view-survey-all-oblique.png").string() + << "\n"; + + writeMetricLine("view-survey-all,bodies=" + std::to_string(placed.size()) + + ",level=" + std::to_string(level) + + ",exagg=" + std::to_string(exagg) + + ",nonBlack=" + std::to_string(nb) + + ",fps=" + std::to_string(fps)); + return nb > 0 ? 0 : 1; + } + + rw->SetWindowName("gpr_poc view-survey-all —— 测绘级 CGCS2000 世界对齐体"); + vtkNew iren; + iren->SetRenderWindow(rw); + vtkNew style; + iren->SetInteractorStyle(style); + ren->ResetCamera(); + vtkCamera* cam = ren->GetActiveCamera(); + cam->Elevation(35.0); + cam->Azimuth(30.0); + ren->ResetCameraClippingRange(); + std::cout << "[view-survey-all] 打开真窗口。左键旋转 / 滚轮缩放 / q 退出。\n"; + iren->Initialize(); + rw->Render(); + iren->Start(); + std::cout << "[view-survey-all] 窗口关闭，退出。\n"; + return 0; +} + int cmdView(int argc, char** argv) { const Args a = parseArgs(argc, argv, 2); if (a.positional.empty()) { @@ -4886,6 +5409,12 @@ void usage() { "[--levels 3] [--coarse F]\n" " gpr_poc build-all --outBase " "[--levels 3] [--coarse F] [--minFreeGB 3]\n" + " gpr_poc build-survey-line --out " + " [--levels 3] [--coarse F] [--cell 0.05] " + "[--radius 0.5] [--gps ]\n" + " gpr_poc build-survey-all --outBase " + "[--levels 3] [--coarse F] [--cell 0.05] [--radius 0.5] " + "[--minFreeGB 3]\n" " gpr_poc load \n" " gpr_poc selftest\n" " gpr_poc offscreen-smoke\n" @@ -4902,6 +5431,8 @@ void usage() { "[--frames 90]\n" " gpr_poc view-all [--preview] " "[--exagg 8] [--level 1] [--spread M] [--shotDir ]\n" + " gpr_poc view-survey-all [--preview] [--exagg 8] " + "[--level 1] [--shotDir ]\n" " gpr_poc polish [--exagg 8] [--frames 90] " "[--localBricks 4]\n"; } @@ -4924,6 +5455,8 @@ int main(int argc, char** argv) { if (cmd == "build-geo") return cmdBuildGeo(argc, argv); if (cmd == "build-line") return cmdBuildLine(argc, argv); if (cmd == "build-all") return cmdBuildAll(argc, argv); + if (cmd == "build-survey-line") return cmdBuildSurveyLine(argc, argv); + if (cmd == "build-survey-all") return cmdBuildSurveyAll(argc, argv); if (cmd == "load") return cmdLoad(argc, argv); if (cmd == "selftest") return cmdSelftest(); if (cmd == "offscreen-smoke") return cmdOffscreenSmoke(); @@ -4936,6 +5469,7 @@ int main(int argc, char** argv) { if (cmd == "fps-budget") return cmdFpsBudget(argc, argv); if (cmd == "view") return cmdView(argc, argv); if (cmd == "view-all") return cmdViewAll(argc, argv); + if (cmd == "view-survey-all") return cmdViewSurveyAll(argc, argv); if (cmd == "polish") return cmdPolish(argc, argv); } catch (const std::exception& e) { std::cerr << "错误: " << e.what() << "\n";