From acbc6d5b4621f1d15a3ec690ad3dd51b54ff05ec Mon Sep 17 00:00:00 2001 From: gaozheng Date: Sun, 7 Jun 2026 19:17:11 +0800 Subject: [PATCH] =?UTF-8?q?spike:=20=E5=87=86=E5=A4=87=20S1=20=E5=86=92?= =?UTF-8?q?=E7=83=9F=E6=9E=84=E5=BB=BA=20+=20voxel=20=E7=A6=BB=E7=BA=BF?= =?UTF-8?q?=E9=AA=8C=E8=AF=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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 拖慢) --- CMakeLists.txt | 9 ++- tools/validate_voxel.py | 121 ++++++++++++++++++++++++++++++++++++++++ vcpkg.json | 9 +-- 3 files changed, 130 insertions(+), 9 deletions(-) create mode 100644 tools/validate_voxel.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d13477..44dcc4c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,7 +25,14 @@ endif() # ===================================================================== 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），随分层逐步启用： # find_package(GDAL CONFIG REQUIRED) diff --git a/tools/validate_voxel.py b/tools/validate_voxel.py new file mode 100644 index 0000000..d69b3c2 --- /dev/null +++ b/tools/validate_voxel.py @@ -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() diff --git a/vcpkg.json b/vcpkg.json index 06f5f6a..f250e56 100644 --- a/vcpkg.json +++ b/vcpkg.json @@ -1,15 +1,8 @@ { "name": "geopro-desktop", "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": [ - "gdal", - "proj", - "eigen3", - "spdlog", - "fmt", - "nlohmann-json", - "openssl", "gtest" ] }