Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
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.
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()
ax3d = fig.add_subplot(111, projection='3d')
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='-.')
axGraph.add_collection(linen)
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])
while len(q) > 0 and not found:
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='-')
ax3d.add_collection(linen)
linen = LineCollection(Vc[c1_to_c2_path,:],color=[1,0.7,0],linewidth=5,linestyle='-')
axGraph.add_collection(linen)
#%% 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:
v2e[e[0]].add(e)
v2e[e[1]].add(e[::-1])
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
visited.add(c)
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=':')
ax3d.add_collection(linen)
linen = LineCollection(Vc[predecessors_a,:],color='r',linewidth=2,linestyle=':')
axGraph.add_collection(linen)
#%% 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.
Find a longitudinal axis of the cylinder (a least-squares linear fit to all your points will yield this).
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$.
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.)
Begin a Dijkstra shortest paths traversal from $C$.
When an edge passes through the plane, color it "left" or color it "right" depending on which side it passed through.
Dijkstra inspects the neighbours of each node and a choice is made.
If a neighbour is unvisited, it is added to the frontier and given its parents colour, if any.
If a neighbour is already visited, another check is performed.
If the neighbour is uncoloured, it is rejected.
If the neighbour's colour matches the colour of the node under consideration, it is rejected.
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.
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.
