第05章 11 动量剖面可视化代码一则

在计算流体力学（CFD）中，动量剖面（Momentum Profiles）通常用于描述流体在流动方向上的动量分布。在 VTK 中，可以通过读取速度场数据，并计算和展示动量剖面来可视化呈现速度场信息。

示例代码

以下是一个示例代码，展示如何使用 VTK 和 C++ 读取速度场数据，并计算和可视化动量剖面。

#include <vtkSmartPointer.h>
#include <vtkXMLImageReader.h>
#include <vtkImageData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkMath.h>
#include <vtkPointLocator.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkLine.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkNamedColors.h>
#include <vtkProperty.h>
#include <vtkCamera.h>int main(int argc, char *argv[])
{if (argc < 2){std::cerr << "Usage: " << argv[0] << " <InputVelocityFile.vti>" << std::endl;return EXIT_FAILURE;}// 读取速度场数据vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();reader->SetFileName(argv[1]);reader->Update();vtkSmartPointer<vtkImageData> imageData = reader->GetOutput();// 获取速度数组vtkSmartPointer<vtkDoubleArray> velocityArray = vtkDoubleArray::SafeDownCast(imageData->GetPointData()->GetVectors());if (!velocityArray){std::cerr << "No velocity vectors found in the input data." << std::endl;return EXIT_FAILURE;}// 定义剖面线的起点和终点double startPoint[3] = {0.0, 0.0, 0.0};double endPoint[3] = {1.0, 0.0, 0.0};// 创建剖面线的点数组vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();// 创建点定位器vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New();locator->SetDataSet(imageData);locator->BuildLocator();// 采样点的数量int numberOfSamples = 100;// 计算剖面线上的点for (int i = 0; i <= numberOfSamples; i++){double t = static_cast<double>(i) / numberOfSamples;double point[3];for (int j = 0; j < 3; j++){point[j] = startPoint[j] * (1.0 - t) + endPoint[j] * t;}points->InsertNextPoint(point);}// 查找最近的图像数据点vtkIdType ptId;double closestPoint[3];vtkSmartPointer<vtkDoubleArray> momentumArray = vtkSmartPointer<vtkDoubleArray>::New();momentumArray->SetNumberOfComponents(1);momentumArray->SetName("Momentum");for (vtkIdType i = 0; i < points->GetNumberOfPoints(); i++){points->GetPoint(i, closestPoint);ptId = locator->FindClosestPoint(closestPoint);double velocity[3];velocityArray->GetTuple(ptId, velocity);double momentum = vtkMath::Dot(velocity, velocity); // 假设密度为1，动量等于速度的平方momentumArray->InsertNextValue(momentum);}// 创建剖面线的PolyDatavtkSmartPointer<vtkPolyData> profileLine = vtkSmartPointer<vtkPolyData>::New();profileLine->SetPoints(points);// 创建线单元vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();for (int i = 0; i < numberOfSamples; i++){line->GetPointIds()->SetId(0, i);line->GetPointIds()->SetId(1, i + 1);lines->InsertNextCell(line);}profileLine->SetLines(lines);profileLine->GetPointData()->AddArray(momentumArray);// 创建LookupTablevtkSmartPointer<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(momentumArray->GetRange());lookupTable->Build();// 创建Mapper和ActorvtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();mapper->SetInputData(profileLine);mapper->SetScalarModeToUsePointData();mapper->SetColorModeToMapScalars();mapper->SelectColorArray("Momentum");mapper->SetScalarRange(momentumArray->GetRange());vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();actor->SetMapper(mapper);actor->GetProperty()->SetLineWidth(3.0);// 创建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. 读取速度场数据：使用 vtkXMLImageDataReader 读取速度场的 VTI 文件。
2. 获取速度数组：从图像数据的点数据中获取速度向量数组。
3. 定义剖面线：定义剖面线的起点和终点，以及采样点的数量。
4. 采样剖面线上的点：在剖面线上均匀采样点，并使用 vtkPointLocator 查找图像数据中最近的点。
5. 计算动量：对于每个采样点，获取对应的速度向量，并计算动量（这里假设密度为1，动量等于速度的平方）。
6. 创建剖面线的PolyData：将采样点和动量数据组织成PolyData，包含点、线单元和动量数组。
7. 创建LookupTable：设置颜色查找表，根据动量值从蓝色到红色渐变。
8. 创建Mapper和Actor：使用PolyData创建Mapper，并设置颜色映射。创建Actor并设置线宽。
9. 渲染和交互：设置Renderer、RenderWindow和Interactor，配置相机，并开始渲染和交互。

---

编译和运行

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

g++ -std=c++11 -o momentum_profiles momentum_profiles.cpp -lvtkCommonCore-9.0 -lvtkCommonDataModel-9.0 -lvtkIOXML-9.0 -lvtkImagingCore-9.0 -lvtkRenderingContext2D-9.0 -lvtkRenderingCore-9.0 -lvtkRenderingFreeType-9.0 -lvtkRenderingOpenGL2-9.0 -lvtkInteractionStyle-9.0
./momentum_profiles /path/to/velocity/file.vti

---

注意事项

1. 速度场数据：确保输入的 VTI 文件包含速度向量数据。
2. 剖面线定义：可以根据需要定义不同的剖面线，例如在不同方向上采样。
3. 动量计算：本示例中假设密度为1，动量等于速度的平方。实际应用中，可能需要根据具体物理模型调整动量的计算公式。
4. 性能优化：对于大规模数据，剖面线的采样点数可能需要调整，以平衡可视化质量和计算效率。

---

输出结果

运行程序后，你将看到一个交互式窗口，窗口中显示剖面线上的动量分布。剖面线的颜色从蓝色到红色表示动量值的变化。可以通过交互式窗口进行缩放和旋转，以更好地观察剖面线的动量分布情况。