基于python的三维射线追踪库-ttcrpy详解（5）

基于python的三维射线追踪库-ttcrpy详解（5）

ttcrpy是加拿大学者伯纳德·吉鲁（Bernard Giroux）于2021年发布的开源python库，详见（https://github.com/groupeLIAMG），参考文献（Giroux B. 2021. ttcrpy: A Python package for traveltime computation and raytracing.
SoftwareX, vol. 16, 100834. doi:10.1016/j.softx.2021.100834
）。

ttcrpy库包含了三种射线追踪方法：快速扫描算法（FSM）、最短路径法（SPM）、动节点最短路径法（DSPM）。包含其二维与三维的实现。

ttcrpy库中给出了2D矩形网格和三角形网格、3D正六面体与四面体网格等网格剖分形式，对于非规则网格，要利用python中的vtk库和pygmsh库生成，本博文研究ttcrpy中三角网格射线追踪。

1、vtk库

Vtk，（visualization toolkit）是一个开源的免费软件系统，主要用于三维计算机图形学、图像处理和可视化。Vtk是在面向对象原理的基础上设计和实现的，它的内核是用C++构建的，包含有大约250,000行代码，2000多个类，还包含有几个转换界面，因此也可以自由的通过Java，Tcl/Tk和Python各种语言使用vtk。

1.1、vtk库的安装

安装：直接在cmd中pip即可。（正常情况下会报错，可以在vtk官网下载支持python的*.whl文件安装）

pip install vtk

csdn有不少博客介绍python vtk的安装方法，本人安装的版本是9.1.0版本（目前最新版本）。

1.2、vtk库的使用

在python vtk官网（地址：https://vtk.org/download/），有一个简单示例，代码如下：

import vtk  
  
cone_a=vtk.vtkConeSource()  
  
coneMapper = vtk.vtkPolyDataMapper()  
coneMapper.SetInputConnection(cone_a.GetOutputPort())  
  
coneActor = vtk.vtkActor()  
coneActor.SetMapper(coneMapper)  
  
  
ren1= vtk.vtkRenderer()  
ren1.AddActor( coneActor )  
ren1.SetBackground( 0.1, 0.2, 0.4 )  
  
renWin = vtk.vtkRenderWindow()  
renWin.AddRenderer( ren1 )  
renWin.SetSize( 600, 600 )  
renWin.Render()  
  
iren=vtk.vtkRenderWindowInteractor()  
iren.SetRenderWindow(renWin)  
  
iren.Initialize()  
iren.Start() 

在spyder下运行结果：

vtk库包含的方法单词较长，ctrl+V比较方便，vtk说明书里面语法讲解的很详细，较易上手，但是熟练掌握要花费很多时间，ttcrpy库只需要部分vtk功能。

2、pygmsh库

Gmsh是一个强大的网格生成工具，使用脚本语言 这是出了名的难写。

pygmsh的目标是将gmsh的强大功能与python的多功能性结合起来，并 从gmsh脚本语言中提供有用的抽象，以便您可以创建复杂的 几何图形更容易。

2.1、pygmsh库安装

直接pip安装。

pip install pygmsh

我安装的版本是7.1.17版本。

2.2、pygmsh库使用

pygmsh官方文档给出了一个简单的示例（见官网https://pypi.org/project/pygmsh/），示例代码如下：

import pygmsh

with pygmsh.geo.Geometry() as geom:
    geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)
    mesh = geom.generate_mesh()

# 将网格文件保存为 vtk 格式
mesh.write("test.vtk")

将生成的vtk文件用paraview打开，显示模型：

如果每次生成vtk模型，都要用paraview打开，将非常麻烦，利用python代码可以将vtk模型显示出来。

2.2.1、二维矩形网格

利用python vtk生成二维矩形网格如下：

对应的python代码如下：

# -*- coding: utf-8 -*-
"""
Created on Thu May  5 14:48:05 2022

@author: 86159
"""

'''
此代码是用vtk绘制一个平面网格
'''

# 导入vtk库
import vtk

# 设置点与网格信息
points = vtk.vtkPoints()
cells = vtk.vtkCellArray()
polydata = vtk.vtkPolyData()
mapper = vtk.vtkPolyDataMapper()

# 图的范围
rangeX = [-10,10]
rangeY = [-10,10]
 
intervalX = 2
intervalY = 2

for gridY in range(rangeY[0], rangeY[1] + intervalY, intervalY):
    lineStart = [rangeX[0], gridY, 1.0]
    lineEnd = [rangeX[1], gridY, 1.0] 
    pointIdStart = points.InsertNextPoint(lineStart)
    pointIdEnd = points.InsertNextPoint(lineEnd)
    singleLineCell = [pointIdStart, pointIdEnd]
    cells.InsertNextCell(2,singleLineCell)
    
    print(lineStart)
    print(lineEnd)
    print(pointIdStart)
    print(pointIdEnd)
    print(singleLineCell)
    
for gridX in range(rangeX[0], rangeX[1] + intervalX, intervalX):
    lineStart = [gridX, rangeY[0], 1.0]
    lineEnd = [gridX, rangeY[1], 1.0] 
    pointIdStart = points.InsertNextPoint(lineStart)
    pointIdEnd = points.InsertNextPoint(lineEnd)
    singleLineCell = [pointIdStart, pointIdEnd]
    cells.InsertNextCell(2,singleLineCell)

 
polydata.SetLines(cells)
polydata.SetPoints(points)
mapper.SetInputData(polydata)
 
actor = vtk.vtkActor()
actor.SetMapper(mapper)
 
ren1 = vtk.vtkRenderer()
ren1.AddActor(actor)
ren1.SetBackground(0.1,0.2,0.4)
 
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren1);
renWin.SetSize(300, 300);
 
#enderWindowInteractor
iren = vtk.vtkRenderWindowInteractor();
iren.SetRenderWindow(renWin);
style = vtk.vtkInteractorStyleTrackballCamera();
iren.SetInteractorStyle(style);
 
renWin.SetSize(800, 800);
renWin.SetWindowName('test ttcrpy vtk')
renWin.Render();
iren.Initialize()
iren.Start()

2.2.2、二维三角形网格

利用python vtk生成二维三角形网格如下：

对应的python代码：

# -*- coding: utf-8 -*-
"""
Created on Fri May  6 17:03:16 2022

@author: 86159
"""

import pygmsh
import vtk
    
with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(
        [
            [0.0, 0.0],
            [1.0, 0.0],
            [1.0, 1.0],
            [0.0, 1.0],
        ],
        mesh_size=0.1,
    )
    mesh = geom.generate_mesh()
    
a = mesh.points
b = mesh.cells_dict['triangle']

points = vtk.vtkPoints()
for ai in a:
    points.InsertNextPoint(ai)
    
cells = vtk.vtkCellArray()
for bi in b:
    cells.InsertNextCell(3, bi)
    
cellColor = vtk.vtkUnsignedCharArray()    
cellColor.SetNumberOfComponents(3)        
for tmp in b:
    cellColor.InsertNextTuple(tmp)
    

cube = vtk.vtkPolyData()                  
cube.SetPoints(points)
cube.SetPolys(cells)

# # 2.创建Mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetColorModeToDefault()
mapper.SetInputData(cube)                        

# 3.创建Actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetEdgeColor(1, 0, 1)
actor.GetProperty().SetEdgeVisibility(1)    # 显示边

# 4.创建Renderer
renderer = vtk.vtkRenderer()
renderer.SetBackground(0.1, 0.2, 0.4)    # 背景白色
renderer.AddActor(actor)           # 将actor加入
renderer.ResetCamera()             # 调整显示

# 5.渲染窗口
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(renderer)
renWin.SetSize(800,800)
renWin.SetWindowName('test ttcrpy vtk')
renWin.Render()

# 6.交互
renWinInteractor = vtk.vtkRenderWindowInteractor()
renWinInteractor.SetRenderWindow(renWin)
renWinInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
renWinInteractor.Start()

2.2.3、一些其他图形

以下图形来自pygmsh官网上给出的示例，将三角形其填充颜色后显示为下列图片，具体的代码可参照pygmsh7.1.17官网，利用python vtk显示图片。

with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(