第05章 10 地形梯度场模拟显示

在 VTK（Visualization Toolkit）中，可以通过计算地形数据的梯度场，并用箭头或线条来表示梯度方向和大小，从而模拟显示地形梯度场。以下是一个示例代码，展示了如何使用 VTK 和 C++ 来计算和显示地形数据的梯度场。

示例代码

#include <vtkSmartPointer.h>
#include <vtkDEMReader.h>
#include <vtkImageData.h>
#include <vtkImageGradient.h>
#include <vtkImageGradientMagnitude.h>
#include <vtkGlyph3D.h>
#include <vtkArrowSource.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkLookupTable.h>
#include <vtkColorTransferFunction.h>
#include <vtkProperty.h>
#include <vtkCamera.h>

int main(int argc, char *argv[])
{
    if (argc < 2)
    {
        std::cerr << "Usage: " << argv[0] << " <DEMFileName>" << std::endl;
        return EXIT_FAILURE;
    }

    // 读取DEM数据
    vtkSmartPointer<vtkDEMReader> reader = vtkSmartPointer<vtkDEMReader>::New();
    reader->SetFileName(argv[1]);
    reader->Update();

    // 获取ImageData
    vtkSmartPointer<vtkImageData> imageData = reader->GetOutput();

    // 计算梯度
    vtkSmartPointer<vtkImageGradient> gradientFilter = vtkSmartPointer<vtkImageGradient>::New();
    gradientFilter->SetInputData(imageData);
    gradientFilter->SetDimensionality(2); // 2D梯度计算
    gradientFilter->Update();

    // 计算梯度大小
    vtkSmartPointer<vtkImageGradientMagnitude> gradientMagnitudeFilter = vtkSmartPointer<vtkImageGradientMagnitude>::New();
    gradientMagnitudeFilter->SetInputConnection(gradientFilter->GetOutputPort());
    gradientMagnitudeFilter->Update();

    // 创建箭头源
    vtkSmartPointer<vtkArrowSource> arrowSource = vtkSmartPointer<vtkArrowSource>::New();

    // 创建Glyph3D映射
    vtkSmartPointer<vtkGlyph3D> glyph3D = vtkSmartPointer<vtkGlyph3D>::New();
    glyph3D->SetInputConnection(gradientFilter->GetOutputPort());
    glyph3D->SetSourceConnection(arrowSource->GetOutputPort());
    glyph3D->SetScaleFactor(0.1); // 设置缩放比例
    glyph3D->SetVectorModeToUseVector();
    glyph3D->Update();

    // 创建颜色查找表
    vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
    lookupTable->SetHueRange(0.667, 0.0); // 从蓝到红的渐变
    lookupTable->SetSaturationRange(1.0, 1.0);
    lookupTable->SetValueRange(1.0, 1.0);
    lookupTable->SetTableRange(0.0, 1.0);
    lookupTable->Build();

    // 创建Mapper和Actor
    vtkSmartPointer<vtkDataSetMapper> mapper = vtkSmartPointer<vtkDataSetMapper>::New();
    mapper->SetInputConnection(glyph3D->GetOutputPort());
    mapper->SetScalarRange(0.0, 1.0);
    mapper->SetLookupTable(lookupTable);

    vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);

    // 创建Renderer, RenderWindow, Interactor
    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    renderer->AddActor(actor);
    renderer->SetBackground(0.1, 0.2, 0.4); // 设置背景色

    vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
    renderWindow->AddRenderer(renderer);
    renderWindow->SetSize(800, 600);

    vtkSmartPointer<vtkRenderWindowInteractor> interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    interactor->SetRenderWindow(renderWindow);

    // 设置相机
    renderer->GetActiveCamera()->SetPosition(0, 0, 1);
    renderer->GetActiveCamera()->SetFocalPoint(0, 0, 0);
    renderer->GetActiveCamera()->SetViewUp(0, 1, 0);
    renderer->ResetCamera();

    // 开始渲染和交互
    renderWindow->Render();
    interactor->Start();

    return EXIT_SUCCESS;
}

代码说明

1. 读取DEM数据：使用 vtkDEMReader 读取 DEM 文件，获取地形数据。
2. 计算梯度：
   • 使用 vtkImageGradient 计算地形数据在 X 和 Y 方向上的梯度。
   • 使用 vtkImageGradientMagnitude 计算梯度的幅值（大小）。
3. 创建箭头源：使用 vtkArrowSource 创建箭头几何体，作为梯度方向的表示。
4. Glyph3D映射：
   • 使用 vtkGlyph3D 将每个点的梯度方向和大小映射到箭头几何体上。
   • 箭头的方向和大小会根据梯度方向和幅值进行调整。
5. 颜色查找表：使用 vtkLookupTable 为梯度场设置颜色映射，可以根据梯度大小从蓝色到红色渐变。
6. 渲染和交互：
   • 使用 vtkRenderer 将梯度场可视化。
   • 使用 vtkRenderWindow 和 vtkRenderWindowInteractor 创建交互式窗口。

编译和运行

确保你已经安装了 VTK 库，并正确配置了开发环境。编译和运行代码时，需要提供 DEM 数据文件作为命令行参数。

g++ -std=c++11 -o terrain_gradient_field terrain_gradient_field.cpp -lvtkCommonCore-9.0 -lvtkCommonDataModel-9.0 -lvtkIOImage-9.0 -lvtkImagingCore-9.0 -lvtkImagingGeneral-9.0 -lvtkRenderingCore-9.0 -lvtkRenderingOpenGL2-9.0 -lvtkInteractionStyle-9.0
./terrain_gradient_field /path/to/dem/file

注意事项

1. DEM 文件格式：确保 DEM 文件是 VTK 支持的格式，例如 GDAL 格式的 DEM 文件。
2. 箭头比例：vtkGlyph3D 的 SetScaleFactor 用于调整箭头的大小。根据数据范围和显示需求，可能需要调整比例因子。
3. 梯度计算：
   • 本示例假设地形数据是 2D 的，因此使用 vtkImageGradient 的 SetDimensionality(2) 来计算 2D 梯度。
   • 如果地形数据是 3D 的，可以调整 SetDimensionality(3)。
4. 性能优化：对于大规模数据，vtkGlyph3D 的性能可能较低。可以考虑采样部分点来显示梯度场，或者使用其他优化方法。

输出结果

运行程序后，你将看到一个交互式窗口，窗口中显示地形数据的梯度场。箭头的方向表示梯度的方向，箭头的大小表示梯度的大小。颜色从蓝色到红色表示梯度大小的变化。