diff --git a/tools/gpr_poc/main.cpp b/tools/gpr_poc/main.cpp index 1532143..9aa47ff 100644 --- a/tools/gpr_poc/main.cpp +++ b/tools/gpr_poc/main.cpp @@ -44,6 +44,9 @@ #include "render/source/ViewAdaptiveVolumeSource.hpp" #include "render/source/WholeVolumeSource.hpp" +// P9: WGS84 → CGCS2000 高斯-克吕格精确平面坐标(vendored 3DGPRViewer,零 Qt)。 +#include "CoordinateTransform.h" + // ---- VTK 离屏渲染 ---- #include #include @@ -3963,16 +3966,19 @@ int cmdViewGallery(const std::string& dir, int frames, } // ============================================================================ -// view-all:全部独立体按真实 GPS 位置/朝向摆进同一 3D 场景一起渲(测区全貌)(Task P7) +// view-all:全部独立体按精确 CGCS2000 坐标/朝向摆进同一 3D 场景一起渲(测区全貌)(P7→P9) // ============================================================================ // -// 对应客户端「选多个 ds 一起生成三维」:每条线是独立 coarse 体(线局部坐标 X=沿测线、 -// Y=通道横向、Z=深度,origin≈0),本命令把它们按各自 .gps 真实位置/航向摆进同一世界框: -// - 公共世界原点 = 全体 .gps 最小经纬(同一世界框); -// - 每条线刚体变换:平移到该线 .gps 起点局部米 + 绕竖直 Z 轴转该线航向角(起→止主方向), -// 使体局部 X(沿测线)对齐真实航向、Y(横向)随之垂直、Z(深度)保持竖直; +// 对应客户端「选多个 ds 一起生成三维」:每条线是独立密实 coarse 体(线局部坐标 X=沿测线、 +// Y=通道横向、Z=深度,origin≈0),本命令把它们按各自 .gps 真实位置/航向摆进同一世界框。 +// +// P9 升级:起点投影从 lonLatToLocalM(简化等距投影)换成 CoordinateTransform::wgs84ToCgcs2000 +// (CGCS2000 高斯-克吕格 3°带,与 P8 测绘级桥接同口径),公共带号取首条线首点经度推断, +// 全线共用同一带号 → 同一参考系;公共世界原点 = 全体 CGCS2000 最小东/北。 +// - 每条线刚体变换:平移到该线 .gps 起点 CGCS2000 局部米 + 绕竖直 Z 轴转该线航向角 +// (起→止主方向,CGCS2000 系),使体局部 X(沿测线)对齐真实航向、Y(横向)垂直、Z 竖直; // - 深度 Z 用 exagg 夸张(只 Z); -// - 每条体加载为一张整卷 vtkImageData(coarse 体小,按 --level 选层整体上纹理,不走 LOD), +// - 每条体加载为一张整卷 vtkImageData(密实体小,按 --level 选层整体上纹理,不走 LOD), // 套上该线 vtkTransform,全加进同一 renderer 一起渲。 // 传函/配色用 P4 默认醒目版(var4:增强灰度 + 实体包络 + 梯度门 + 光照),逐体按 2/98 分位标定。 // @@ -4089,8 +4095,13 @@ int cmdViewAll(int argc, char** argv) { geopro::io::gpr::GpsTrack track; }; std::vector gpsList; - double minLat = std::numeric_limits::infinity(); - double minLon = std::numeric_limits::infinity(); + // P9:CGCS2000 公共参考。带号由首条可用线首点经度推断(3°带),全线共用同一带号 → + // 同一高斯-克吕格平面参考系;公共原点取全体投影后最小东/北(数值减到较小,避免大坐标精度损失)。 + constexpr double kDeg2Rad = 3.14159265358979323846 / 180.0; + int cgcsZone = 0; + bool zoneSet = false; + double minEast = std::numeric_limits::infinity(); + double minNorth = std::numeric_limits::infinity(); for (const std::string& nm : storeNames) { const std::size_t us = nm.find_last_of('_'); if (us == std::string::npos) { @@ -4120,9 +4131,17 @@ int cmdViewAll(int argc, char** argv) { std::cerr << "[view-all] 跳过 " << nm << ":.gps 轨迹点 <2(无法定位/航向)\n"; continue; } + // 带号:首条可用线首点经度推断(与 P8 Gpr3dvSurveyVolumeBridge 同口径),全线共用。 + if (!zoneSet) { + cgcsZone = static_cast(std::lround(tr.pts.front().lon / 3.0)); + zoneSet = true; + } for (const auto& p : tr.pts) { - minLat = std::min(minLat, p.lat); - minLon = std::min(minLon, p.lon); + double cx = 0.0, cy = 0.0; // cx=北(northing), cy=东(easting,含带号) + CoordinateTransform::wgs84ToCgcs2000(p.lat * kDeg2Rad, p.lon * kDeg2Rad, + cgcsZone, cx, cy); + minEast = std::min(minEast, cy); + minNorth = std::min(minNorth, cx); } gpsList.push_back({nm, num, gpsPath, std::move(tr)}); } @@ -4130,8 +4149,10 @@ int cmdViewAll(int argc, char** argv) { std::cerr << "[view-all] 错误: 无任何线同时具备体与可用 .gps\n"; return 1; } - std::cout << "[view-all] 公共世界原点(最小经纬) lat0=" << minLat - << " lon0=" << minLon << " (共 " << gpsList.size() << " 线参与)\n"; + std::cout << "[view-all] CGCS2000 带号=" << cgcsZone + << " 公共世界原点(最小东/北)=(" << std::fixed << minEast << ", " + << minNorth << ")" << std::defaultfloat << " (共 " << gpsList.size() + << " 线参与)\n"; if (spread <= 0.0) { std::cout << "[view-all] 提示: --spread=0 用纯真实 GPS 位置;本工区为同一条路重复多趟、" "横向仅约数十米,真实位置下多趟高度重叠会叠成一条带。" @@ -4142,12 +4163,16 @@ int cmdViewAll(int argc, char** argv) { std::vector placements; int lineIdx = 0; for (const LineGps& lg : gpsList) { - // 轨迹 → 局部米(绕公共原点)。 + // 轨迹 → CGCS2000 局部米(投影到公共带号后减公共原点)。XY 约定 x=东、y=北, + // 与下游航向 atan2(hy,hx) / 刚体平移 Translate(x,y) 一致。 std::vector trackM; trackM.reserve(lg.track.pts.size()); - for (const auto& p : lg.track.pts) - trackM.push_back( - geopro::io::gpr::lonLatToLocalM(p.lat, p.lon, minLat, minLon)); + for (const auto& p : lg.track.pts) { + double cx = 0.0, cy = 0.0; // cx=北, cy=东 + CoordinateTransform::wgs84ToCgcs2000(p.lat * kDeg2Rad, p.lon * kDeg2Rad, + cgcsZone, cx, cy); + trackM.push_back(geopro::io::gpr::XY{cy - minEast, cx - minNorth}); + } const geopro::io::gpr::XY& start = trackM.front(); // 航向:起→止主方向(直接首尾向量,稳健于逐点抖动)。