如何使用 meshlib 在三维网格表面计算两点间的测地线路径

本文详解如何利用 meshlib 的 `computegeodesicpath` 函数，在任意三角网格表面上精确计算两点间的最短测地线路径，并重点说明如何将笛卡尔坐标点正确转换为所需的 `meshtripoint` 格式。

在三维几何处理中，测地线（geodesic path）是指曲面上连接两点的局部最短路径，其长度沿表面度量而非欧氏空间直线距离。meshlib（特别是其 Python 绑定 mrmeshpy）提供了高效的测地线计算能力，但关键难点在于：输入端点必须是 MeshTriPoint 类型——即“附着于网格三角形上的点”，而非普通三维坐标。直接构造 MeshTriPoint 是不可行的；它必须通过投影（projection） 从世界坐标系中“落点”到网格表面后生成。

核心流程如下：

1. 构建或加载三角网格（mm.Mesh）；
2. 定义起点与终点的笛卡尔坐标（mm.Vector3f）；
3. 调用 mm.findProjection(point, mesh) 获取最近表面点信息，从中提取 .mtp 字段（即 MeshTriPoint）；
4. 传入 computeGeodesicPath(mesh, startMtp, stopMtp, ...) 计算路径；
5. 遍历返回的 EdgePoint 序列，用 mesh.edgePoint(ep) 提取实际三维坐标。

以下为适配您原始圆柱体示例的完整可运行代码（已修正坐标生成、投影逻辑及依赖导入）：

import numpy as np
import meshlib.mrmeshpy as mm  # 注意：应使用 mrmeshpy，非 mrmeshnumpy（后者不包含 findProjection）

def cyl2cart(rho, phi, z):
    return rho * np.cos(phi), rho * np.sin(phi), z

# 生成圆柱面网格（N×M 规则采样）
N, M = 101, 37
radius = 5.0
phi = np.linspace(0, 2*np.pi, M, endpoint=False)
z = np.linspace(-10, 10, N)

# 构造 UV 网格：x[i,j], y[i,j], z[i,j]
x = radius * np.outer(np.ones(N), np.cos(phi))
y = radius * np.outer(np.ones(N), np.sin(phi))
z_grid = np.outer(z, np.ones(M))

# 从规则网格生成三角网格（注意：meshFromUVPoints 要求 (u,v) 参数化，此处 u=z, v=phi）
mesh = mm.meshFromUVPoints(x, y, z_grid)

# 定义两个笛卡尔端点（确保在圆柱表面附近，便于投影收敛）
xp1, yp1, zp1 = cyl2cart(radius, -10/180.*np.pi, -3)
xp2, yp2, zp2 = cyl2cart(radius, 60/180.*np.pi, 8)

start_vec = mm.Vector3f(xp1, yp1, zp1)
stop_vec = mm.Vector3f(xp2, yp2, zp2)

# ✅ 关键步骤：投影到网格，获取 MeshTriPoint
start_proj = mm.findProjection(start_vec, mesh)
stop_proj = mm.findProjection(stop_vec, mesh)

if not start_proj.valid() or not stop_proj.valid():
    raise RuntimeError("Projection failed: one or both points are too far from the mesh surface.")

start_mtp = start_proj.mtp
stop_mtp = stop_proj.mtp

# 计算测地线路径（推荐使用 DijkstraBiDir 提升效率与精度）
path = mm.computeGeodesicPath(
    mesh, 
    start_mtp, 
    stop_mtp, 
    mm.GeodesicPathApprox.DijkstraBiDir
)

# 提取路径点坐标（list of (x,y,z) tuples）
geodesic_points = []
for ep in path:
    p3d = mesh.edgePoint(ep)  # EdgePoint → Vector3f
    geodesic_points.append((p3d.x, p3d.y, p3d.z))

geodesic_array = np.array(geodesic_points)

# 可视化（Plotly 示例）
import plotly.graph_objects as go

verts = mm.getNumpyVerts(mesh)
faces = mm.getNumpyFaces(mesh.topology)

fig = go.Figure(data=[
    go.Mesh3d(
        x=verts[:, 0], y=verts[:, 1], z=verts[:, 2],
        i=faces[:, 0], j=faces[:, 1], k=faces[:, 2],
        opacity=0.7, color='lightblue', name='Cylinder Mesh'
    ),
    go.Scatter3d(
        x=geodesic_array[:, 0], 
        y=geodesic_array[:, 1], 
        z=geodesic_array[:, 2],
        mode='lines+markers',
        line=dict(color='red', width=4),
        marker=dict(size=3, color='red'),
        name='Geodesic Path'
    ),
    go.Scatter3d(
        x=[xp1, xp2], y=[yp1, yp2], z=[zp1, zp2],
        mode='markers', marker=dict(size=6, color=['green', 'orange'], symbol='circle'),
        name='Start & End Points'
    )
])

fig.update_layout(
    title="Geodesic Path on Cylinder Surface",
    scene=dict(aspectmode='data'),
    margin=dict(l=0, r=0, b=0, t=40)
)
fig.show()

⚠️ 重要注意事项：

• 模块选择：务必使用 meshlib.mrmeshpy（含完整几何算法），而非 mrmeshnumpy（仅提供基础数据桥接，无 findProjection）；
• 投影鲁棒性：findProjection 要求输入点距网格足够近（默认搜索半径为网格包围盒尺寸的 10%）。若投影失败，请先调用 mm.closestPointOnMesh(point, mesh) 检查距离，或手动将点沿法向拉近表面；
• 网格质量：meshFromUVPoints 生成的是结构化四边形网格（自动三角化），适用于参数化曲面；对非结构化/扫描点云，需先用 mm.makeMeshFromPoints 或外部工具（如 Poisson 重建）生成流形三角网格；
• 性能提示：对于高分辨率网格，DijkstraBiDir 比单向 Dijkstra 更快；若需更高精度（如亚像素级），可启用 mm.GeodesicPathApprox.Exact（计算开销显著增加）。

掌握 findProjection → MeshTriPoint 这一桥梁，即可将任意三维点精准锚定至网格，解锁 meshlib 全套表面分析能力——从测地线、测地线距离场，到曲率估计与参数化映射。