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

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.

1 Like

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.

3 Likes

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.

3 Likes

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

2 Likes

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.
Inkscape has Hershey Text which can be changed to svg which might do what you need.

I just ported the Medial Axis from Straight Skeleton to Java.

Code
``````import static processing.core.PApplet.round;

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.stream.Stream;

import processing.core.PVector;

/**
* Algorithm:
* <p>
* 1. calculate the bisecting vector at each vertex of a chosen polygon. Make it
* a ray.<br>
* 2. find the intersection points with the previous and the next neighboring
* ray.<br>
* 3. find all pairs of rays that intersect sooner than either intersects its
* other neighbor.<br>
* 4. sort them by distance to their intersection point.<br>
* 5. 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.<br>
* 6. repeat from step 2 until there are only 2 rays left.<br>
* 7. connect the originating points of the last 2 rays to complete the
* skeleton<br>
*
*/
public class SolubSkeleton {

/**
* https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/7
*/

final float tol;

public ArrayList<PVector> vertices;
public ArrayList<Ray> rays;
public ArrayList<Bone> bones;
public ArrayList<Bone> branches;

public SolubSkeleton(ArrayList<PVector> vertices, float tol) {
this.vertices = vertices;
this.tol = tol;
rays = new ArrayList<>();
bones = new ArrayList<>();
branches = new ArrayList<>();
}

public void run() {
init();
reduce();
}

/**
* Starts by computing the rays originating from the corners of the polygon.
* Each ray object also stores the PVectors of its neighboring polygon edges
* (previous and next).
*/
private void init() {
int lv = vertices.size();

for (int i = 0; i < vertices.size(); i++) {
PVector pv = vertices.get(Math.floorMod(i - 1, lv)); // previous vertex
PVector cv = vertices.get(i); // current vertex
PVector nv = vertices.get((i + 1) % lv); // next vertex

PVector b = bisector(pv, cv, cv, nv);
PVector ept = b.setMag(100000).add(cv); // end-point of the ray (far from start-point)

Ray r = new Ray(cv, ept, new Edge(pv, cv), new Edge(cv, nv));
}

}

private void reduce() {

/*
* Check for ray-ray collision
*/
while (rays.size() > 2) {

ArrayList<Intersection> intersections = new ArrayList<>(); // Nearest intersection points from parent
int lr = rays.size();

for (int i = 0; i < lr; i++) {
Ray pRay = rays.get(Math.floorMod(i - 1, lr)); // previous ray
Ray cr = rays.get(i); // current ray
Ray nRay = rays.get((i + 1) % lr); // next ray

float minD = 100000; // min-distance
PVector X = null; // intersection point (optional)

// check current ray for intersection with prev and next

PVector x = intersect(cr, pRay); // check with prev
if (x != null) { // rays intersected
float d = x.dist(cr.startPoint);
if (d < minD) {
minD = d;
X = x;
}
}

x = intersect(cr, nRay); // check with next
if (x != null) { // rays intersected
float d = x.dist(cr.startPoint);
if (d < minD) {
minD = d;
X = x;
}
}

}

ArrayList<Candidate> candidates = new ArrayList<>();

for (int i = 0; i < intersections.size(); i++) {

/**
* zip(intersections, intersections[1:]+intersections[:1]) is equivalent to
* iterating over consecutive pairs of Iterations
*/

Intersection i1 = intersections.get(i);
Intersection i2 = intersections.get((i + 1) % intersections.size()); // wrap around to 0 for last pair
if ((i1.X() != null && i2.X() != null)) { // if D not null and equivalent

if (i1.X().equals(i2.X())) {
// sum of distances from selected rays (s-pt) to their intersection point
float dSum = i1.minD + i2.minD;
}
}
}

if (candidates.isEmpty()) {
// making sure that NO consecutive rays intersect each other before either
// intersects its other neighbor
if (new HashSet<>(intersections).size() == intersections.size()) {
// similar loop to above
for (int id = 0; id < intersections.size(); id++) {
// iterate over consecutive pairs of Iterations
Intersection i1 = intersections.get(id);
Intersection i2 = intersections.get((id + 1) % intersections.size()); // wrap around to 0
if (i1.X() != null && i2.X() != null) {
// sending to candidates anyway
float dSum = i1.minD + i2.minD; // sum of distances from selected rays (s-pt) to their
// intersection point
}
}

}
}

// sort rays by distance from their s-pt to their intersection pt
ArrayList<Candidate> srays = new ArrayList<>(candidates);
Collections.sort(srays);

ArrayList<Candidate> selectedRays = new ArrayList<>();
for (Candidate r : srays) {
if (r.minD == srays.get(0).minD) {
}
}

if (!selectedRays.isEmpty()) {
int offset = 0;

for (Candidate ray : selectedRays) {
int id = ray.id;
PVector X = ray.intersectionPoint;
Ray r1 = rays.get(id - offset);
Ray r2 = rays.get((id + 1 - offset) % lr);

// stores bones (segments from parent rays to intersection point 'X')
Stream.of(r1, r2).forEach(r -> {
if (!vertices.contains(r.startPoint)) {
float d1 = X.dist(r1.prevEdge.p2);
float d2 = X.dist(r2.nextEdge.p1);
if ((d1 + d2) / 2 > tol) {
} else {
}
} else {
}
});

// compute direction of the new ray
PVector b = bisector(r1, r2);

// check if ray points in the right direction
int si = vertices.indexOf(r2.nextEdge.p1); // #start-pt index
int ei = vertices.indexOf(r1.prevEdge.p2); // #end-pt index

ArrayList<PVector> slice = new ArrayList<>();
slice.add(X); // moved from '[X] + slice' to here
if (ei > si) {
for (int i = si; i < ei; i++) {
}
} else {
for (int i = si; i < vertices.size(); i++) {
}
for (int i = 0; i < ei; i++) {
}
}

ept.mult(-1);
}

// delete parents rays from array list and insert their 'child' instead
rays.set(id - offset, new Ray(X, ept, r1.prevEdge, r2.nextEdge));
rays.remove((id + 1 - offset) % lr);
offset++;
lr = rays.size();

}
} else {
System.err.println(
"Error: no rays have been found for reduction. A shared intersection point is probably missing.");
break;
}
}
bones.add(new Bone(rays.get(0).startPoint, rays.get(1).startPoint)); // connect start-points of the last 2 rays
}

/**
* Computes the bisector of a corner between edge p1-p2 and edge p3-p4.
*/
private static PVector bisector(PVector p1, PVector p2, PVector p3, PVector p4) {

PVector dir1 = PVector.sub(p2, p1).normalize(); // direction between 1st point and 2nd point of edge 1
PVector dir2 = PVector.sub(p4, p3).normalize(); // direction between 1st point and 2nd point of edge 2

return new PVector(dsum.y, -dsum.x).normalize();
}

private static PVector bisector(Ray r1, Ray r2) {
return bisector(r1.prevEdge.p1, r1.prevEdge.p2, r2.nextEdge.p1, r2.nextEdge.p2);
}

/**
* Checks if 2 lines are intersecting. Optional: returns location of
* intersection point.
*/
private static PVector intersect(PVector p1, PVector p2, PVector p3, PVector p4) {

float 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));
float 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));

float secX, secY;
if (uA >= 0 && uA <= 1 && uB >= 0 && uB <= 1) {

secX = p1.x + (uA * (p2.x - p1.x));
secY = p1.y + (uA * (p2.y - p1.y));
return new PVector(round(secX), round(secY));
}
return null; // no intersection
}

private static PVector intersect(Ray r1, Ray r2) {
return intersect(r1.startPoint, r1.endPoint, r2.startPoint, r2.endPoint);
}

private static boolean contains(ArrayList<PVector> verts, PVector pt) {

int N = verts.size();
boolean isInside = false;

for (int i = 0; i < N; i++) {
PVector v1 = verts.get((i + 1) % N);
PVector v2 = verts.get(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 = !isInside;
}
}
}

return isInside;
}

/**
* Simple container.
*
* Rays originate from the corners of the polygon, moving inwards.
*/
public class Ray {

public PVector startPoint; // the vertex ray originated from
public PVector endPoint; // live end point of ray
public Edge prevEdge; // edge from previous vertex to startpoint of this ray
public Edge nextEdge; // from to next vertex

public Ray(PVector startPoint, PVector endPoint, Edge e1, Edge e2) {
this.startPoint = startPoint;
this.endPoint = endPoint;
this.prevEdge = e1;
this.nextEdge = e2;
}
}

/**
* Simple container
*/
public class Edge {

public PVector p1; // from
public PVector p2; // to

public Edge(PVector p1, PVector p2) {
this.p1 = p1;
this.p2 = p2;
}
}

/**
* Simple container
*/
class Intersection {

float minD;
PVector intersectionPoint;

public Intersection(PVector intersectionPoint, float minD) {
this.minD = minD;
this.intersectionPoint = intersectionPoint;
}

/**
* Where the rays cross
*/
PVector X() {
return intersectionPoint;
}
}

/**
* Simple container
*/
class Candidate implements Comparable<Candidate> {

int id;
PVector intersectionPoint;
float minD;

public Candidate(int id, PVector intersectionPoint, float minD) {
this.id = id;
this.intersectionPoint = intersectionPoint;
this.minD = minD;
}

@Override
public int compareTo(Candidate o) {
if (minD > o.minD) { // this bigger
return 1;
}
if (minD < o.minD) { // this smaller
return -1;
}
return 0; // equal
}

}

public class Bone {

public PVector sp1; // startPoint of ray 1
public PVector sp2; // startPoint of ray 2

public Bone(PVector sp1, PVector sp2) {
this.sp1 = sp1;
this.sp2 = sp2;
}
}
}
``````
1 Like

can’t you just do this?

Code
``````ArrayList<rct> rcts = new ArrayList<rct>();
float sx, sy, w, h;
void setup() {
size(600,600);
}
void draw() {
background(0);
for(int i = 0; i < rcts.size(); i++) rcts.get(i).display();
for(int i = 0; i < rcts.size(); i++) rcts.get(i).displayLine();
fill(255,50);
//if(mousePressed) rect(sx,sy,abs(sx-mouseX),abs(sy-mouseY));
}
void mousePressed() {
sx = mouseX;
sy = mouseY;
}
void mouseReleased() {
float ex = mouseX, ey = mouseY;
if(ex < sx) {
float temp = sx;
sx = ex;
ex = temp;
}
if(ey < sy) {
float temp = sy;
sy = ey;
ey = temp;
}
w = abs(sx-ex);
h = abs(sy-ey);
}

class rct {
float x,y,w,h;
boolean verticle = false;
rct(float x_, float y_, float w_, float h_) {
x = x_;
y = y_;
w = w_;
h = h_;
verticle = h > w;
}
void display() {
noStroke();
fill(255);
rect(x,y,w,h);
}
void displayLine() {
float a = 0;
if(!verticle) {
a = h/2;
} else {
a = w/2;
}
stroke(255,0,0);
line(x+a,y+a,x+w-a,y+h-a);
}
}
``````

Now do it like that for the airplane .

2 Likes

@micycle – Thank you very much, that is so kind of you to have ported the whole sketch to Java. Much respect for taking the time to get your head around the algorithm and make the conversion from Python!

I should stress however that the code I wrote was only meant to serve as an example to illustrate the topological differences between a Voronoi diagram-based medial axis and a Straight Skeleton-based one. It does not follow any of the logic of the conventional algorithms dedicated to 2D skeletonization. It only handles very basic shapes, is prone to errors and barely works overall.

Good news is that there are now 2 reliable libraries for Processing:

• Trace Skeleton (already suggested in this other thread). Works on pixels (still images, webcam, etc) and outputs a medial axis. Never tried it but looks really powerful.

• CampSkeleton, now available inside Hemesh as an external package. Works on polygons (without holes only) and outputs a straight skeleton.

For retrieving the medial axis of a polygon, one can combine the latter (Hemesh) with a skeleton pruning method like the “Discrete Curve Evolution” algorithm. A naive implementation based on this paper (Python mode):

``````from collections import defaultdict

verts = [PVector(162, 387),PVector(186, 401),PVector(337, 368),PVector(434, 267), PVector(442, 241),
PVector(416, 250), PVector(402, 237), PVector(429, 200), PVector(417, 144), PVector(377, 139),
PVector(301, 212), PVector(254, 203), PVector(225, 236), PVector(203, 236), PVector(193, 195),
PVector(193, 269), PVector(176, 281), PVector(140, 276), PVector(134, 262), PVector(125, 273),
PVector(144, 315)][::-1]

temp = verts[:] # temporatory store vertices

G = defaultdict(set) # graph (store connection info of skeleton graph)
P = sum(p1.dist(p2) for p1, p2 in zip(verts, verts[1:] + verts[0:])) # perimeter
N = len(verts) # total number of vertices

def setup():
size(600, 500)
background(255)

global N

relevances = [relevance(i) for i in xrange(N)]

#### COMPUTE STRAIGHT SKELETON ####

# Convert vertices to WB_Points (must be in clockwise order)
polyverts = [WB_Point(p.x, p.y) for p in verts]

# Set a very large offset for each vertex (similar to a wavefront so far projected it will output the whole skeleton)
offsets =  * N

# Some HE_Mesh magic trick
pyr = WB_PyramidFactory().setPoints(polyverts)
mesh = HE_Mesh(pyr.createOffset(0, offsets))

# Store Skeleton edges
edges = set()
for i1, i2 in mesh.getEdgesAsInt():
v1 = mesh.getVertex(i1).getPosition()
v2 = mesh.getVertex(i2).getPosition()
if v1 in polyverts and v2 in polyverts:
continue

#### DCE ####

# As long as the perimeter of the Discrete Curve is above 90% of the perimeter of the original perimeter
while sum(p1.dist(p2) for p1, p2 in zip(verts, verts[1:] + verts[0:])) > P * .9:

# vertex index with lowest relevance
idmin = relevances.index(min(relevances))

# remove whole branch starting at 'idmin' vertex
p = verts[idmin]
try:
cur = p
while True:
n = next(iter(G[cur]))

G[n] -= set([cur])
del G[cur]

if len(G[n]) > 1:
break

cur = n

except:
print "no"
pass

# remove corresponding vertex and relevance value
del verts[idmin]
del relevances[idmin]
N -= 1

# compute new relevance value of neighboring vertices
relevances[idmin%N] = relevance(idmin%N)
relevances[(idmin-1)%N] = relevance((idmin-1)%N)

# remove branches from concave vertices (not sure, concave in the original polygon (temp) or concave in the DCE (verts) ?)
for i in xrange(N):
cv = verts[i]
pv = verts[(i-1)%N]
nv = verts[(i+1)%N]
if not isConvex(pv, cv, nv):
try:
cur = cv
while True:
n = next(iter(G[cur]))
G[n] -= set([cur])
del G[cur]

if len(G[n]) > 1:break

cur = n

except: pass

#### RENDER ####

# Shape (polygon)
pushStyle()
strokeWeight(3)
for v1, v2 in zip(polyverts, polyverts[1:] + polyverts[0:]):
line(v1.x, v1.y, v2.x, v2.y)
popStyle()

# Skeleton
for v1, v2 in edges:
line(v1.x, v1.y, v2.x, v2.y)

# Discrete Curve
pushStyle()
strokeWeight(2)
stroke(0, 100, 255)
for v1, v2 in zip(verts, verts[1:] + verts[0:]):
line(v1.x, v1.y, v2.x, v2.y)
popStyle()

# Medial Axis
pushStyle()
strokeWeight(2)
stroke(255, 30, 30)
for p1 in G:
for p2 in G[p1]:
line(p1.x, p1.y, p2.x, p2.y)
popStyle()

noLoop()

def relevance(i):

pv = verts[(i-1)%N]
cv = verts[i]
nv = verts[(i+1)%N]

a = atan2(nv.y-cv.y, nv.x-cv.x) - atan2(cv.y-pv.y, cv.x-pv.x)
l1 = PVector.sub(cv, pv).mag()
l2 = PVector.sub(nv, cv).mag()

return (abs(a)*l1*l2) / (l1+l2) / P

def isConvex(p1, p2, p3):

return (p2.x - p1.x) * (p3.y - p2.y) - (p3.x - p2.x) * (p2.y - p1.y) > 0
``````

2 Likes

I’ve made a sketch which does just this a while back. It was originally a canny attempt so essentially its calculating the skeleton of the edges, but is also able to calculate contours.

Theres currently a severe lack of comments on it and some of the other functions are nothing more than work in progress, but I’m happy to add comments and answer any questions you may have.

As for the algorithm used I just create a sobel of an image then for each black point I travel up down left right and diagonal in every direction until I find a white pixel store the steps required and then place a black pixel where the step count in opposite directions are equal to each other.

I think your code is better than you give it credit for – provided there’s no holes!

Although, yeah… This happens somewhat often…

I’d actually seen camp-skeleton, but I had problems using it, so I looked around and found your Python approach instead, and decided to try to port it – thankfully you included input data, so I could validate and debug it quite easily. It’s nice when you have an example to compare it to!

I’d also come across Stefan Huber’s work work on the matter, but straight-skeletons ain’t an easy problem (maybe that’s why there are almost no solutions out there!).

Two points:

A slightly different topic: would you mind sharing how you tackled generating the mitered offset of a polygon (given its skeleton); particularly about culling the correct parts of offset segments in a concave cell.  I’ve recently done my own medial axis implementation in Processing using JTS. I use the technique of: densifying the polygon’s lines; computing voronoi diagram; culling lines where at least one coordinate is outside the polygon.

Of course, the resulting medial axis is comprised of many small disconnected (but overlapping) line segments, and can have small prongs, which are ripe for pruning. It looks like that paper solves both of these problems. Could you briefly explain (in English) how the axis-line merging and pruning works in it?

I must admit, I laughed. (that output is horrific)

Yes! I remember of this work. It’s probably one of these pages everyone stumble upon when doing some research on straight skeletons.

I don’t remember posting about mitered offsets here, are you referring to this animation I made to illustrate a past SO question of mine ? If so, I just:

• converted each polygonal edge to a ray
• offsetted that parallel ray inward
• looked for intersections with the boundaries of the corresponding skeleton face

Example with Hemesh:

``````add_library('hemesh')

verts = [PVector(144, 315), PVector(125, 273), PVector(134, 262), PVector(140, 276), PVector(176, 281), #05
PVector(193, 269), PVector(193, 195), PVector(203, 236), PVector(225, 236), PVector(254, 203), #10
PVector(301, 212), PVector(377, 139), PVector(417, 144), PVector(429, 200), PVector(402, 237), #15
PVector(416, 250), PVector(442, 241), PVector(434, 267), PVector(337, 368), PVector(186, 401), #20
PVector(162, 387)]

def setup():
size(600, 500)
background('#FFFFFF')

#### COMPUTE STRAIGHT SKELETON ####

#Convert vertices to WB_Points (must be in clockwise order)
polyverts = [WB_Point(p.x, p.y) for p in verts]

#Set a very large offset for each vertex (similar to a wavefront so far projected it will output the whole skeleton)
offsets =  * len(verts)

#Some HE_Mesh magic trick
pyr = WB_PyramidFactory().setPoints(polyverts)
mesh = HE_Mesh(pyr.createOffset(0, offsets))

#Store Skeleton edges (just for rendering)
edges = set()
for i1, i2 in mesh.getEdgesAsInt():
v1 = mesh.getVertex(i1).getPosition()
v2 = mesh.getVertex(i2).getPosition()

#make sure to cull edges belonging to the polygonal shape (select skeleton only)
if v1 in polyverts and v2 in polyverts:
continue

#### COMPUTE MITERED OFFSETS ####

contours = []

#Retrieve inner faces of the skeleton
for i, face_id in enumerate(mesh.getFacesAsInt()):
start_index = None
face_verts = []

#Get each vertex of visited face & store them
for idx, n in enumerate(face_id):
v = mesh.getVertex(n).getPosition()
v = PVector(v.x, v.y) # convert to PVector
face_verts.append(v)

if v == verts[i]:
start_index = idx

#Re-order vertices & edges so they all start from the original polygon edges
face_verts = face_verts[start_index:] + face_verts[:start_index]
face_edges = zip(face_verts, face_verts[1:] + face_verts[:1])

#Start-point & end-point of the 1st edge
spt, ept = face_edges

#Convert 1st edge to ray & offset multiple times
for gap in xrange(0, 100, 5):

dir = PVector.sub(ept, spt).normalize()
theta = dir.heading() + HALF_PI #angle perpendicular to the edge

#Offset ray parallely
vec = PVector.fromAngle(theta).setMag(gap)

#Look for intersections with the boundaries of current "inner" face
intersections = []
for p3, p4 in face_edges[1:]:
X = intersect(p1, p2, p3, p4)
if X:
intersections.append(X)

#If any, store them grouped by 2
if intersections:
it = iter(intersections)
contours.extend(zip(it, it))

#### RENDER ####

#Polygonal shape
pushStyle()
strokeWeight(3)
for v1, v2 in zip(verts, verts[1:] + verts[0:]):
line(v1.x, v1.y, v2.x, v2.y)
popStyle()

#Straight-Skeleton
for v1, v2 in edges:
line(v1.x, v1.y, v2.x, v2.y)

#Offsets
for v1, v2 in contours:
line(v1.x, v1.y, v2.x, v2.y)

noLoop()

def intersect(p1, p2, p3, p4):

'''Returns a PVector.
Checks if 2 lines are intersecting. If so, 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)) + .01)
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)) + .01)

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(secX, secY)
`````` Indeed, that would be the Voronoi diagram based technique I mentioned and tested above. It sure works well for certain cases!

Although I vaguely remember it was quite easy to implement I unfortunately can’t recall the details of this method. Also I must apologize as I’m now even starting to doubt that the paper I linked previously is the appropriate one. From my recollection the creator of this algorithm is a certain Mr. Latecki.
Really sorry, it’s all I got now. I’m sure you can find plenty of better implementations on GitHub for insights.

4 Likes

Hi there,
After scouring the web for this topic I was excited to find this post in the forum and amazed at the level of detail and help, particularly yourself (@micycle ) & @solub… so i feel a little bad for asking for further help. I tried using the JAVA provided but get some errors (and tried porting the python to use in a processing sketch… umm let’s just say that that’s a day of my life i’m never getting back) I’m a perpetual programming/processing newb but generally find my way using the processing console & simple language/guidance in the p5 reference.

Anyway, my question: it possible to get this running in the processing console so it can be used with other processing sketches? For example I get a token unexpected error at this line (“r->” is very weird to me)

Although I’m not familiar with Streams, Collections & HashSets I feel that if i knew what that line of code was doing i might be able to wrangle it into something closer to basic processing language.
If anyone has an explanation of that line or tip to get this working in processing I’d be super grateful!

``````Stream.of(r1, r2).forEach(r -> {
``````

@ReinoudvanLaar, this was just an easy way to iterate a pair of objects: create a stream consisting of the pair, then iterate over it (the alternative is creating an array and adding each element to it separately).

Use Processing 4 – it should support the lambda (unlike Processing 3) since this is a Java 8+ feature.

Alternatively, a better straight skeleton implementation is wrapped by my PGS library – with that you can create straight skeletons in Processing using a single line of code.

@micycle Thanks for your reply & advice, I’m looking forward to checking out your library I’ve revisited this topic recently after unearthing a straight-skeleton implementation that is ~50x faster than camp-skeleton (however, it can produce malformed outputs). It will be integrated into the next version of PGS.

Here are some illustrations of what can be done with PGS (namely the `PGS_Contour.straightSkeleton()`, `PGS_Coloring.colorMesh()` and `PGS_Contour.offsetCurvesInward()` methods).

1 Like