spike: 准备 S1 冒烟构建 + voxel 离线验证
- tools/validate_voxel.py: 两交叉剖面 IDW+maxDist 成体素, 切片可视化; 实证 105k 体素仅 15.9% 有约束(可信体需≥3线/3D网格) - CMakeLists: VTK find_package 指定 COMPONENTS(VTK9 必需, 否则 VTK_LIBRARIES 空) - vcpkg.json: 依赖按层递增, 当前仅 gtest(免 S1 被 GDAL/PROJ 拖慢)
This commit is contained in:
parent
b219dfeae1
commit
acbc6d5b46
|
|
@ -25,7 +25,14 @@ endif()
|
||||||
# =====================================================================
|
# =====================================================================
|
||||||
|
|
||||||
find_package(Qt6 REQUIRED COMPONENTS Core Gui Widgets Network Sql Concurrent)
|
find_package(Qt6 REQUIRED COMPONENTS Core Gui Widgets Network Sql Concurrent)
|
||||||
find_package(VTK REQUIRED) # 来自 VTK_DIR(预编 install 前缀);含 GUISupportQt
|
# VTK 9 必须指定 COMPONENTS,否则 VTK_LIBRARIES 为空、链接不到任何 VTK。
|
||||||
|
# 来自 VTK_DIR(external/vtk-install)。随渲染层增补模块(Volume/Filters 等)。
|
||||||
|
find_package(VTK REQUIRED COMPONENTS
|
||||||
|
GUISupportQt
|
||||||
|
RenderingOpenGL2
|
||||||
|
InteractionStyle
|
||||||
|
FiltersSources
|
||||||
|
)
|
||||||
|
|
||||||
# 非 Qt 依赖(vcpkg),随分层逐步启用:
|
# 非 Qt 依赖(vcpkg),随分层逐步启用:
|
||||||
# find_package(GDAL CONFIG REQUIRED)
|
# find_package(GDAL CONFIG REQUIRED)
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,121 @@
|
||||||
|
#!/usr/bin/env python3
|
||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
"""
|
||||||
|
dd_voxel 离线验证 / 渲染基准:把两条交叉剖面(剖面原数据1+2)三维空间插值(IDW + maxDist 裁剪)
|
||||||
|
成体素,切片可视化。用于:
|
||||||
|
1) 提前验证 IDW + 包络裁剪算法(与 Phase 1 core/algo/IdwInterpolator 同逻辑);
|
||||||
|
2) 把"可信体"问题看实——两条 1D 测线只在其所在竖直面附近有约束,体内大部分靠插值/留空;
|
||||||
|
3) 产出 VTK 体绘制/切片的地面真值对照图。
|
||||||
|
|
||||||
|
每个测点三维坐标 = (projX, projY, z=ylist[高程]),值 = vlist。
|
||||||
|
运行: python tools/validate_voxel.py 输出: docs/_validate/voxel_*.png
|
||||||
|
依赖: numpy, scipy, matplotlib
|
||||||
|
"""
|
||||||
|
import json
|
||||||
|
import os
|
||||||
|
import sys
|
||||||
|
|
||||||
|
import numpy as np
|
||||||
|
from scipy.spatial import cKDTree
|
||||||
|
import matplotlib
|
||||||
|
matplotlib.use("Agg")
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
|
||||||
|
DATA = os.path.join(os.path.dirname(__file__), "..", "docs", "剖面网格数据的色阶数据2等文件")
|
||||||
|
OUT = os.path.join(os.path.dirname(__file__), "..", "docs", "_validate")
|
||||||
|
|
||||||
|
POWER = 2.0 # IDW 幂
|
||||||
|
K = 8 # 取最近 K 个点
|
||||||
|
MAX_DIST = 4.0 # 超过此距离(米)→ 留空(NaN),约束插值域(设计 §10)
|
||||||
|
DXY = 1.0 # 水平网格步长(米)
|
||||||
|
DZ = 0.5 # 垂向步长(米)
|
||||||
|
|
||||||
|
|
||||||
|
def load_points(name):
|
||||||
|
d = json.load(open(os.path.join(DATA, name), encoding="utf-8"))["data"]
|
||||||
|
px = np.asarray(d["projectXList"], float)
|
||||||
|
py = np.asarray(d["projectYList"], float)
|
||||||
|
z = np.asarray(d["ylist"], float) # 高程
|
||||||
|
v = np.asarray(d["vlist"], float)
|
||||||
|
return np.column_stack([px, py, z]), v
|
||||||
|
|
||||||
|
|
||||||
|
def idw_voxel(pts, vals, spec):
|
||||||
|
ox, oy, oz, nx, ny, nz, dx, dy, dz = spec
|
||||||
|
tree = cKDTree(pts)
|
||||||
|
# 网格中心点
|
||||||
|
gx = ox + (np.arange(nx) + 0.5) * dx
|
||||||
|
gy = oy + (np.arange(ny) + 0.5) * dy
|
||||||
|
gz = oz + (np.arange(nz) + 0.5) * dz
|
||||||
|
GX, GY, GZ = np.meshgrid(gx, gy, gz, indexing="ij")
|
||||||
|
Q = np.column_stack([GX.ravel(), GY.ravel(), GZ.ravel()])
|
||||||
|
dist, idx = tree.query(Q, k=K)
|
||||||
|
dist = np.maximum(dist, 1e-9)
|
||||||
|
w = 1.0 / dist ** POWER
|
||||||
|
vol = (w * vals[idx]).sum(1) / w.sum(1)
|
||||||
|
vol[dist[:, 0] > MAX_DIST] = np.nan # 最近点都超距 → 留空
|
||||||
|
return vol.reshape(nx, ny, nz), (gx, gy, gz)
|
||||||
|
|
||||||
|
|
||||||
|
def main():
|
||||||
|
os.makedirs(OUT, exist_ok=True)
|
||||||
|
p1, v1 = load_points("剖面原数据1.txt")
|
||||||
|
p2, v2 = load_points("剖面原数据2.txt")
|
||||||
|
pts = np.vstack([p1, p2])
|
||||||
|
vals = np.concatenate([v1, v2])
|
||||||
|
|
||||||
|
# 局部化(减最小,规避大坐标;同 LocalFrame 思路)
|
||||||
|
origin = pts.min(0)
|
||||||
|
ptsL = pts - origin
|
||||||
|
ext = ptsL.max(0)
|
||||||
|
nx = int(ext[0] / DXY) + 1
|
||||||
|
ny = int(ext[1] / DXY) + 1
|
||||||
|
nz = int(ext[2] / DZ) + 1
|
||||||
|
spec = (0, 0, 0, nx, ny, nz, DXY, DXY, DZ)
|
||||||
|
print(f"points={len(vals)} grid={nx}x{ny}x{nz}={nx*ny*nz} cells "
|
||||||
|
f"horiz_extent={ext[0]:.1f}x{ext[1]:.1f}m depth={ext[2]:.1f}m")
|
||||||
|
|
||||||
|
vol, (gx, gy, gz) = idw_voxel(ptsL, vals, spec)
|
||||||
|
filled = np.isfinite(vol).sum()
|
||||||
|
print(f"filled cells={filled} ({100*filled/vol.size:.1f}%) "
|
||||||
|
f"-> 其余为留空(两测线包络外无约束) | v[{np.nanmin(vol):.0f},{np.nanmax(vol):.0f}]")
|
||||||
|
|
||||||
|
vmin, vmax = np.nanpercentile(vol, [2, 98])
|
||||||
|
|
||||||
|
# 1) 中深度水平切片:应显示两条交叉的带(十字支撑)
|
||||||
|
kz = nz // 2
|
||||||
|
plt.figure(figsize=(6, 5))
|
||||||
|
plt.imshow(vol[:, :, kz].T, origin="lower", extent=[0, ext[0], 0, ext[1]],
|
||||||
|
cmap="turbo", vmin=vmin, vmax=vmax)
|
||||||
|
plt.scatter(ptsL[:, 0], ptsL[:, 1], s=1, c="k", alpha=0.15)
|
||||||
|
plt.title(f"voxel 水平切片 z={gz[kz]:.1f} (应见交叉十字支撑)")
|
||||||
|
plt.xlabel("E(m)"); plt.ylabel("N(m)"); plt.colorbar(shrink=.8)
|
||||||
|
plt.tight_layout(); plt.savefig(f"{OUT}/voxel_hslice.png", dpi=90); plt.close()
|
||||||
|
|
||||||
|
# 2) 沿 E 的垂直切片(取 profile1 附近 N):应重现剖面竖直结构
|
||||||
|
jy = ny // 2
|
||||||
|
plt.figure(figsize=(8, 3))
|
||||||
|
plt.imshow(vol[:, jy, :].T, origin="lower", extent=[0, ext[0], 0, ext[2]],
|
||||||
|
cmap="turbo", vmin=vmin, vmax=vmax, aspect="auto")
|
||||||
|
plt.title(f"voxel 垂直切片 N={gy[jy]:.1f}")
|
||||||
|
plt.xlabel("E(m)"); plt.ylabel("elev(m)"); plt.colorbar(shrink=.8)
|
||||||
|
plt.tight_layout(); plt.savefig(f"{OUT}/voxel_vslice.png", dpi=90); plt.close()
|
||||||
|
|
||||||
|
# 3) 非留空区的三维散点(直观看"十字片"支撑)
|
||||||
|
xi, yi, zi = np.where(np.isfinite(vol))
|
||||||
|
sub = slice(None, None, max(1, len(xi) // 6000))
|
||||||
|
fig = plt.figure(figsize=(7, 5)); ax = fig.add_subplot(111, projection="3d")
|
||||||
|
sc = ax.scatter(gx[xi[sub]], gy[yi[sub]], gz[zi[sub]],
|
||||||
|
c=vol[xi[sub], yi[sub], zi[sub]], cmap="turbo",
|
||||||
|
vmin=vmin, vmax=vmax, s=3)
|
||||||
|
ax.set_title("voxel 非留空区(两测线竖直面附近=十字片)")
|
||||||
|
ax.set_xlabel("E"); ax.set_ylabel("N"); ax.set_zlabel("elev")
|
||||||
|
fig.colorbar(sc, shrink=.6)
|
||||||
|
plt.tight_layout(); plt.savefig(f"{OUT}/voxel_support.png", dpi=90); plt.close()
|
||||||
|
|
||||||
|
print("written:", sorted(f for f in os.listdir(OUT) if f.startswith("voxel")))
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
sys.stdout.reconfigure(encoding="utf-8")
|
||||||
|
main()
|
||||||
|
|
@ -1,15 +1,8 @@
|
||||||
{
|
{
|
||||||
"name": "geopro-desktop",
|
"name": "geopro-desktop",
|
||||||
"version": "0.1.0",
|
"version": "0.1.0",
|
||||||
"description": "Geopro 3.0 desktop client (Qt6 + VTK9) - M1. 方案②-修订: 仅非 Qt 依赖走 vcpkg; Qt/VTK/ADS/QtKeychain 对接官方 MSVC Qt(见 docs/ENV_SETUP_Windows.md)",
|
"description": "Geopro 3.0 desktop client (Qt6 + VTK9) - M1. 方案②-修订: Qt/VTK/ADS/QtKeychain 对接官方 MSVC Qt(不走 vcpkg); 仅非 Qt 依赖走 vcpkg, 按层递增。",
|
||||||
"dependencies": [
|
"dependencies": [
|
||||||
"gdal",
|
|
||||||
"proj",
|
|
||||||
"eigen3",
|
|
||||||
"spdlog",
|
|
||||||
"fmt",
|
|
||||||
"nlohmann-json",
|
|
||||||
"openssl",
|
|
||||||
"gtest"
|
"gtest"
|
||||||
]
|
]
|
||||||
}
|
}
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue