基于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库中包含了矩形网格和三角形网格两种网格剖分形式，对于三角网格，要利用python中的vtk库和pygmsh库，本文研究ttcrpy中三角网格射线追踪。
文章目录
1、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库
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(
[
[0.0, 0.0],
[1.0, -0.2],
[1.1, 1.2],
[0.1, 0.7],
],
mesh_size=0.1,
)
mesh = geom.generate_mesh()
# ellpsoid with holes
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_max = 0.1
ellipsoid = geom.add_ellipsoid([0.0, 0.0, 0.0], [1.0, 0.7, 0.5])
cylinders = [
geom.add_cylinder([-1.0, 0.0, 0.0], [2.0, 0.0, 0.0], 0.3),
geom.add_cylinder([0.0, -1.0, 0.0], [0.0, 2.0, 0.0], 0.3),
geom.add_cylinder([0.0, 0.0, -1.0], [0.0, 0.0, 2.0], 0.3),
]
geom.boolean_difference(ellipsoid, geom.boolean_union(cylinders))
mesh = geom.generate_mesh()
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_min = 0.1
geom.characteristic_length_max = 0.1
rectangle = geom.add_rectangle([-1.0, -1.0, 0.0], 2.0, 2.0)
disk1 = geom.add_disk([-1.2, 0.0, 0.0], 0.5)
disk2 = geom.add_disk([+1.2, 0.0, 0.0], 0.5)
disk3 = geom.add_disk([0.0, -0.9, 0.0], 0.5)
disk4 = geom.add_disk([0.0, +0.9, 0.0], 0.5)
flat = geom.boolean_difference(
geom.boolean_union([rectangle, disk1, disk2]),
geom.boolean_union([disk3, disk4]),
)
geom.extrude(flat, [0, 0, 0.3])
mesh = geom.generate_mesh()
with pygmsh.geo.Geometry() as geom:
poly = geom.add_polygon(
[
[0.0, 0.0],
[2.0, 0.0],
[3.0, 1.0],
[1.0, 2.0],
[0.0, 1.0],
],
mesh_size=0.3,
)
field0 = geom.add_boundary_layer(
edges_list=[poly.curves[0]],
lcmin=0.05,
lcmax=0.2,
distmin=0.0,
distmax=0.2,
)
field1 = geom.add_boundary_layer(
nodes_list=[poly.points[2]],
lcmin=0.05,
lcmax=0.2,
distmin=0.1,
distmax=0.4,
)
geom.set_background_mesh([field0, field1], operator="Min")
mesh = geom.generate_mesh()
from math import pi
import pygmsh
with pygmsh.geo.Geometry() as geom:
poly = geom.add_polygon(
[
[+0.0, +0.5],
[-0.1, +0.1],
[-0.5, +0.0],
[-0.1, -0.1],
[+0.0, -0.5],
[+0.1, -0.1],
[+0.5, +0.0],
[+0.1, +0.1],
],
mesh_size=0.05,
)
geom.twist(
poly,
translation_axis=[0, 0, 1],
rotation_axis=[0, 0, 1],
point_on_axis=[0, 0, 0],
angle=pi / 3,
)
mesh = geom.generate_mesh()
3、ttrpy中tmesh三角网射线追踪
3.1、导入对应的库
需要导入的库有time、numpy、matplotlib.pyplot、pygmsh、ttcrpy中的tmesh、numpy、vtk等。
import time
import numpy as np
import matplotlib.pyplot as plt
import pygmsh
from ttcrpy import tmesh
from numpy.random import default_rng
import vtk
3.2、定义vtk中显示网格的函数
此部分代码如下，基本上复制上述代码即可。
def draw_vtk(mesh):
'''
在vtk中显示网格模型
'''
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(1, 0, 0) # 背景白色
renderer.AddActor(actor) # 将actor加入
renderer.ResetCamera() # 调整显示
# 5.渲染窗口
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(renderer)
renWin.SetSize(600,600)
renWin.SetWindowName('meshio example')
renWin.Render()
# 6.交互
renWinInteractor = vtk.vtkRenderWindowInteractor()
renWinInteractor.SetRenderWindow(renWin)
renWinInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
renWinInteractor.Start()
3.3、设置网格
设置一个大小为1×1×1的三维立体网格，三角形单元的边长为0.042，显示图形
with pygmsh.geo.Geometry() as geom:
geom.add_box(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
mesh_size = 0.042)
msh = geom.generate_mesh()
# 显示模型
draw_vtk(msh)
模型显示结果如下所示：
3.4、设置慢度模型
根据ttcrpy官方文档中示例，设置模型参数如下：
def f(z):
return 1.0 / (1.5 + 4.5*z)
slowness = f(msh.points[:,2])
rng = default_rng(1966)
src = rng.uniform(0.01, 0.99, (10, 3)) # 10 src points
rcv = rng.uniform(0.05, 0.95, (20, 3)) # 20 receivers
# src & rcv must have the same number of rows, with rows corresponding to src-rcv pairs
src = np.kron(np.ones((20, 1)), src)
rcv = np.kron(rcv, np.ones((10, 1)))
3.5、射线追踪并计算走时
运行以下代码，输出计算得到的走时。
tm = tmesh.Mesh3d(msh.points,msh.cells_dict['tetra'],
n_threads = 1, cell_slowness=0, method='DSPM',
gradient_method = 1, tt_from_rp = 1,
process_vel = False,
maxit = 20, min_dist = 1e-5,
n_secondary = 2, n_tertiary = 2,
radius_factor_tertiary=3, translate_grid = False)
tm.set_slowness(slowness)
tstart = time.time()
tt1 = tm.raytrace(src, rcv)
tend1 = time.time() - tstart
print('time serial = ', tend1)
计算结果如下
计算耗时：time serial = 13.381624937057495
走时曲线
以上就是完整过程了，详见官方文档。
4、参考资料
Giroux B. 2021. ttcrpy: A Python package for traveltime computation and raytracing.
SoftwareX, vol. 16, 100834. doi: 10.1016/j.softx.2021.100834
https://github.com/groupeLIAMG/ttcr
