feat(gpr3dv): 移植精确坐标/轨迹/世界网格(CGCS2000)+测绘级逐线世界对齐建体
复制 CoordinateTransform/TrajectoryCalculator/CScanGridder/PosParser(逐字节一致)进 external/gpr3dviewer;新增 Gpr3dvSurveyVolumeBridge 按 CGCS2000+逐道GPS轨迹建世界对齐体; gpr_poc 加 build-survey-all/view-survey-all(各体自带世界origin,精确就位跟GPS弯)。
This commit is contained in:
gaozheng
parent
ccd3040a67
commit
f713c15366
|
|
@ -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
|
||||
)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,336 @@
|
|||
#include "CScanGridder.h"
|
||||
|
||||
#include <QHash>
|
||||
#include <QPair>
|
||||
#include <QtGlobal>
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
|
||||
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<float> gaussianKernel(double sigma)
|
||||
{
|
||||
sigma = std::max(0.1, sigma);
|
||||
const int radius = std::max(1, static_cast<int>(std::ceil(sigma * 3.0)));
|
||||
QVector<float> 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<float>(v);
|
||||
sum += v;
|
||||
}
|
||||
if (sum > 0.0) {
|
||||
for (float &v : kernel) v = static_cast<float>(v / sum);
|
||||
}
|
||||
return kernel;
|
||||
}
|
||||
|
||||
static void smoothMasked(CScanGridResult &grid, double sigma)
|
||||
{
|
||||
if (!grid.valid || grid.width <= 0 || grid.height <= 0) return;
|
||||
|
||||
const QVector<float> kernel = gaussianKernel(sigma);
|
||||
const int radius = kernel.size() / 2;
|
||||
const int count = grid.width * grid.height;
|
||||
|
||||
QVector<float> tmpValues(count, 0.0f);
|
||||
QVector<float> 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<float>(sum / weight);
|
||||
tmpWeights[outIdx] = static_cast<float>(weight);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
QVector<float> outValues(count, 0.0f);
|
||||
QVector<unsigned char> 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<float>(sum / weight);
|
||||
outMask[outIdx] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
grid.values = std::move(outValues);
|
||||
grid.validMask = std::move(outMask);
|
||||
}
|
||||
|
||||
static float percentile(QVector<float> 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<int>(std::floor(pos));
|
||||
const int hi = static_cast<int>(std::ceil(pos));
|
||||
if (lo == hi) return values[lo];
|
||||
const double t = pos - lo;
|
||||
return static_cast<float>(values[lo] * (1.0 - t) + values[hi] * t);
|
||||
}
|
||||
|
||||
static void updateDisplayRange(CScanGridResult &grid, double lowPercent, double highPercent)
|
||||
{
|
||||
QVector<float> 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<SurveyLine> &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<CScanSamplePoint> samples;
|
||||
double minX = std::numeric_limits<double>::max();
|
||||
double minY = std::numeric_limits<double>::max();
|
||||
double maxX = std::numeric_limits<double>::lowest();
|
||||
double maxY = std::numeric_limits<double>::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<int>(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<QVector3D> &traj = line.channelTrajectories[c];
|
||||
const int n = std::min(traceCount, static_cast<int>(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<float>(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<double>(pos.y()));
|
||||
minY = std::min(minY, static_cast<double>(pos.x()));
|
||||
maxX = std::max(maxX, static_cast<double>(pos.y()));
|
||||
maxY = std::max(maxY, static_cast<double>(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<int>(std::ceil((maxX - minX) / cellSize)) + 1;
|
||||
int height = static_cast<int>(std::ceil((maxY - minY) / cellSize)) + 1;
|
||||
while ((width > options.maxGridWidth || height > options.maxGridHeight
|
||||
|| static_cast<qint64>(width) * height > options.maxGridCells) && cellSize < 10.0) {
|
||||
cellSize *= 1.5;
|
||||
width = static_cast<int>(std::ceil((maxX - minX) / cellSize)) + 1;
|
||||
height = static_cast<int>(std::ceil((maxY - minY) / cellSize)) + 1;
|
||||
}
|
||||
|
||||
if (width <= 0 || height <= 0 || width > options.maxGridWidth || height > options.maxGridHeight
|
||||
|| static_cast<qint64>(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<BinKey, QVector<int>> bins;
|
||||
bins.reserve(samples.size());
|
||||
for (int i = 0; i < samples.size(); ++i) {
|
||||
const CScanSamplePoint &p = samples[i];
|
||||
const BinKey key{static_cast<int>(std::floor((p.x - minX) / binSize)),
|
||||
static_cast<int>(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<int>(std::ceil(options.searchRadiusM / binSize)));
|
||||
|
||||
for (int y = 0; y < height; ++y) {
|
||||
const double cy = minY + y * cellSize;
|
||||
const int by = static_cast<int>(std::floor((cy - minY) / binSize));
|
||||
for (int x = 0; x < width; ++x) {
|
||||
const double cx = minX + x * cellSize;
|
||||
const int bx = static_cast<int>(std::floor((cx - minX) / binSize));
|
||||
|
||||
double sum = 0.0;
|
||||
double weightSum = 0.0;
|
||||
double bestDistSq = std::numeric_limits<double>::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<float>(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;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
#ifndef CSCANGRIDDER_H
|
||||
#define CSCANGRIDDER_H
|
||||
|
||||
#include "GPRDataModel.h"
|
||||
|
||||
#include <QVector>
|
||||
|
||||
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<float> values;
|
||||
QVector<unsigned char> validMask;
|
||||
QVector<int> nearestLine;
|
||||
QVector<int> nearestChannel;
|
||||
QVector<int> nearestTrace;
|
||||
float displayMin = 0.0f;
|
||||
float displayMax = 1.0f;
|
||||
};
|
||||
|
||||
class CScanGridder
|
||||
{
|
||||
public:
|
||||
static CScanGridResult build(const QVector<SurveyLine> &surveyLines, CScanGridOptions options);
|
||||
};
|
||||
|
||||
#endif // CSCANGRIDDER_H
|
||||
|
|
@ -0,0 +1,219 @@
|
|||
#include "CoordinateTransform.h"
|
||||
#include <cmath>
|
||||
|
||||
// 静态辅助：预计算的子午线弧长系数
|
||||
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<int>(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;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
#ifndef COORDINATETRANSFORM_H
|
||||
#define COORDINATETRANSFORM_H
|
||||
|
||||
#include <cmath>
|
||||
|
||||
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
|
||||
|
|
@ -0,0 +1,73 @@
|
|||
#include "PosParser.h"
|
||||
#include <QFile>
|
||||
#include <QTextStream>
|
||||
#include <QDebug>
|
||||
#include <QRegularExpression>
|
||||
|
||||
bool PosParser::loadPosFile(const QString &posFilePath, QVector<QVector3D> &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<float>(x), static_cast<float>(y), static_cast<float>(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;
|
||||
}
|
||||
|
|
@ -0,0 +1,17 @@
|
|||
#ifndef POSPARSER_H
|
||||
#define POSPARSER_H
|
||||
|
||||
#include <QString>
|
||||
#include <QVector>
|
||||
#include <QVector3D>
|
||||
|
||||
class PosParser {
|
||||
public:
|
||||
// 加载 .pos 文件（4列：道号 X Y Z）
|
||||
static bool loadPosFile(const QString &posFilePath, QVector<QVector3D> &outPositions);
|
||||
|
||||
// 加载 center.ccc 项目中心坐标（2列：纬度 经度）
|
||||
static bool loadCenterCcc(const QString &cccFilePath, double &outLat, double &outLon);
|
||||
};
|
||||
|
||||
#endif // POSPARSER_H
|
||||
|
|
@ -0,0 +1,242 @@
|
|||
#include "TrajectoryCalculator.h"
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
|
||||
namespace {
|
||||
|
||||
double planarDistance(const QVector3D &a, const QVector3D &b)
|
||||
{
|
||||
const double dx = static_cast<double>(a.x()) - b.x();
|
||||
const double dy = static_cast<double>(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<float>(1.0 - t) + b * static_cast<float>(t);
|
||||
}
|
||||
|
||||
QVector3D catmullRom(const QVector3D &p0, const QVector3D &p1, const QVector3D &p2, const QVector3D &p3, double t)
|
||||
{
|
||||
const float tf = static_cast<float>(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<QVector3D> &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<QVector3D> &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<double> 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<float>(chX),
|
||||
static_cast<float>(chY),
|
||||
static_cast<float>(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<QVector3D> TrajectoryCalculator::resampleTrajectoryLinear(const QVector<QVector3D> &input, int targetCount)
|
||||
{
|
||||
QVector<QVector3D> 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<double>(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<int>(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<QVector3D> TrajectoryCalculator::resampleTrajectorySpline(const QVector<QVector3D> &input, int targetCount)
|
||||
{
|
||||
if (input.size() < 4) return resampleTrajectoryLinear(input, targetCount);
|
||||
QVector<QVector3D> output;
|
||||
if (targetCount <= 0) return output;
|
||||
if (targetCount == 1) return QVector<QVector3D>{input.first()};
|
||||
|
||||
output.reserve(targetCount);
|
||||
const double scale = static_cast<double>(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<int>(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<QVector3D> &input,
|
||||
const TrajectoryFilterOptions &options,
|
||||
int targetCount)
|
||||
{
|
||||
TrajectoryFilterResult result;
|
||||
const int n = input.size();
|
||||
result.keepMask = QVector<bool>(n, true);
|
||||
result.distanceOutlierMask = QVector<bool>(n, false);
|
||||
result.speedOutlierMask = QVector<bool>(n, false);
|
||||
result.angleOutlierMask = QVector<bool>(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<QVector3D> 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;
|
||||
}
|
||||
|
|
@ -0,0 +1,49 @@
|
|||
#ifndef TRAJECTORYCALCULATOR_H
|
||||
#define TRAJECTORYCALCULATOR_H
|
||||
|
||||
#include "GPRDataModel.h"
|
||||
#include "SurveyGeometry.h"
|
||||
#include <QVector3D>
|
||||
#include <QStringList>
|
||||
|
||||
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<QVector3D> outputPositions;
|
||||
QVector<bool> keepMask;
|
||||
QVector<bool> distanceOutlierMask;
|
||||
QVector<bool> speedOutlierMask;
|
||||
QVector<bool> angleOutlierMask;
|
||||
QStringList warnings;
|
||||
};
|
||||
|
||||
class TrajectoryCalculator {
|
||||
public:
|
||||
// 计算第 traceIdx 个 RTK 点的行进方向角（弧度）
|
||||
// 返回 theta_y = atan2(dY, dX)，其中 X=北向, Y=东向
|
||||
static double headingAt(int traceIdx, const QVector<QVector3D> &gpsPositions);
|
||||
|
||||
// 根据 RTK 轨迹和天线几何参数，计算所有通道的绝对坐标
|
||||
// 结果写入 line.channelTrajectories 和 line.data.traces[].position
|
||||
static bool computeTrajectories(SurveyLine &line, const SurveyGeometry &geom);
|
||||
|
||||
static QVector<QVector3D> resampleTrajectoryLinear(const QVector<QVector3D> &input, int targetCount);
|
||||
static QVector<QVector3D> resampleTrajectorySpline(const QVector<QVector3D> &input, int targetCount);
|
||||
static TrajectoryFilterResult filterAndInterpolateTrajectory(const QVector<QVector3D> &input,
|
||||
const TrajectoryFilterOptions &options,
|
||||
int targetCount);
|
||||
};
|
||||
|
||||
#endif // TRAJECTORYCALCULATOR_H
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,270 @@
|
|||
#include "io/gpr/Gpr3dvSurveyVolumeBridge.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <chrono>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <stdexcept>
|
||||
|
||||
#include <QString>
|
||||
#include <QVector>
|
||||
#include <QVector3D>
|
||||
|
||||
#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<long double>(v < 0 ? -v : v);
|
||||
++count;
|
||||
}
|
||||
}
|
||||
return count > 0 ? static_cast<double>(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<double, std::milli>(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<int>(std::lround(lon0Deg / 3.0));
|
||||
const double centralMeridianDeg = CoordinateTransform::centralMeridianFromZone(zone);
|
||||
geom.cgcsZone = zone;
|
||||
geom.centralMeridianDeg = centralMeridianDeg;
|
||||
|
||||
QVector<QVector3D> gpsCgcs;
|
||||
gpsCgcs.reserve(static_cast<int>(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<float>(cx), static_cast<float>(cy),
|
||||
static_cast<float>(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<SurveyLine> 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<double>::infinity();
|
||||
double vmax = -std::numeric_limits<double>::infinity();
|
||||
// 缓存每层网格(避免二次计算)：网格层数 = nzOut，每层 nx*ny float，明星路约
|
||||
// 数百×数千×数百 → 单线峰值可控(逐线建，建完即释放)。
|
||||
std::vector<CScanGridResult> layers;
|
||||
layers.reserve(static_cast<std::size_t>(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<double>(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<std::size_t>(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<double>(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<int>(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<std::int64_t>(nx) * ny * nz;
|
||||
metricsOut->loadMs = loadMs;
|
||||
metricsOut->pipelineMs = pipelineMs;
|
||||
metricsOut->trajMs = trajMs;
|
||||
metricsOut->gridMs = gridMs;
|
||||
}
|
||||
|
||||
return built;
|
||||
}
|
||||
|
||||
} // namespace geopro::io::gpr
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
#ifndef GEOPRO_IO_GPR_GPR3DVSURVEYVOLUMEBRIDGE_HPP
|
||||
#define GEOPRO_IO_GPR_GPR3DVSURVEYVOLUMEBRIDGE_HPP
|
||||
|
||||
#include <string>
|
||||
|
||||
#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
|
||||
|
|
@ -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<double>(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<double>(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 <lineDir> <linePrefix> "
|
||||
"--out <storeDir> [--levels 3] [--coarse F] [--cell 0.05] "
|
||||
"[--radius 0.5] [--gps <path>]\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/<lineName>/。
|
||||
// 各线 .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 <lineDir> --outBase <baseDir> "
|
||||
"[--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<NN>.iprh，截 _A<NN> 得前缀)，与 build-all 同口径。
|
||||
std::set<std::string> 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<unsigned char>(name[d])))
|
||||
++d;
|
||||
if (d == pos + 2) continue;
|
||||
prefixSet.insert(name.substr(0, pos));
|
||||
}
|
||||
if (prefixSet.empty()) {
|
||||
std::cerr << "[build-survey-all] 未在 " << lineDir
|
||||
<< " 发现任何测线(*_A<NN>.iprh)\n";
|
||||
return 1;
|
||||
}
|
||||
std::vector<std::string> 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<LineBuildResult> 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<vtkVolume> 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<vtkVolumeProperty> 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<vtkSmartVolumeMapper> mapper;
|
||||
mapper->SetInputData(img);
|
||||
mapper->SetRequestedRenderMode(vtkSmartVolumeMapper::GPURenderMode);
|
||||
mapper->SetAutoAdjustSampleDistances(0);
|
||||
mapper->SetInteractiveAdjustSampleDistances(0);
|
||||
|
||||
auto volume = vtkSmartPointer<vtkVolume>::New();
|
||||
volume->SetMapper(mapper);
|
||||
volume->SetProperty(prop);
|
||||
|
||||
// 体已落近零世界框；仅 Z 夸张(局部深度)。
|
||||
auto xf = vtkSmartPointer<vtkTransform>::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 <storesDir> [--preview] "
|
||||
"[--exagg 8] [--level 1] [--shotDir <dir>]\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<std::string> 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<vtkImageData> img;
|
||||
geopro::data::StoreMeta meta;
|
||||
};
|
||||
std::vector<SurveyPlaced> placed;
|
||||
double refX = std::numeric_limits<double>::infinity();
|
||||
double refY = std::numeric_limits<double>::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<vtkImageData> 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<vtkRenderWindow>::New();
|
||||
if (!preview) rw->SetSize(winW, winH);
|
||||
vtkNew<vtkRenderer> ren;
|
||||
ren->SetBackground(kViewDefaultVariant.bg[0], kViewDefaultVariant.bg[1],
|
||||
kViewDefaultVariant.bg[2]);
|
||||
rw->AddRenderer(ren);
|
||||
|
||||
auto capWin = vtkSmartPointer<CapturingOutputWindow>::New();
|
||||
vtkOutputWindow::SetInstance(capWin);
|
||||
|
||||
int added = 0;
|
||||
for (const SurveyPlaced& sp : placed) {
|
||||
double vmin = 0, vmax = 0;
|
||||
vtkSmartPointer<vtkVolume> 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<vtkRenderWindowInteractor> iren;
|
||||
iren->SetRenderWindow(rw);
|
||||
vtkNew<vtkInteractorStyleTrackballCamera> 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 <lineDir> --outBase <baseDir> "
|
||||
"[--levels 3] [--coarse F] [--minFreeGB 3]\n"
|
||||
" gpr_poc build-survey-line <lineDir> <linePrefix> --out "
|
||||
"<storeDir> [--levels 3] [--coarse F] [--cell 0.05] "
|
||||
"[--radius 0.5] [--gps <path>]\n"
|
||||
" gpr_poc build-survey-all <lineDir> --outBase <baseDir> "
|
||||
"[--levels 3] [--coarse F] [--cell 0.05] [--radius 0.5] "
|
||||
"[--minFreeGB 3]\n"
|
||||
" gpr_poc load <storeDir>\n"
|
||||
" gpr_poc selftest\n"
|
||||
" gpr_poc offscreen-smoke\n"
|
||||
|
|
@ -4902,6 +5431,8 @@ void usage() {
|
|||
"[--frames 90]\n"
|
||||
" gpr_poc view-all <storesDir> <gpsDir> [--preview] "
|
||||
"[--exagg 8] [--level 1] [--spread M] [--shotDir <dir>]\n"
|
||||
" gpr_poc view-survey-all <storesDir> [--preview] [--exagg 8] "
|
||||
"[--level 1] [--shotDir <dir>]\n"
|
||||
" gpr_poc polish <storeDir> [--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";
|
||||
|
|
|
|||
Loading…
Reference in New Issue