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.

