# Find shortest path around a cylinder represented by 3d triangular mesh

Suppose I have a 3d triangular mesh with the topology of a finite cylinder. Let $$C$$ be a vertex on that mesh.

How can I find the shortest path from $$C$$ to itself that goes around the cylinder? By shortest, I mean the sum of Euclidean distances between consecutive vertices on the path. By "around the cylinder" I mean that such a path would essentially divide the cylinder into two new cylinders.

I have also asked this question on Mathematics SE.

• A path cannot divide a 3D object into two halves, it has to be a surface. – Maxim Umansky Oct 8 at 15:28
• It's simple: Due to the fact that circles that their center lies on the axis of the cylinder and are in the planes that their normal vector are parallel to the cylinder's axis are the geodesic curves of the cylinder, the shortest path from C to itself is a circle that pass the C and it's center is in the cylinder axis and the normal of vector of the plane that contains the circle is parallel to the cylinder axis. – Alone Programmer Oct 8 at 15:45
• If we need to find the shortest path between grid nodes on a 2D manifold (say the surface of the cylinder) then this seems to be the case for the Shortest Path on a Graph problem, there are several algorithms for it, en.wikipedia.org/wiki/Shortest_path_problem. – Maxim Umansky Oct 9 at 3:40
• I think the tricky part here is that its not just a travelling-salesman-problem or shortest-path problem, but that there is a constraint. The hard part will be to find a viable definition of the constraint! "Goes around once" may be nontrivial to code:-) – MPIchael Oct 9 at 8:26

OK, after thinking about it for a while, I came up with an answer.

Step 1: Find the caps of the cylinder, in other words two closed disjoint paths along the graph's borders.

Step 2: Find a path along the face graph from one cap to the other.

Step 3: Create a new sub-graph by removing all edges which lie along the path found in step 2. Keep track of the edges removed as they will be used later. Since the path is from cap to cap, it cuts the cylinder so that the resulting sub-graph has a planar topology.

Step 4: Use the Dijkstra algorithm to find the shortest path from $$C$$ to every other point on the sub-graph from step 3.

Step 5: Let $$e$$ be an edge among the edges removed in step 3. Let $$p_1$$ and $$p_2$$ be the vertices connected by that edge. Let $$D(C,p_1)$$ and $$D(C,p_2)$$ be the distances along the shortest paths found in step 4, from $$C$$ to $$p_1$$ and $$p_2$$ respectively. Let $$D(e)$$ be the length of edge $$e$$. Among the edges removed in step 3, find the one that minimizes $$D(C,p_1) + D(C,p_2) + D(e)$$

The shortest path is obtained by joining the shortest path from $$C$$ to $$p_1$$ as found in step 4, the edge $$e$$ and then the shortest path from $$p_2$$ to $$C$$.

I spent the last day writing an example in python. The example first generates a graph with cylindrical topology, and then finds the shortest path from C to itself around the cylinder according to the algorithm described:


import numpy as np

from matplotlib import pyplot as plt
from mpl_toolkits import mplot3d

debugPlot = True
forceEdgePoint = False

#%% Generate random points and duplicate them along the y axis
N = 100
V_orig = np.random.rand(N,2)

V = np.concatenate((V_orig,V_orig+np.array([0,1]),V_orig+np.array([0,2])),axis=0)

from scipy import spatial

#%% Perform delaunay triangulation
tri = spatial.Delaunay(V)

if debugPlot:
plt.figure()
plt.triplot(V[:,0],V[:,1],tri.simplices.copy(),'.-b')
plt.gca().set_aspect('equal')

#%% Create a cylinder topology on the points in the middle
S = tri.simplices.copy()
#remove everything except for the middle:
S = S[np.any((S>=N) & (S < 2*N),axis=1),:]
if debugPlot:
plt.triplot(V[:,0],V[:,1],S,'+-r')
plt.gca().set_aspect('equal')

Sc = np.mod(S,N)
#remove duplicates
imin = np.argmin(Sc,axis=1)
reindex = np.vstack((imin,imin+1,imin+2)).transpose() % 3
Sc = Sc[np.arange(Sc.shape[0]).reshape(-1,1),reindex]
Sc = np.unique(Sc,axis=0)
#remove degenerates
Sc = Sc[~((Sc[:,0] == Sc[:,1]) | (Sc[:,1] == Sc[:,2]) | (Sc[:,0] == Sc[:,2])),:]
Vc = V[:N,:]
plt.figure()
if debugPlot:
plt.triplot(Vc[:,0],Vc[:,1],Sc,'.-g',linewidth=5)
axGraph = plt.gca()

#%% Remove the cylinder caps
u = np.sort(Vc[:,0])
capsLowerBound = u[int(len(u) * 0.1)]
capsUpperBound = u[int(len(u) * 0.9)]
caps = (Vc[:,0] < capsLowerBound) | (Vc[:,0] > capsUpperBound)
caps_i = np.nonzero(caps)[0]
#Vc = Vc[~caps,:]
Sc = Sc[~np.any(np.isin(Sc,caps_i),axis=1),:]
plt.triplot(Vc[:,0],Vc[:,1],Sc,'o-b',linewidth=3)
axGraph = plt.gca()

#%% Plot in 3d a cylinder generated from the 2d coordinates generated above
Z = Vc[:,0]
X = np.cos(Vc[:,1]*2*np.pi)
Y = np.sin(Vc[:,1]*2*np.pi)
XYZ = np.vstack((X,Y,Z)).transpose()
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
h = ax3d.plot_trisurf(X,Y,Sc,Z)
#h.set_facecolor(None)
h.set_edgecolor('k')

#%% Prepare halfedges
halfEdges = np.concatenate((np.array([Sc[:,0],Sc[:,1]]),
np.array([Sc[:,1],Sc[:,2]]),
np.array([Sc[:,2],Sc[:,0]])),axis=1).transpose()
halfEdge2Tri = {tuple(x):(y % Sc.shape[0]) for y,x in enumerate(halfEdges)}

#%% find non unique half-edges (there should be none)
iSEdges = np.lexsort(np.flipud(halfEdges.transpose()))
halfEdges[iSEdges,:]
iNonUnique = np.nonzero(np.all(np.diff(halfEdges[iSEdges,:],axis=0)==0,axis=1))[0]
nnuEdges = halfEdges[iSEdges,:][iNonUnique,:]
nnuVc = Vc[nnuEdges,:]
from matplotlib.collections import LineCollection
linen = LineCollection(nnuVc,color='r',linewidth=2,linestyle='-.')
if len(nnuEdges)>0:
raise Exception('Found non unique half-edges')
halfEdgesSet = set(list(tuple(x) for x in halfEdges))

#%% Find borders of the cylinder (caps)
borderHalfEdges = {x for x in halfEdgesSet if x[::-1] not in halfEdgesSet}
borderHalfEdgePerVertex = dict()
for e in borderHalfEdges:
borderHalfEdgePerVertex[e[0]] = e

def FindClosedBorderPath(halfEdgeSubset):
try:
e = next(iter(halfEdgeSubset))
except StopIteration:
return []
capPath = [e[0],e[1]]
halfEdgeSubset.remove(e)
while capPath[-1] != capPath[0]:
e = borderHalfEdgePerVertex[capPath[-1]]
capPath.append(e[1])
halfEdgeSubset.remove(e)
return capPath

capPath1 = FindClosedBorderPath(borderHalfEdges)
capPath2 = FindClosedBorderPath(borderHalfEdges)
capPath3 = FindClosedBorderPath(borderHalfEdges)
if len(capPath3) != 0:
raise Exception("For some reason there are 3 borders")

#%% Find a path from  one cap to the other along the face graph
c1 = (capPath1[0],capPath1[1])
c2 = (capPath2[0],capPath2[1])
q = [c1[::-1]]
predecessor = dict()
found = False

def GenEdgesForTri(t):
yield (t[0],t[1])
yield (t[1],t[2])
yield (t[2],t[0])

c = q.pop(0)
try:
t = Sc[halfEdge2Tri[c[::-1]]]
except KeyError:
continue
for e in GenEdgesForTri(t):
if e not in predecessor:
q.append(e)
predecessor[e] = c
if e == c2:
found = True
break
assert(found)

c1_to_c2_path = [c2]
n = c2
while n!=c1[::-1]:
n = predecessor[n]
c1_to_c2_path.append(n)
c1_to_c2_path.reverse()

#%%plot the newly found path
if debugPlot:
c1_c2_path_xyz = XYZ[c1_to_c2_path,:]
linen = mplot3d.art3d.Line3DCollection(XYZ[c1_to_c2_path,:],color=[1,0.7,0],linewidth=5,linestyle='-')
linen = LineCollection(Vc[c1_to_c2_path,:],color=[1,0.7,0],linewidth=5,linestyle='-')

#%% Remove all triangles that are part of the path and save it as a new subgraph
subgraphHalfEdges = set(halfEdgesSet) #copy
removedEdges=set(c1_to_c2_path) | {x[::-1] for x in c1_to_c2_path}
for e in removedEdges:
try:
subgraphHalfEdges.remove(e)
except KeyError:
pass
#%% Choose C
if forceEdgePoint:
iSorted = np.argsort(Vc[Sc.flat,0])
iC = Sc.flat[iSorted[0]]
else:
iC = np.random.choice(Sc.flat)
ax3d.plot([X[iC]],[Y[iC]],[Z[iC]],'*m',markersize=30)
axGraph.plot([V[iC,0]],[V[iC,1]],'*m',markersize = 20)

#%% Make Dijkstra on subgraph
import bisect

def DistFunc(e):
d = XYZ[e[0],:] - XYZ[e[1],:]
return np.sqrt(d.dot(d))

from collections import defaultdict
v2e = defaultdict(lambda : set())
for e in subgraphHalfEdges:

distances = defaultdict(lambda : np.inf)
visited = set()
predecessor = {}
q = [(0,iC)]
distances[iC]=0
while len(q) > 0:
c = q.pop(0)[1]
if c in visited:
continue
thisDist = distances[c]
for e in v2e[c]:
v = e[1]
if v in visited:
continue
lenE = DistFunc(e)
vCurDist = distances[v]
if vCurDist > thisDist + lenE:
distances[v] = thisDist + lenE
predecessor[v] = c
bisect.insort(q,(thisDist+lenE,v))
#%% plot the result
if debugPlot:
l=list(distances.items())
xyz = XYZ[list(x[0] for x in l),:]
d = np.array([x[1] for x in l])
c = d / np.max(d)
ax3d.scatter(xyz[:,0],xyz[:,1],c=c,zs=xyz[:,2])
predecessors_a = np.array(list(predecessor.items()))
linen = mplot3d.art3d.Line3DCollection(XYZ[predecessors_a,:],color='w',linewidth=2,linestyle=':')
linen = LineCollection(Vc[predecessors_a,:],color='r',linewidth=2,linestyle=':')

#%% Find the path through the removed edges that would make the shortest total path
candidates = []
for e in removedEdges:
d = DistFunc(e)
candidates.append( (d + distances[e[0]] + distances[e[1]], e))
candidates.sort()
chosenE = candidates[0][1]

subpaths = []
for i in range(2):
n = chosenE[i]
path_C_to_p = [n]
while n!=iC:
n = predecessor[n]
path_C_to_p.append(n)
subpaths.append(path_C_to_p)

shortestPath = list(reversed(subpaths[0])) + subpaths[1]
shortestPathXyz = XYZ[shortestPath,:]
shortestPathV = Vc[shortestPath,:]
axGraph.plot(shortestPathV[:,0],shortestPathV[:,1],'-c',linewidth=3)
ax3d.plot(shortestPathXyz[:,0],shortestPathXyz[:,1],shortestPathXyz[:,2],'-c',linewidth=3)


In the following figures, the orange lines are the edges removed in step 3, the magenta star is $$C$$, and the cyan line is the shortest path. The 3d embedding of the graph as a cylinder was made for clarity, it is in no way necessary to solve the problem, although the edge lengths are taken from that 3d embedding.

• Mapping the flood-fill back to a graph representation is the challenge here. – Richard Oct 9 at 16:04
• @Richard How about if I perform the flood fill on the face graph? – iliar Oct 9 at 16:25
• I don't see how you determine where the seam/cut is that the flood-fill makes. The frontier of the flood-fill is always meeting itself. I think what worries me is that this approach is tied solely to local properties of the mesh, but tries to extract a non-local quantity. It's difficult for me to intuit that there isn't an irregular mesh that would break things. – Richard Oct 9 at 16:33
• @Richard I'm not sure, perhaps you're right so I modified the algorithm to use a path in the face graph from one cylinder cap to the other, this one is sure to cut the cylinder in two along the 'axis'. – iliar Oct 10 at 18:59
1. Find a longitudinal axis of the cylinder (a least-squares linear fit to all your points will yield this).

2. Construct a plane passing through this axis. Any orientation should be fine, but let's say it is perpendicular to a normal passing from the axis to $$C$$.

3. Reorientate so the axis is vertical and the plane is divided into "left" and "right" halves by the axis. (This step is more conceptual than mathematical.)

4. Begin a Dijkstra shortest paths traversal from $$C$$.

5. When an edge passes through the plane, color it "left" or color it "right" depending on which side it passed through.

6. Dijkstra inspects the neighbours of each node and a choice is made.

1. If a neighbour is unvisited, it is added to the frontier and given its parents colour, if any.

2. If a neighbour is already visited, another check is performed.

1. If the neighbour is uncoloured, it is rejected.

2. If the neighbour's colour matches the colour of the node under consideration, it is rejected.

3. If the neighbour's colour differs from the node under consideration, then the edge between the node under consideration and this neighbour completes the shortest circumferential path.

7. Use parent pointers to backtrace the path.

The left-right plane intersections should be simple geometric checks and you shouldn't need to worry about floating-point issues since it doesn't matter if the checks are perfect or not: there's a quarter-cylinder worth of wiggle room.

• Your answer would be good in case the object in question is an actual cylinder. Perhaps I wasn't clear, but I want an answer that's independent of the actual geometry, as the object in question merely has the topology of a cylinder but its geometry may be distorted. – iliar Oct 10 at 18:33
• @iliar: Can you clarify "distorted"? Can the cylinder be coiled up like a snake, made toroidal, &c? – Richard Oct 10 at 21:35
• I think it's enough for it to be wide and bent or skewed. – iliar Oct 10 at 22:44