pythonpath3dsurface

Geodesic path between two points on surface using meshlib


I would like to use Python to get the path coordinates of the shortest (geodesic) path on a surface in 3D between two points. The solution is somehow already shown in the accepted answer in this post, however, I am not able to reproduce the result of Fedor, since I don't understand the format that is expected by meshlib to provide the point locations.

The geodesic path calculation is done with the computeGeodesicPath function which expects MeshTriPoints as input. Unfortunately, I don't understand how this format works and how I can convert to it from cartesian point coordinates. Here is an example for a simple cylindrical 3D surface, but I want it to be able to work also on more complex surfaces.

import numpy as np
import meshlib.mrmeshnumpy as mnp

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

# cylinder coordinates
N = 101
radius = 5

rho = radius*np.ones(N)
phi = np.linspace(0,2*np.pi,37) # 10 degree steps
z = np.linspace(-10,10,N)

# surface grid coordinates
x = np.outer(rho,np.cos(phi))
y = np.outer(rho,np.sin(phi))
z = np.outer(z,np.ones(len(phi)))

# create mesh
mesh = mnp.meshFromUVPoints(x, y, z)

# define two points on the surface
# point1
rhop = radius
phip = -10/180.*np.pi
zp = -3
xp1,yp1,zp1 = cyl2cart(rhop, phip, zp)
# point2
rhop = radius
phip = 60/180.*np.pi
zp = 8
xp2,yp2,zp2 = cyl2cart(rhop, phip, zp)

# This part is not working because I don't understand how to provide the point coordinates
# p1 = mpy.MeshTriPoint(??)
# p2 = mpy.MeshTriPoint(??)
# mpy.computeGeodesicPath(mesh, p1, p2)

#%% optional plot with plotly
import plotly.graph_objects as go
# extract numpy arrays
verts = mnp.getNumpyVerts(mesh)
faces = mnp.getNumpyFaces(mesh.topology)
# prepare data for plotly
vertsT = np.transpose(verts)
facesT = np.transpose(faces)
# draw
fig = go.Figure(data=[
    go.Mesh3d(
        x = vertsT[0],
        y = vertsT[1],
        z = vertsT[2],
        i = facesT[0],
        j = facesT[1],
        k = facesT[2]
    )
])
fig.add_trace(go.Scatter3d(x=[xp1,xp2],y=[yp1,yp2],z=[zp1,zp2],mode='markers',marker={'size':8,'symbol':'circle'}))

fig.update_layout(scene=dict(aspectmode='data'))
fig.show(renderer='browser')

Solution

  • mpy.MeshTriPoint is a point that is attached to mesh object. To get mpy.MeshTriPoint from cartesian point, project world point to mesh and find corresponding MeshTriPoint:

    from meshlib import mrmeshpy as mm
    
    mesh = mm.makeSphere( mm.SphereParams() ) # simple mesh for example
    
    startPoint = mm.Vector3f(-1,-1,-1) # cartesian start coordinate
    stopPoint = mm.Vector3f(1,1,1) # cartesian stop coordinate
    
    # find corresponding MeshTriPoints
    startMtp = mm.findProjection( startPoint, mesh).mtp
    stopMtp = mm.findProjection( stopPoint, mesh).mtp
    
    path = mm.computeGeodesicPath(mesh,startMtp,stopMtp,mm.GeodesicPathApprox.DijkstraBiDir)
    
    # print path
    for p in path:
         pc = mesh.edgePoint(p)
         print(pc.x,pc.y,pc.z)