 # Straight Skeleton - or how to draw a center line in a polygon or shape?

Hello,

for a project I want to draw a “center line” into a shape/polygon. The picture below shows the letter “H” outlined (this is the data I have) and then the center line in red (this is the data I am looking for). I am looking for a way to calculate the red line. My code stores the coordinates of the shape (black line) in a PShape. So what I am looking for is a function or library which takes the PShape and returns the red line as PShape or even as a list of coordinates.

I found many links in google to the topic “Straight Skeleton” but I am wondering if someone have had a similar challenge in the past and found a solution in Processing already. Any help would be much appreciated.

Thanks
Mike

1 Like

Hi @Maxs,

Quick questions:

• Are you working with fonts only or with any kind of shapes ?
• Are you using Processing Java or Python mode ?

Quick answers: Computing the straight skeleton of a polygon is not trivial and you probably want to use a library for that. The problem is that there are very few. In Python mode I am using Polyskel, for Java you could try CampSkeleton (haven’t tested this one).

2 Likes

Thanks for your fast reply. I am working mostly with fonts. However some of the fonts are quite complex. I am using Processing in Java mode. I will check the libraries you mentioned.

If I would use fonts only would that change anything regarding the solution? Maybe with the geomerative lib?

If (and only if) your fonts are `sans serif` (i.e inner and outer edges are always equally spaced) you could offset (not scale) a letter inward until it forms a line (a very thin polygon to be exact).

Regarding “Geomerative”, I think the library gives you access to the inner and outer edges of a font (its `path`) but not to its skeleton.

If it’s just a shape and you have it’s coordinates then you’re going to want to calculate the points (unknown) that fall exactly between any two other points (known). This will be very easy for rectangular shapes.

Then you can attempt to determine the space inside the shape given an offset and possibly draw the line.

For this ‘H’ there are three rectangles you should discover algorithmically and with then you could draw that red interior.

P.S.

If you want do it for less uniform shapes you’ll have to generalize a big. Calculating which points are closest to each other on the same axis (x or y line).

Thanks solub and Nathan for your excellent input. I will work on this in the next few days. If I come up with something I will share here. @solub, the idea to offset is very interesting and a good workaround for the time being. @Nathan thanks for your response, especially the last sentence is very helpful and resonates most.

@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.x, seg.y, seg.x, seg.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:

#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.x , b.y, b.x , b.y)
popStyle()

#Draw Medial-Axis
pushStyle()
stroke(255, 70, 122)
strokeWeight(1.8)
for b in bones:
line(b.x , b.y, b.x , b.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, cr, r, r)
if x:
d = x.dist(cr)
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) #rays sorted by distance from their s-pt to their intersection pt
selectedrays =  [r for r in srays if r == srays] #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, r2):
if r not in self.vertices:
d1 = X.dist(r1)
d2 = X.dist(r2)
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, r1, r2, r2) #s-pt & e-pt of 1st edge, s-pt & e-pt of 2nd edge

#check if ray points in the right direction
si = self.vertices.index(r2) #s-pt index
ei = self.vertices.index(r1) #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, r2)
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, self.rays))

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.

1 Like

Dear solub,

first of all, sorry! Sorry that I did not come back to you, especially seeing your excellent research. For whatever reason I did not get any notification that there was an update to my post. I will need a bit of time to read thru it… and try to implement option 1. I will come back to you. Thank you so much.
Mike

Hi @Maxs, no problem !

I went ahead and tried both methods with font shapes.

It turns out the first method is likely to be the least suitable in your case. Two main reasons for this:

• Smooth center lines require a lot of sample points along the contour of a letter. As a result, computing a Voronoi diagram for a whole word becomes very expensive.

• Branch-like segments can not be removed. The center line will always be attached to some sort of connection lines

Test Voronoi
``````add_library('geomerative')

def setup():
size(800, 800, P2D)
background('#FFFFFF')
strokeWeight(1.8)
noFill()
smooth(8)

RG.init(this)

#Select word + font + size then draw
txt = RG.getText("T", "Roboto-Black.ttf", 800, CENTER) # <-- put the font file inside the sketch folder
txt.translate(width/2, height/1.17)
txt.draw()

#Convert sample points (along the contour of the font) to PVectors
vertices = [PVector(p.x, p.y) for p in txt.getPoints()]

#Feed these points to the algorithm to retrieve the medial axis
axis = MedialAxis(vertices)
segments = axis.segments

#Draw the segments forming the medial axis
stroke(255, 70, 122)
for seg in segments:
line(seg.x, seg.y, seg.x, seg.y)

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:

#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]):
e = {v1, v2}
if e not in self.edges:
self.edges.append(e)

@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:
if self.contain(self.vertices, p1) and self.contain(self.vertices, p2):
for v1, v2 in zip(self.vertices, self.vertices[1:] + self.vertices[:1]):
if self.intersect(p1, p2, v1, v2):
break
else:
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'''

try:
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))

except:
return

if uA >= 0 and uA <= 1 and uB >= 0 and uB <= 1:

return 1
``````

On the contrary the Straight-Skeleton approach doesn’t require lots of vertices and is much faster. However:

• You will need to treat fonts with inner contours (like the letter ‘O’ or ‘B’) as polygons with holes. This is a bit of a problem since:

• Geomerative doesn’t make a distinction between inner and outer vertices (all points come in the same array). As a result you will have to split the array manually for each letter.

• Geomerative doesn’t provide the main vertices of a font. You will need to retrieve them by culling the control points that are not on its contour.

• not all S-Skeleton libraries handle polygons with holes.

• Fonts have a lot of nearly-parallel lines and most S-Skeleton libraries don’t handle well shapes with symmetrical edges. You may need to add some noise to the vertices in that case.

• Most S-Skeleton libraries don’t make the distinction between “branches” and inner “bones” (the center line). You may need to remove these branches afterwards with a line-line collision algorithm.

• The center line will not be as smooth as you would expect.

1 Like

Thanks solub,

I you rock! It is so pleasurable to read your reply and I am sure that others benefit from it, too. In fact the reason why I need this at all is, that I want to generate gCode (CNC) so that a tool can “write” the letters with just one stroke. I found some TTF (fonts) which are very thin, however after scaling there is still a chance that the tool passes some line twice.

I am so thankful, that you provide python code and since I am using processing in java mode I have to think about calculating the Straight-Skeleton “offline” and then import the result. For my application there is no time criticality so handing over the data to Python and then read the results in processing is really a way to go.

I have to go on a trip the next few days and I will definitely take your script into the plane. Thanks for your exceptional support. I hope I can come up with some results and share with you. Mike

1 Like

If you are looking to carve text with a cnc machine, you may wish to consider “Hershey Fonts”. It is not related to chocolate bars. The fonts are easy to use lines for computers. Do a quick look at https://en.wikipedia.org/wiki/Hershey_fonts
Search the internet for software which has already been written to read these fonts. This likely already exists for many things.