@Maxs How is your project going ? Have you managed to retrieve the medial axis from your fonts ?
I have been doing some digging and came up with 2 different approaches that could be of interest to you. The first one (based on this SO thread) is certainly the simplest and consists in computing a Voronoi diagram from a polygon’s vertices and isolate the edges whose both endpoints lie within the polygon.
Algorithm:
- compute the Voronoi diagram from the vertices of a chosen polygon (a “letter” in your case)
- remove the duplicate Voronoi edges (adjacent Voronoi cells have similar edges)
- find all the Voronoi edges whose both endpoints are within the polygon
- from this selection, remove any Voronoi edge that intersect a polygon edge (important!)
You will end-up with an array of edges (or segments) that, stuck together, form a kind of center line.
Below an implementation in Processing Python using the Voronoi() class from the Toxiclibs library.
Medial Axis from Voronoi diagram
add_library('toxiclibs')
verts = [PVector(901.12, 300),PVector(853.48, 279.76),PVector(768.04, 270.04),PVector(701.56, 261.04),
PVector(684.4, 253.72),PVector(684.4, 250.48),PVector(701.32, 234.76),PVector(701.2, 231.28),
PVector(685.84, 223.48),PVector(631.12, 220.6),PVector(570.52, 218.32),PVector(534.88, 200.68),
PVector(434.2, 50.56),PVector(362.68, 50.56),PVector(320.32, 218.44),PVector(306.16, 213.52),
PVector(258.28, 140.68),PVector(191.44, 140.68),PVector(164.8, 253.12),PVector(233.68, 262.96),
PVector(245.92, 265.12),PVector(252.64, 268.72),PVector(278.2, 268.12),PVector(265.0, 276.04),
PVector(244.24, 284.68), PVector(244.24, 315.32), PVector(265.0, 323.96),PVector(278.2, 331.88),
PVector(252.64, 331.28),PVector(245.92, 334.88),PVector(233.68, 337.04),PVector(164.8, 346.88),
PVector(191.44, 459.32),PVector(258.28, 459.32),PVector(306.16, 386.47998),PVector(320.32, 381.56),
PVector(362.68, 549.44),PVector(434.2, 549.44),PVector(534.88, 399.32),PVector(570.52, 381.68),
PVector(631.12, 379.4),PVector(685.84, 376.52002),PVector(701.2, 368.72),PVector(701.32, 365.24),
PVector(684.4, 349.52002),PVector(684.4, 346.28),PVector(701.56, 338.96),PVector(768.04, 329.96),
PVector(853.48, 320.24)]
def setup():
size(1000, 600, P2D)
background('#FFFFFF')
smooth(8)
axis = MedialAxis(verts)
segments = axis.segments
edges = set(axis.edges).difference(segments)
#Draw polygon's silhouette
pushStyle()
stroke(50, 100, 255)
strokeWeight(3)
for p1, p2 in zip(verts, verts[1:] + verts[:1]):
line(p1.x, p1.y, p2.x, p2.y)
popStyle()
#Draw vertices
pushStyle()
stroke(50)
strokeWeight(6)
for i, p in enumerate(verts):
point(p.x, p.y)
popStyle()
#Draw Voronoi cells
pushStyle()
strokeWeight(.3)
stroke(40, 235, 180)
for p1, p2 in edges:
line(p1.x, p1.y, p2.x, p2.y)
popStyle()
#Draw Medial-Axis
pushStyle()
stroke(255, 70, 122)
strokeWeight(1.8)
for seg in segments:
line(seg[0].x, seg[0].y, seg[1].x, seg[1].y)
popStyle()
noLoop()
class MedialAxis(object):
def __init__(self, verts):
self.vertices = verts
self.segments = []
self.edges = []
self.run()
def run(self):
self.initialize
self.cull
@property
def initialize(self):
'''Computes a Voronoi diagram.
Stores the vertices of each region in a separate sub-list.'''
#instantiate Toxiclibs Voronoi() class
voronoi = Voronoi()
#convert vertices to Vec2D format
for v in self.vertices:
voronoi.addPoint(Vec2D(v.x, v.y))
#store unique edges of the Voronoi diagram
for i, poly in enumerate(voronoi.getRegions()):
#convert vertices from Vec2D() to PVector()
vertices = [PVector(v.x(), v.y()) for v in poly.vertices]
for v1, v2 in zip(vertices, vertices[1:] + vertices[:1]):
if (v1, v2) not in self.edges and (v2, v1) not in self.edges:
self.edges.append((v1, v2))
@property
def cull(self):
'''Finds medial-axis segments.
Stores Voronoi cell's edge if it lies entirely inside the main polygon.'''
for p1, p2 in self.edges:
isValid = False
if self.contain(verts, p1) and self.contain(verts, p2):
isValid = True
for v1, v2 in zip(self.vertices, self.vertices[1:] + self.vertices[:1]):
if self.intersect(p1, p2, v1, v2):
isValid = False
break
if isValid:
self.segments.append((p1, p2))
def contain(self, verts, pt):
'''Returns a boolean.
Determines whether a point lies inside a shape/polygon or not.'''
N = len(verts)
isInside = 0
for i in xrange(N):
v1 = verts[(i+1)%N]
v2 = verts[i]
if (v2.y > pt.y) != (v1.y > pt.y):
if pt.x < (v1.x - v2.x) * (pt.y - v2.y) / (v1.y - v2.y) + v2.x:
isInside = not isInside
return isInside
def intersect(self, p1, p2, p3, p4):
'''Returns a Boolean.
Checks if 2 lines are intersecting. Option: returns location of intersection point'''
uA = ((p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y+1))
uB = ((p2.x-p1.x)*(p1.y-p3.y) - (p2.y-p1.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y+1))
if uA >= 0 and uA <= 1 and uB >= 0 and uB <= 1:
return 1
The second approach, and undoubtedly the most the complex one, consists in computing the Straight Skeleton of a polygon and remove the “branches” that connect the inner bones to the vertices of that polygon.
Algorithm:
- calculate the bisecting vector at each vertex of a chosen polygon. Make it a ray.
- find the intersection points with the previous and the next neighboring ray.
- find all pairs of rays that intersect sooner than either intersects its other neighbor.
- sort them by distance to their intersection point.
- take the pair with the shortest distance and replace it with a single new ray originating from their intersection point. That child ray must be aiming towards the average direction of its neighboring polygon edges.
- repeat from step 2 until there are only 2 rays left.
- connect the originating points of the last 2 rays to complete the skeleton
Below a naive implementation in Processing Python:
Medial Axis from Straight Skeleton
verts = [PVector(901.12, 300),PVector(853.48, 279.76),PVector(768.04, 270.04),PVector(701.56, 261.04),
PVector(684.4, 253.72),PVector(684.4, 250.48),PVector(701.32, 234.76),PVector(701.2, 231.28),
PVector(685.84, 223.48),PVector(631.12, 220.6),PVector(570.52, 218.32),PVector(534.88, 200.68),
PVector(434.2, 50.56),PVector(362.68, 50.56),PVector(320.32, 218.44),PVector(306.16, 213.52),
PVector(258.28, 140.68),PVector(191.44, 140.68),PVector(164.8, 253.12),PVector(233.68, 262.96),
PVector(245.92, 265.12),PVector(252.64, 268.72),PVector(278.2, 268.12),PVector(265.0, 276.04),
PVector(244.24, 284.68), PVector(244.24, 315.32), PVector(265.0, 323.96),PVector(278.2, 331.88),
PVector(252.64, 331.28),PVector(245.92, 334.88),PVector(233.68, 337.04),PVector(164.8, 346.88),
PVector(191.44, 459.32),PVector(258.28, 459.32),PVector(306.16, 386.47998),PVector(320.32, 381.56),
PVector(362.68, 549.44),PVector(434.2, 549.44),PVector(534.88, 399.32),PVector(570.52, 381.68),
PVector(631.12, 379.4),PVector(685.84, 376.52002),PVector(701.2, 368.72),PVector(701.32, 365.24),
PVector(684.4, 349.52002),PVector(684.4, 346.28),PVector(701.56, 338.96),PVector(768.04, 329.96),
PVector(853.48, 320.24)]
def setup():
size(1000, 600, P2D)
background('#FFFFFF')
smooth(8)
skeleton = Skeleton(verts, 20) #(vertices, threshold)
branches = skeleton.branches
bones = skeleton.bones
#Draw polygon's silhouette
pushStyle()
stroke(50, 100, 255)
strokeWeight(3)
for p1, p2 in zip(verts, verts[1:] + verts[:1]):
line(p1.x, p1.y, p2.x, p2.y)
popStyle()
#Draw vertices
pushStyle()
stroke(30)
strokeWeight(6)
for i, p in enumerate(verts):
point(p.x, p.y)
popStyle()
#Draw Branches
pushStyle()
strokeWeight(.3)
stroke(40, 235, 180)
for i, b in enumerate(branches):
line(b[0].x , b[0].y, b[1].x , b[1].y)
popStyle()
#Draw Medial-Axis
pushStyle()
stroke(255, 70, 122)
strokeWeight(1.8)
for b in bones:
line(b[0].x , b[0].y, b[1].x , b[1].y)
popStyle()
noLoop()
class Skeleton(object):
def __init__(self, verts, tol = 0):
self.vertices = verts
self.branches = []
self.bones = []
self.rays = []
self.tol = tol
self.run()
def run(self):
'''Runs the algorithm.'''
self.initialize
self.reduce
@property
def initialize(self):
'''Starts by computing the rays originating from the corners of the polygon.
Each ray object also stores the coordinates of its neighboring polygon edges (previous and next).'''
lv = len(self.vertices)
for i in xrange(lv):
pv = self.vertices[(i-1)%lv] #previous vertex
cv = self.vertices[i] #current vertex
nv = self.vertices[(i+1)%lv] #next vertex
b = self.bisector(pv, cv, cv, nv)
ept = b.setMag(100000).add(cv) #end-point of the ray (far from start-point)
self.rays.append((cv, ept, (pv, cv), (cv, nv))) #start-point, end-point, previous edge, next edge
@property
def reduce(self):
'''Finds consecutive rays intersecting each other before either intersects its other neighbor.
Replaces these pairs of rays by a single one originating from their intersection point and
aiming towards the average direction of its neighboring polygon edges.'''
while len(self.rays) > 2:
intersections = [] # Nearest intersection points from parent
lr = len(self.rays)
for i in xrange(lr):
pr = self.rays[(i-1)%lr] #previous ray
cr = self.rays[i] #current ray
nr = self.rays[(i+1)%lr] #next ray
mind, X = 100000, None #min-distance, intersection point
for r in (pr, nr):
x = self.intersect(cr[0], cr[1], r[0], r[1])
if x:
d = x.dist(cr[0])
if d < mind:
mind = d
X = x
intersections.append((X, mind))
candidates = []
for id, ((X1, d1), (X2, d2)) in enumerate(zip(intersections, intersections[1:]+intersections[:1])):
if (X1 and X2) != None and X1 == X2:
dsum = d1 + d2 #sum of distances from selected rays (s-pt) to their intersection point
candidates.append((id, X1, dsum))
if not candidates:
#making sure that NO consecutive rays intersect each other before either intersects its other neighbor
if len(set(intersections)) == len(intersections):
for id, ((X1, d1), (X2, d2)) in enumerate(zip(intersections, intersections[1:]+intersections[:1])):
if (X1 and X2) != None:
#sending to candidates anyway
dsum = d1 + d2 #sum of distances from selected rays (s-pt) to their intersection point
candidates.append((id, X1, dsum))
srays = sorted(candidates, key=lambda t: t[2]) #rays sorted by distance from their s-pt to their intersection pt
selectedrays = [r for r in srays if r[2] == srays[0][2]] #in cases when several rays are equally distant from their intersection pt
if selectedrays:
offset = 0
for ray in selectedrays:
id, X, _ = ray
r1 = self.rays[id-offset]
r2 = self.rays[(id+1-offset)%lr]
#stores bones (segments from parent rays to intersection point 'X')
for r in (r1[0], r2[0]):
if r not in self.vertices:
d1 = X.dist(r1[2][1])
d2 = X.dist(r2[3][0])
if (d1 + d2) / 2 > self.tol:
self.bones.append((r, X))
else:
self.branches.append((r, X))
else:
self.branches.append((r, X))
#compute direction of the new ray
b = self.bisector(r1[2][0], r1[2][1], r2[3][0], r2[3][1]) #s-pt & e-pt of 1st edge, s-pt & e-pt of 2nd edge
ept = X.copy().add(b.setMag(100000))
#check if ray points in the right direction
si = self.vertices.index(r2[3][0]) #s-pt index
ei = self.vertices.index(r1[2][1]) #e-pt index
slice = self.vertices[si:ei] if ei > si else self.vertices[si:] + self.vertices[:ei]
if not self.contain([X] + slice, X.copy().add(b.setMag(1))):
ept.mult(-1)
#delete parents rays from array list and insert their 'child' instead
self.rays[id-offset] = (X, ept, r1[2], r2[3])
del self.rays[(id+1-offset)%lr]
offset += 1
lr = len(self.rays)
else:
print "Error: no rays have been found for reduction. A shared intersection point is probably missing."
break
#connect start-points of the last 2 rays
self.bones.append((self.rays[0][0], self.rays[1][0]))
def bisector(self, p1, p2, p3, p4):
'''Returns a PVector.
Computes the bisector of a corner between edge p1-p2 and edge p3-p4.'''
dir1 = (p2 - p1).normalize() #direction between 1st point and 2nd point of edge 1
dir2 = (p4 - p3).normalize() #direction between 1st point and 2nd point of edge 2
dsum = dir1 + dir2
return PVector(dsum.y, -dsum.x).normalize()
def intersect(self, p1, p2, p3, p4):
'''Returns a Boolean.
Checks if 2 lines are intersecting. Option: returns location of intersection point'''
uA = ((p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y))
uB = ((p2.x-p1.x)*(p1.y-p3.y) - (p2.y-p1.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y))
if uA >= 0 and uA <= 1 and uB >= 0 and uB <= 1:
secX = p1.x + (uA * (p2.x-p1.x))
secY = p1.y + (uA * (p2.y-p1.y))
return PVector(round(secX), round(secY))
def contain(self, verts, pt):
'''Returns a boolean.
Determine whether a point lies inside a shape/polygon or not.'''
N = len(verts)
isInside = 0
for i in xrange(N):
v1 = verts[(i+1)%N]
v2 = verts[i]
if (v2.y > pt.y) != (v1.y > pt.y):
if pt.x < (v1.x - v2.x) * (pt.y - v2.y) / (v1.y - v2.y) + v2.x:
isInside = not isInside
return isInside
Please note that this version is made simple for clarity and doesn’t address all the special cases (polygon with holes, with a large number of nearly-parallel edges, those that exhibit symmetrical conditions…) that a more robust library would probably handle with ease.
Yet this is probably one of the few implementations that makes possible to select the bones above a chosen threshold distance ( tol in the code) from the polygon’s vertices. These inner bones are the very segments that form the medial axis.
.
