第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();// 获取ImageDatavtkSmartPointer<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和ActorvtkSmartPointer<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, InteractorvtkSmartPointer<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 的性能可能较低。可以考虑采样部分点来显示梯度场，或者使用其他优化方法。

---

输出结果

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