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

**URL:** <https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208>\
**Category:** Processing\
**Created:** [January 20, 2020, 4:30pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208 "2020-01-20T16:30:24Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Maxs](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/maxs/32/7862_2.png) [@Maxs](https://discourse.processing.org/u/Maxs)\
**Post date:** [January 20, 2020, 4:30pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/1 "2020-01-20T16:30:24Z")

</div>

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.

 ![hskelet](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/a/a298276ae51329b69349d78175cd8f952c92d861.png)

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

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [January 20, 2020, 4:56pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/2 "2020-01-20T16:56:30Z")

</div>

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](https://github.com/Botffy/polyskel), for Java you could try [CampSkeleton](https://github.com/twak/campskeleton) (haven’t tested this one).

---

<div class="post-metadata">

**Author:** ![Maxs](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/maxs/32/7862_2.png) [@Maxs](https://discourse.processing.org/u/Maxs)\
**Post date:** [January 20, 2020, 5:23pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/3 "2020-01-20T17:23:02Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [January 20, 2020, 5:35pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/4 "2020-01-20T17:35:43Z")

</div>

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_.

---

<div class="post-metadata">

**Author:** ![Nathan](https://avatars.discourse-cdn.com/v4/letter/n/c5a1d2/32.png) [@Nathan](https://discourse.processing.org/u/Nathan)\
**Post date:** [January 21, 2020, 12:42am UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/5 "2020-01-21T00:42:29Z")

</div>

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.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/0/099bab58a2140038e0d10decfa3d43262a211106.jpeg)

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).

---

<div class="post-metadata">

**Author:** ![Maxs](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/maxs/32/7862_2.png) [@Maxs](https://discourse.processing.org/u/Maxs)\
**Post date:** [January 21, 2020, 7:23am UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/6 "2020-01-21T07:23:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [January 27, 2020, 8:13pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/7 "2020-01-27T20:13:06Z")

</div>

@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](https://stackoverflow.com/questions/1069523/find-medial-axis-of-a-polygon-using-c-sharp) 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**
>
> ```auto
> 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**
>
> ```auto
> 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.

 ![Annotation 2020-01-27 192412](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/d/de05617f4c03ca08146eeca3b3f18ed81f386bed.jpeg)

---

<div class="post-metadata">

**Author:** ![Maxs](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/maxs/32/7862_2.png) [@Maxs](https://discourse.processing.org/u/Maxs)\
**Post date:** [February 11, 2020, 1:29pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/8 "2020-02-11T13:29:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [February 11, 2020, 7:24pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/9 "2020-02-11T19:24:10Z")

</div>

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

 ![Annotation 2020-02-11 193029](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/b/b027357529d35501a662f043ada411fff1861bc1.png)

> **Test Voronoi**
>
> ```auto
> add_library('geomerative')
> add_library('toxiclibs')
> 
> 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[0].x, seg[0].y, seg[1].x, seg[1].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:
> 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]):
> 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:

- 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.

 ![Annotation 2020-02-11 194115](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/d/dbeb851e322ff25ea5a15a7f3fd64bc1c27d4e9a.png)

---

<div class="post-metadata">

**Author:** ![Maxs](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/maxs/32/7862_2.png) [@Maxs](https://discourse.processing.org/u/Maxs)\
**Post date:** [February 13, 2020, 9:23am UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/10 "2020-02-13T09:23:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![magpie](https://avatars.discourse-cdn.com/v4/letter/m/45deac/32.png) [@magpie](https://discourse.processing.org/u/magpie)\
**Post date:** [May 26, 2020, 6:17pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/11 "2020-05-26T18:17:46Z")

</div>

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](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.  
For example [https://www.evilmadscientist.com/2011/hershey-text-an-inkscape-extension-for-engraving-fonts/](https://www.evilmadscientist.com/2011/hershey-text-an-inkscape-extension-for-engraving-fonts/) seems to have comments about pen plotters vs cnc for the font.  
Inkscape has Hershey Text which can be changed to svg which might do what you need.

---

<div class="post-metadata">

**Author:** ![micycle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/micycle/32/201_2.png) [@micycle](https://discourse.processing.org/u/micycle)\
**Post date:** [February 2, 2021, 3:01pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/12 "2021-02-02T15:01:32Z")

</div>

I just ported the _Medial Axis from Straight Skeleton_ to Java.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/2/2d8eab82d3ae66fe97aa5da8bb69164b85785869.png) ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/9/9ff410c4129632cbe7e7f137aeb033c065694a0a.png)

> **Code**
>
> ```auto
> 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));
> rays.add(r);
> }
> 
> }
> 
> 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;
> }
> }
> 
> intersections.add(new Intersection(X, minD));
> }
> 
> 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;
> candidates.add(new Candidate(i, i1.X(), dSum));
> }
> }
> }
> 
> 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
> candidates.add(new Candidate(id, i1.X(), dSum));
> }
> }
> 
> }
> }
> 
> // 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) {
> selectedRays.add(r);
> }
> }
> 
> 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) {
> bones.add(new Bone(r.startPoint, X));
> } else {
> branches.add(new Bone(r.startPoint, X));
> }
> } else {
> branches.add(new Bone(r.startPoint, X));
> }
> });
> 
> // compute direction of the new ray
> PVector b = bisector(r1, r2);
> PVector ept = X.copy().add(b.setMag(100000));
> 
> // 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++) {
> slice.add(vertices.get(i)); // vertices[si:ei]
> }
> } else {
> for (int i = si; i < vertices.size(); i++) {
> slice.add(vertices.get(i)); // vertices[si:]
> }
> for (int i = 0; i < ei; i++) {
> slice.add(vertices.get(i)); // self.vertices[:ei]
> }
> }
> 
> if (!contains(slice, X.copy().add(b.setMag(1)))) {
> 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
> 
> PVector dsum = PVector.add(dir1, dir2);
> 
> 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;
> }
> }
> }
> 
> ```

---

<div class="post-metadata">

**Author:** ![CodeMasterX](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/codemasterx/32/2942_2.png) [@CodeMasterX](https://discourse.processing.org/u/CodeMasterX)\
**Post date:** [February 2, 2021, 4:53pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/13 "2021-02-02T16:53:41Z")

</div>

can’t you just do this?

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/b/b2913d1935c89220b39c37d354a615ab34715f30.png)

> **Code**
>
> ```auto
> 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);
> rcts.add(new rct(sx,sy,w,h));
> }
> 
> 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);
> }
> }
> 
> ```

---

<div class="post-metadata">

**Author:** ![micycle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/micycle/32/201_2.png) [@micycle](https://discourse.processing.org/u/micycle)\
**Post date:** [February 2, 2021, 5:40pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/14 "2021-02-02T17:40:36Z")

</div>

Now do it like that for the airplane 😉.

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [February 2, 2021, 9:06pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/15 "2021-02-02T21:06:33Z")

</div>

@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](https://skeleton-tracing.netlify.app/) (already suggested in this other [thread](https://skeleton-tracing.netlify.app/)). Works on **pixels** (still images, webcam, etc) and outputs a _medial axis_. Never tried it but looks really powerful.

- [CampSkeleton](https://github.com/twak/campskeleton), now available inside [Hemesh](https://github.com/wblut/HE_Mesh) 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](http://www.solism.ca/bib/solisetal2012skeleton.pdf) (Python mode):

```auto
from collections import defaultdict
add_library('hemesh')

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 = [10000] * 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
        G[PVector(v1.x, v1.y)].add(PVector(v2.x, v2.y))
        G[PVector(v2.x, v2.y)].add(PVector(v1.x, v1.y))
        edges.add((v1, v2))
        
        
    #### 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

```

 ![Discrete Curve Evolution](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/6/69ad0f52f5d79cf34bf4133ef9b3ec15e5d75f8e.png)

---

<div class="post-metadata">

**Author:** ![paulgoux](https://avatars.discourse-cdn.com/v4/letter/p/b9bd4f/32.png) [@paulgoux](https://discourse.processing.org/u/paulgoux)\
**Post date:** [February 2, 2021, 9:59pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/16 "2021-02-02T21:59:36Z")

</div>

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.

> [@OpenCV update, Heaviside Step Function, deriving parametric functions of shapes, Sobel filter](https://discourse.processing.org/t/opencv-update-heaviside-step-function-deriving-parametric-functions-of-shapes-sobel-filter/21931/10):
>
> They seem to use a wide range to techniques depending on the image type and the maths is hidden inside the Wolfram language. The HSF is used to limit the parts of individual curves to include in the calculation so is useful if you are trying to describe the shape using a single pair of parametric curves i.e. x=f(t) and y=g(t) which is what they have in the Einstein face drawing. In your code you are creating the twitter outline using a series of arcs and using the ca array to control the port…

---

<div class="post-metadata">

**Author:** ![micycle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/micycle/32/201_2.png) [@micycle](https://discourse.processing.org/u/micycle)\
**Post date:** [February 2, 2021, 11:13pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/17 "2021-02-02T23:13:11Z")

</div>

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

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/2/2d25d986fb3246ea73ffb22f0b7d22d5a9e8a93f.jpeg)

Although, yeah… This happens somewhat often…

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/c/c90022731f4c2841d2d98865b6eba141d6f86409.jpeg)

* * *

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](https://www.sthu.org/research/straightskeleton/) 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.

![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/f/f3e58420f2482231b6aea2493ab2a636fa32677f.png) ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/8/85426b4fbf5952dada842b03fec6b170bad1505d.png)

* * *

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?

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/3/3894c49e4b0e0cd945350069175b5d87166d3e87.png)

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [February 3, 2021, 2:54pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/18 "2021-02-03T14:54:50Z")

</div>

> [@micycle](#):
>
> Although, yeah… This happens somewhat often…

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

> [@micycle](#):
>
> I’d also come across Stefan Huber’s [work](https://www.sthu.org/research/straightskeleton/) work on the matter

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

> [@micycle](#):
>
> would you mind sharing how you tackled generating the mitered offset of a polygon (given its skeleton)

I don’t remember posting about mitered offsets here, are you referring to [this animation](https://streamable.com/intlj) 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:

```auto
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 = [1000] * 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
        
        edges.add((v1, v2))
        
    
    #### 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[0]

        #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)
            p1 = spt.copy().add(dir.mult(10000.0)).add(vec) #end-pt of ray
            p2 = spt.copy().add(dir.mult(-1.0)).add(vec) #start-pt of ray
                                
            #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)

```

![Mitered Offsets](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/2/26cfe1d463d99eab4e0ad88ea793f90c9219248d.png)

> [@micycle](#):
>
> I’ve recently done my own medial axis implementation […] densifying the polygon’s lines; computing voronoi diagram; culling lines where at least one coordinate is outside the polygon.

Indeed, that would be the Voronoi diagram based technique I mentioned and tested above. It sure works well for certain cases!

> [@micycle](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![ReinoudvanLaar](https://avatars.discourse-cdn.com/v4/letter/r/fbc32d/32.png) [@ReinoudvanLaar](https://discourse.processing.org/u/ReinoudvanLaar)\
**Post date:** [July 30, 2021, 7:05pm UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/19 "2021-07-30T19:05:43Z")

</div>

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)

> [@micycle](#):
>
> `Stream.of(r1, r2).forEach(r -> {`

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!

---

<div class="post-metadata">

**Author:** ![micycle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/micycle/32/201_2.png) [@micycle](https://discourse.processing.org/u/micycle)\
**Post date:** [August 7, 2021, 10:11am UTC](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208/20 "2021-08-07T10:11:06Z")

</div>

```auto
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](https://github.com/micycle1/PGS) – with that you can create straight skeletons in Processing using a single line of code.

[Next page](https://discourse.processing.org/t/straight-skeleton-or-how-to-draw-a-center-line-in-a-polygon-or-shape/17208.md?page=2)