feat(render): dd_voxel 体绘制(IDW->vtkImageData->GPU RayCast) + 交互切片

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

src/app/CMakeLists.txt

@@ -4,7 +4,9 @@
 find_package(VTK REQUIRED COMPONENTS
   GUISupportQt
   RenderingOpenGL2
+  RenderingVolumeOpenGL2
   InteractionStyle
+  InteractionWidgets
   FiltersGeometry
   FiltersModeling
 )

src/app/main.cpp

@@ -35,10 +35,23 @@
 #include "CameraPreset.hpp"
 #include "Scene.hpp"
 #include "actors/GridContourActor.hpp"
+#include "actors/VoxelActor.hpp"
+
+#include "algo/IInterpolator.hpp"
+#include "algo/IdwInterpolator.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
 
 #include <QVTKOpenGLStereoWidget.h>
 #include <vtkGenericOpenGLRenderWindow.h>
+#include <vtkImageData.h>
+#include <vtkImagePlaneWidget.h>
+#include <vtkRenderWindowInteractor.h>
 #include <vtkRenderer.h>
+#include <vtkSmartPointer.h>
+#include <vtkVolume.h>
 
 namespace {
 
@@ -71,6 +84,93 @@ std::string readPem(const std::string& path)
     return ss.str();
 }
 
+// dd_voxel 构建结果：体绘制 volume + 其底层 vtkImageData（供切片 widget）+ 网格信息。
+struct VoxelBuildResult {
+    vtkSmartPointer<vtkVolume> volume;
+    vtkSmartPointer<vtkImageData> image;
+    int nx = 0, ny = 0, nz = 0;
+    double vmin = 0.0, vmax = 0.0;
+};
+
+// dd_voxel：合并两条交叉剖面散点 → 局部化 → IDW 体素 → GPU 体绘制。
+// 与 tools/validate_voxel.py 同逻辑（power=2, maxDist≈4, dxy=1m, dz=0.5m）。
+VoxelBuildResult buildVoxelFromScatters(
+    const std::vector<geopro::core::ScatterField>& scatters,
+    const geopro::core::ColorScale& cs)
+{
+    using geopro::core::GridSpec;
+    using geopro::core::IdwInterpolator;
+    using geopro::core::PointSet;
+
+    constexpr double kDxy = 1.0;       // 水平步长(米)
+    constexpr double kDz = 0.5;        // 垂向步长(米)
+    constexpr double kPower = 2.0;     // IDW 幂
+    constexpr double kMaxDist = 4.0;   // 超距留空(NaN)
+
+    // 合并两剖面点为一组 (X=projX, Y=projY, Z=z)。
+    std::vector<double> X, Y, Z, V;
+    for (const auto& s : scatters) {
+        const size_t n = s.v.size();
+        for (size_t p = 0; p < n; ++p) {
+            // projX/projY 为 GIS 平面坐标；z 为高程；v 为值。容错缺失字段。
+            if (p < s.projX.size() && p < s.projY.size() && p < s.z.size()) {
+                X.push_back(s.projX[p]);
+                Y.push_back(s.projY[p]);
+                Z.push_back(s.z[p]);
+                V.push_back(s.v[p]);
+            }
+        }
+    }
+    if (V.empty()) return VoxelBuildResult{};
+
+    // 局部化：三轴各减最小值（规避大坐标；ox=oy=oz=0）。
+    const double minX = *std::min_element(X.begin(), X.end());
+    const double minY = *std::min_element(Y.begin(), Y.end());
+    const double minZ = *std::min_element(Z.begin(), Z.end());
+    PointSet pts;
+    double extX = 0.0, extY = 0.0, extZ = 0.0;
+    for (size_t p = 0; p < V.size(); ++p) {
+        const double lx = X[p] - minX, ly = Y[p] - minY, lz = Z[p] - minZ;
+        pts.x.push_back(lx);
+        pts.y.push_back(ly);
+        pts.z.push_back(lz);
+        pts.v.push_back(V[p]);
+        extX = std::max(extX, lx);
+        extY = std::max(extY, ly);
+        extZ = std::max(extZ, lz);
+    }
+
+    GridSpec spec;
+    spec.ox = spec.oy = spec.oz = 0.0;
+    spec.dx = spec.dy = kDxy;
+    spec.dz = kDz;
+    spec.nx = static_cast<int>(extX / kDxy) + 1;
+    spec.ny = static_cast<int>(extY / kDxy) + 1;
+    spec.nz = static_cast<int>(extZ / kDz) + 1;
+    spec.power = kPower;
+    spec.maxDist = kMaxDist;
+
+    const auto vol = IdwInterpolator().interpolate(pts, spec);
+
+    // 算 vmin/vmax（忽略 NaN）。
+    double vmin = std::numeric_limits<double>::max();
+    double vmax = std::numeric_limits<double>::lowest();
+    for (double v : vol.data()) {
+        if (std::isnan(v)) continue;
+        vmin = std::min(vmin, v);
+        vmax = std::max(vmax, v);
+    }
+    if (vmin > vmax) { vmin = 0.0; vmax = 1.0; }  // 全 NaN 兜底
+
+    VoxelBuildResult r;
+    r.image = vtkSmartPointer<vtkImageData>::New();
+    r.volume = geopro::render::buildVoxel(
+        vol, cs, spec.ox, spec.oy, spec.oz, spec.dx, spec.dy, spec.dz, vmin, vmax, r.image);
+    r.nx = spec.nx; r.ny = spec.ny; r.nz = spec.nz;
+    r.vmin = vmin; r.vmax = vmax;
+    return r;
+}
+
 // 在给定 QMainWindow 上构建 M1 工作台：ADS 三栏 + 对象树 → 渲染联动 + 属性面板。
 // repo 生命周期须覆盖到事件循环结束（由调用方保证）。
 // 相机模式：默认二维俯视。
@@ -118,18 +218,30 @@ void buildWorkbench(QMainWindow& window, geopro::data::LocalSampleRepository& re
     cameraGroup->addAction(act2D);
     cameraGroup->addAction(act3D);
     act2D->setChecked(true);  // 默认二维
+    viewToolBar->addSeparator();
+    auto* actVoxel = viewToolBar->addAction(QStringLiteral("dd_voxel"));
     centerLayout->addWidget(viewToolBar);
     centerLayout->addWidget(vtkWidget, 1);
 
+    // 交互切片 widget：dd_voxel 时对体素 vtkImageData 加可拖切面（dd_slice）。
+    // 需 interactor（render window 已设），切回非体素视图时 Off。用 shared_ptr 跨 lambda 共享。
+    auto sliceWidget = std::make_shared<vtkSmartPointer<vtkImagePlaneWidget>>();
+
+    auto hideSlice = [sliceWidget]() {
+        if (sliceWidget->Get()) (*sliceWidget)->Off();
+    };
+
     QObject::connect(act2D, &QAction::triggered, vtkWidget,
-                     [cameraMode, rendererPtr, renderWindowPtr]() {
+                     [cameraMode, rendererPtr, renderWindowPtr, hideSlice]() {
                          *cameraMode = CameraMode::Top2D;
+                         hideSlice();  // 切回剖面视图：关闭体素切片
                          geopro::render::applyTop2D(rendererPtr);
                          renderWindowPtr->Render();
                      });
     QObject::connect(act3D, &QAction::triggered, vtkWidget,
-                     [cameraMode, rendererPtr, renderWindowPtr]() {
+                     [cameraMode, rendererPtr, renderWindowPtr, hideSlice]() {
                          *cameraMode = CameraMode::Free3D;
+                         hideSlice();
                          geopro::render::applyFree3D(rendererPtr);
                          renderWindowPtr->Render();
                      });
@@ -182,6 +294,58 @@ void buildWorkbench(QMainWindow& window, geopro::data::LocalSampleRepository& re
     QObject::connect(tree, &QTreeWidget::itemClicked, tree,
                      [renderDataset](QTreeWidgetItem* it, int) { renderDataset(it); });
 
+    // 当前体素 vtkImageData（dd_slice 切面 widget 需复用，保活到下次重建）。
+    auto voxelImage = std::make_shared<vtkSmartPointer<vtkImageData>>();
+
+    // dd_voxel：读两条交叉剖面 → IDW 体素 → GPU 体绘制 + 默认 3D 视角 + 可拖切片。
+    QObject::connect(actVoxel, &QAction::triggered, vtkWidget,
+                     [&repo, scene, rendererPtr, renderWindowPtr, cameraMode, act3D, propLabel,
+                      sliceWidget, voxelImage]() {
+                         const auto scatters = repo.loadVoxelScatters();
+                         const auto cs = repo.loadColorScale("grid1");
+                         auto result = buildVoxelFromScatters(scatters, cs);
+                         if (!result.volume) {
+                             propLabel->setText(QStringLiteral("dd_voxel: 无可用散点数据"));
+                             return;
+                         }
+
+                         // 关闭旧切片（避免引用已释放的旧 image）。
+                         if (sliceWidget->Get()) (*sliceWidget)->Off();
+
+                         scene->clear();
+                         rendererPtr->AddViewProp(result.volume);
+                         *voxelImage = result.image;
+
+                         // 体素默认 3D 视角（同步工具条勾选态）。
+                         *cameraMode = CameraMode::Free3D;
+                         act3D->setChecked(true);
+                         geopro::render::applyFree3D(rendererPtr);
+                         rendererPtr->ResetCamera();
+                         renderWindowPtr->Render();  // 先渲一帧确保 interactor 就绪
+
+                         // 交互切片 widget：对体素 vtkImageData 加 X 轴可拖切面（dd_slice）。
+                         vtkRenderWindowInteractor* interactor = renderWindowPtr->GetInteractor();
+                         if (interactor) {
+                             if (!sliceWidget->Get())
+                                 *sliceWidget = vtkSmartPointer<vtkImagePlaneWidget>::New();
+                             (*sliceWidget)->SetInteractor(interactor);
+                             (*sliceWidget)->SetInputData(result.image);
+                             (*sliceWidget)->SetPlaneOrientationToXAxes();
+                             (*sliceWidget)->SetSliceIndex(result.nx / 2);
+                             (*sliceWidget)->On();
+                         }
+
+                         renderWindowPtr->Render();
+                         propLabel->setText(
+                             QStringLiteral("dd_voxel 体素\n网格: %1 x %2 x %3\nvmin / vmax: %4 / %5\n"
+                                            "(拖动切面查看 dd_slice)")
+                                 .arg(result.nx)
+                                 .arg(result.ny)
+                                 .arg(result.nz)
+                                 .arg(result.vmin)
+                                 .arg(result.vmax));
+                     });
+
     // 默认渲染第一个 DS，让窗口一打开就有图。
     if (auto* first = tree->topLevelItemCount() > 0 ? tree->topLevelItem(0) : nullptr) {
         // 递归找到首个带 dsId 的项。

src/data/repo/LocalSampleRepository.cpp

@@ -20,6 +20,9 @@ constexpr const char* kDsId = "grid1";
 constexpr const char* kGridFile = u8"剖面网格数据1.txt";
 constexpr const char* kColorScaleFile = u8"剖面网格数据的色阶数据1.txt";
 constexpr const char* kScatterFile = u8"剖面原数据1.txt";
+// dd_voxel：两条交叉剖面散点（约 78° 相交），合并为体素插值输入。
+constexpr const char* kVoxelScatterFile1 = u8"剖面原数据1.txt";
+constexpr const char* kVoxelScatterFile2 = u8"剖面原数据2.txt";
 constexpr const char* kAnomalyFile = u8"剖面网格数据1——对应的异常圈定数据.txt";
 
 // 校验 dsId，未知则抛错（输入边界验证）。
@@ -89,4 +92,11 @@ std::vector<Anomaly> LocalSampleRepository::loadAnomalies(const std::string& dsI
     return parseAnomalies(readFile(kAnomalyFile));
 }
 
+std::vector<ScatterField> LocalSampleRepository::loadVoxelScatters() {
+    std::vector<ScatterField> out;
+    out.push_back(parseScatter(readFile(kVoxelScatterFile1)));
+    out.push_back(parseScatter(readFile(kVoxelScatterFile2)));
+    return out;
+}
+
 }  // namespace geopro::data

src/data/repo/LocalSampleRepository.hpp

@@ -17,6 +17,10 @@ public:
     geopro::core::ColorScale loadColorScale(const std::string& dsId) override;
     std::vector<geopro::core::Anomaly> loadAnomalies(const std::string& dsId) override;
 
+    // dd_voxel 输入：读两条交叉剖面散点（剖面原数据1.txt + 剖面原数据2.txt），
+    // 返回两者，供体素插值合并。具体类专有方法（不进 IDatasetRepository 接口）。
+    std::vector<geopro::core::ScatterField> loadVoxelScatters();
+
 private:
     // 用 QFile（路径 = fromUtf8(dir)+fromUtf8(name)）读全文件并转 UTF-8 std::string。
     // 读失败抛 std::runtime_error。

src/render/CMakeLists.txt

@@ -1,6 +1,6 @@
-find_package(VTK REQUIRED COMPONENTS CommonCore CommonDataModel FiltersGeometry FiltersModeling RenderingOpenGL2 InteractionStyle)
+find_package(VTK REQUIRED COMPONENTS CommonCore CommonDataModel FiltersGeometry FiltersModeling RenderingOpenGL2 RenderingVolumeOpenGL2 InteractionStyle InteractionWidgets)
 add_library(geopro_render STATIC
-  Scene.cpp ColorLutBuilder.cpp CameraPreset.cpp actors/GridContourActor.cpp)
+  Scene.cpp ColorLutBuilder.cpp CameraPreset.cpp actors/GridContourActor.cpp actors/VoxelActor.cpp)
 target_include_directories(geopro_render PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 target_link_libraries(geopro_render PUBLIC geopro_core ${VTK_LIBRARIES})
 target_compile_features(geopro_render PUBLIC cxx_std_17)

src/render/actors/VoxelActor.cpp

@@ -0,0 +1,99 @@
+#include "actors/VoxelActor.hpp"
+
+#include <cmath>
+#include <limits>
+
+#include <vtkColorTransferFunction.h>
+#include <vtkDoubleArray.h>
+#include <vtkGPUVolumeRayCastMapper.h>
+#include <vtkNew.h>
+#include <vtkPiecewiseFunction.h>
+#include <vtkPointData.h>
+#include <vtkVolumeProperty.h>
+
+namespace geopro::render {
+
+namespace {
+
+// 颜色/不透明度传递函数采样级数。
+constexpr int kTransferSamples = 64;
+// 体绘制最大不透明度（[vmin,vmax] 线性 0→kMaxOpacity）。
+constexpr double kMaxOpacity = 0.15;
+
+// NaN/留空格的哨兵：落在 [vmin,vmax] 之外，传递函数把它映射为完全透明。
+double sentinel(double vmin) { return vmin - 1.0; }
+
+}  // namespace
+
+vtkSmartPointer<vtkVolume> buildVoxel(const geopro::core::ScalarVolume& vol,
+                                      const geopro::core::ColorScale& cs,
+                                      double ox, double oy, double oz,
+                                      double dx, double dy, double dz,
+                                      double vmin, double vmax,
+                                      vtkSmartPointer<vtkImageData>& outImage)
+{
+    const int nx = vol.nx(), ny = vol.ny(), nz = vol.nz();
+
+    // vmin/vmax 退化兜底，避免传递函数区间为零。
+    if (vmin >= vmax) vmax = vmin + 1.0;
+    const double blank = sentinel(vmin);
+
+    auto img = vtkSmartPointer<vtkImageData>::New();
+    img->SetDimensions(nx, ny, nz);
+    img->SetOrigin(ox, oy, oz);
+    img->SetSpacing(dx, dy, dz);
+
+    vtkNew<vtkDoubleArray> sc;
+    sc->SetName("v");
+    sc->SetNumberOfTuples(static_cast<vtkIdType>(nx) * ny * nz);
+    // 点序 i 最快、j 次之、k 最慢（匹配 vtkImageData 与 ScalarVolume::idx）。
+    for (int k = 0; k < nz; ++k)
+        for (int j = 0; j < ny; ++j)
+            for (int i = 0; i < nx; ++i) {
+                const double v = vol.at(i, j, k);
+                const vtkIdType id = (static_cast<vtkIdType>(k) * ny + j) * nx + i;
+                sc->SetValue(id, std::isnan(v) ? blank : v);  // NaN → 哨兵
+            }
+    img->GetPointData()->SetScalars(sc);
+    outImage = img;
+
+    // 颜色传递函数：在 [vmin,vmax] 按 ColorScale 采样若干 RGB 点。
+    vtkNew<vtkColorTransferFunction> color;
+    for (int t = 0; t < kTransferSamples; ++t) {
+        const double val = vmin + (vmax - vmin) * t / (kTransferSamples - 1);
+        const auto c = cs.colorAt(val);
+        color->AddRGBPoint(val, c.r / 255.0, c.g / 255.0, c.b / 255.0);
+    }
+
+    // 不透明度传递函数：哨兵 → 0（透明）；[vmin,vmax] 线性递增到 kMaxOpacity。
+    vtkNew<vtkPiecewiseFunction> opacity;
+    opacity->AddPoint(blank, 0.0);
+    opacity->AddPoint(vmin, 0.0);
+    opacity->AddPoint(vmax, kMaxOpacity);
+
+    vtkNew<vtkGPUVolumeRayCastMapper> mapper;
+    mapper->SetInputData(img);
+
+    vtkNew<vtkVolumeProperty> prop;
+    prop->SetColor(color);
+    prop->SetScalarOpacity(opacity);
+    prop->SetInterpolationTypeToLinear();
+    prop->ShadeOff();
+
+    auto volume = vtkSmartPointer<vtkVolume>::New();
+    volume->SetMapper(mapper);
+    volume->SetProperty(prop);
+    return volume;
+}
+
+vtkSmartPointer<vtkVolume> buildVoxel(const geopro::core::ScalarVolume& vol,
+                                      const geopro::core::ColorScale& cs,
+                                      double ox, double oy, double oz,
+                                      double dx, double dy, double dz,
+                                      double vmin, double vmax)
+{
+    vtkSmartPointer<vtkImageData> ignored;
+    return buildVoxel(vol, cs, ox, oy, oz, dx, dy, dz, vmin, vmax, ignored);
+}
+
+}  // namespace geopro::render

src/render/actors/VoxelActor.hpp

@@ -0,0 +1,26 @@
+#pragma once
+#include <vtkSmartPointer.h>
+#include <vtkVolume.h>
+#include <vtkImageData.h>
+#include "model/Field.hpp"
+#include "model/ColorScale.hpp"
+namespace geopro::render {
+
+// 把 core 规则标量体（IDW 输出，含 NaN 留空）转 vtkImageData，再建 GPU 光线投射体绘制。
+// 颜色按 ColorScale 在 [vmin,vmax] 采样；NaN/留空格 → 不透明度 0（透明）。
+// 返回 vtkVolume（由调用方加入 renderer）。
+vtkSmartPointer<vtkVolume> buildVoxel(const geopro::core::ScalarVolume& vol,
+                                      const geopro::core::ColorScale& cs,
+                                      double ox, double oy, double oz,
+                                      double dx, double dy, double dz,
+                                      double vmin, double vmax);
+
+// 同上，但额外通过 out 参数暴露内部 vtkImageData，供 main 建交互切片 widget。
+vtkSmartPointer<vtkVolume> buildVoxel(const geopro::core::ScalarVolume& vol,
+                                      const geopro::core::ColorScale& cs,
+                                      double ox, double oy, double oz,
+                                      double dx, double dy, double dz,
+                                      double vmin, double vmax,
+                                      vtkSmartPointer<vtkImageData>& outImage);
+
+}  // namespace geopro::render

tests/CMakeLists.txt

@@ -57,8 +57,10 @@ endif()
 
 # render 层：ColorLutBuilder（core ColorScale -> vtkLookupTable）。
 # 需 vtkLookupTable（VTK::CommonCore）；geopro_render 已 PUBLIC 传递其余 VTK 组件。
-find_package(VTK REQUIRED COMPONENTS CommonCore)
+find_package(VTK REQUIRED COMPONENTS CommonCore CommonDataModel)
 target_sources(geopro_tests PRIVATE render/test_color_lut.cpp)
+# dd_voxel：buildVoxel(ScalarVolume->vtkImageData->GPU 体绘制) 构建不崩 + dims 正确。
+target_sources(geopro_tests PRIVATE render/test_voxel_build.cpp)
 target_link_libraries(geopro_tests PRIVATE geopro_render ${VTK_LIBRARIES})
 vtk_module_autoinit(TARGETS geopro_tests MODULES ${VTK_LIBRARIES})
 

tests/render/test_voxel_build.cpp

@@ -0,0 +1,56 @@
+#include <gtest/gtest.h>
+
+#include <cmath>
+#include <limits>
+
+#include <vtkImageData.h>
+
+#include "actors/VoxelActor.hpp"
+#include "model/ColorScale.hpp"
+#include "model/Field.hpp"
+
+using namespace geopro::core;
+
+// buildVoxel 用小 ScalarVolume 构建：不崩、返回非空 volume、image dims/spacing/origin 正确。
+TEST(VoxelBuild, BuildsVolumeFromSmallScalarVolume) {
+    // 2x3x4 小体，填一些值 + 一个 NaN（验证哨兵透明分支）。
+    ScalarVolume vol(2, 3, 4);
+    double n = 0.0;
+    for (int k = 0; k < 4; ++k)
+        for (int j = 0; j < 3; ++j)
+            for (int i = 0; i < 2; ++i) vol.at(i, j, k) = n++;
+    vol.at(1, 1, 1) = std::numeric_limits<double>::quiet_NaN();
+
+    ColorScale cs;
+    cs.addStop(0.0, Rgba{0, 0, 255, 255});
+    cs.addStop(24.0, Rgba{255, 0, 0, 255});
+
+    vtkSmartPointer<vtkImageData> image;
+    auto volume = geopro::render::buildVoxel(
+        vol, cs, /*ox*/ 1.0, /*oy*/ 2.0, /*oz*/ 3.0,
+        /*dx*/ 1.0, /*dy*/ 1.0, /*dz*/ 0.5, /*vmin*/ 0.0, /*vmax*/ 24.0, image);
+
+    ASSERT_NE(volume.GetPointer(), nullptr);
+    ASSERT_NE(image.GetPointer(), nullptr);
+
+    int dims[3];
+    image->GetDimensions(dims);
+    EXPECT_EQ(dims[0], 2);
+    EXPECT_EQ(dims[1], 3);
+    EXPECT_EQ(dims[2], 4);
+
+    double sp[3];
+    image->GetSpacing(sp);
+    EXPECT_DOUBLE_EQ(sp[0], 1.0);
+    EXPECT_DOUBLE_EQ(sp[1], 1.0);
+    EXPECT_DOUBLE_EQ(sp[2], 0.5);
+
+    double org[3];
+    image->GetOrigin(org);
+    EXPECT_DOUBLE_EQ(org[0], 1.0);
+    EXPECT_DOUBLE_EQ(org[1], 2.0);
+    EXPECT_DOUBLE_EQ(org[2], 3.0);
+
+    // NaN 格应被替换为哨兵（vmin-1 = -1），不应残留 NaN。
+    EXPECT_FALSE(std::isnan(image->GetScalarComponentAsDouble(1, 1, 1, 0)));
+}