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

Hello,

for a project I want to draw a “center line” into a shape/polygon. The picture below shows the letter “H” outlined (this is the data I have) and then the center line in red (this is the data I am looking for). I am looking for a way to calculate the red line.

hskelet

My code stores the coordinates of the shape (black line) in a PShape. So what I am looking for is a function or library which takes the PShape and returns the red line as PShape or even as a list of coordinates.

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

Thanks
Mike

1 Like

Hi @Maxs,

Quick questions:

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

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

3 Likes

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

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

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

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

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[0].x, seg[0].y, seg[1].x, seg[1].y)
    popStyle()
    
    
    noLoop()


class MedialAxis(object):
    def __init__(self, verts):
        self.vertices = verts
        self.segments = []
        self.edges = []
        self.run()
    
        
    def run(self):
        self.initialize
        self.cull
        
    
        
    @property
    def initialize(self):
        
        '''Computes a Voronoi diagram.
           Stores the vertices of each region in a separate sub-list.'''
        
        #instantiate Toxiclibs Voronoi() class
        voronoi = Voronoi() 
        
        #convert vertices to Vec2D format
        for v in self.vertices:
            voronoi.addPoint(Vec2D(v.x, v.y))
        
        #store unique edges of the Voronoi diagram
        for i, poly in enumerate(voronoi.getRegions()):
            
            #convert vertices from Vec2D() to PVector()
            vertices = [PVector(v.x(), v.y()) for v in poly.vertices] 
            
            for v1, v2 in zip(vertices, vertices[1:] + vertices[:1]):
                if (v1, v2) not in self.edges and (v2, v1) not in self.edges:
                    self.edges.append((v1, v2))
     
       
       
    @property    
    def cull(self):
        
        '''Finds medial-axis segments.
           Stores Voronoi cell's edge if it lies entirely inside the main polygon.'''
        
        for p1, p2 in self.edges:
            isValid = False
            if self.contain(verts, p1) and self.contain(verts, p2):
                isValid = True
                for v1, v2 in zip(self.vertices, self.vertices[1:] + self.vertices[:1]):
                    if self.intersect(p1, p2, v1, v2):
                        isValid = False
                        break
            
            if isValid:
                self.segments.append((p1, p2))
            


    def contain(self, verts, pt):
        
        '''Returns a boolean.
           Determines whether a point lies inside a shape/polygon or not.'''
        
        N = len(verts)
        isInside = 0
        
        for i in xrange(N):
            v1 = verts[(i+1)%N]
            v2 = verts[i]
            
            if (v2.y > pt.y) != (v1.y > pt.y):
                if pt.x < (v1.x - v2.x) * (pt.y - v2.y) / (v1.y - v2.y) + v2.x:
                    isInside = not isInside
                    
        return isInside
        
        
        
    def intersect(self, p1, p2, p3, p4):
            
        '''Returns a Boolean.
           Checks if 2 lines are intersecting. Option: returns location of intersection point'''
        
        uA = ((p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y+1)) 
        uB = ((p2.x-p1.x)*(p1.y-p3.y) - (p2.y-p1.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y+1))
        
        if uA >= 0 and uA <= 1 and uB >= 0 and uB <= 1:
            
            return 1

The second approach, and undoubtedly the most the complex one, consists in computing the Straight Skeleton of a polygon and remove the “branches” that connect the inner bones to the vertices of that polygon.

Algorithm:

  • calculate the bisecting vector at each vertex of a chosen polygon. Make it a ray.
  • find the intersection points with the previous and the next neighboring ray.
  • find all pairs of rays that intersect sooner than either intersects its other neighbor.
  • sort them by distance to their intersection point.
  • take the pair with the shortest distance and replace it with a single new ray originating from their intersection point. That child ray must be aiming towards the average direction of its neighboring polygon edges.
  • repeat from step 2 until there are only 2 rays left.
  • connect the originating points of the last 2 rays to complete the skeleton

Below a naive implementation in Processing Python:

Medial Axis from Straight Skeleton
verts = [PVector(901.12, 300),PVector(853.48, 279.76),PVector(768.04, 270.04),PVector(701.56, 261.04),
         PVector(684.4, 253.72),PVector(684.4, 250.48),PVector(701.32, 234.76),PVector(701.2, 231.28),
         PVector(685.84, 223.48),PVector(631.12, 220.6),PVector(570.52, 218.32),PVector(534.88, 200.68),
         PVector(434.2, 50.56),PVector(362.68, 50.56),PVector(320.32, 218.44),PVector(306.16, 213.52),
         PVector(258.28, 140.68),PVector(191.44, 140.68),PVector(164.8, 253.12),PVector(233.68, 262.96),
         PVector(245.92, 265.12),PVector(252.64, 268.72),PVector(278.2, 268.12),PVector(265.0, 276.04),
         PVector(244.24, 284.68), PVector(244.24, 315.32), PVector(265.0, 323.96),PVector(278.2, 331.88),
         PVector(252.64, 331.28),PVector(245.92, 334.88),PVector(233.68, 337.04),PVector(164.8, 346.88),
         PVector(191.44, 459.32),PVector(258.28, 459.32),PVector(306.16, 386.47998),PVector(320.32, 381.56),
         PVector(362.68, 549.44),PVector(434.2, 549.44),PVector(534.88, 399.32),PVector(570.52, 381.68),
         PVector(631.12, 379.4),PVector(685.84, 376.52002),PVector(701.2, 368.72),PVector(701.32, 365.24),
         PVector(684.4, 349.52002),PVector(684.4, 346.28),PVector(701.56, 338.96),PVector(768.04, 329.96),
         PVector(853.48, 320.24)]



def setup():
    size(1000, 600, P2D)
    background('#FFFFFF')
    smooth(8)
    
    skeleton = Skeleton(verts, 20) #(vertices, threshold)
    branches = skeleton.branches
    bones = skeleton.bones


    #Draw polygon's silhouette
    pushStyle()
    stroke(50, 100, 255)
    strokeWeight(3)
    for p1, p2 in zip(verts, verts[1:] + verts[:1]):
        line(p1.x, p1.y, p2.x, p2.y)
    popStyle()
    
    
    #Draw vertices
    pushStyle()
    stroke(30)
    strokeWeight(6)
    for i, p in enumerate(verts):
        point(p.x, p.y)
    popStyle()
    
    
    #Draw Branches
    pushStyle()
    strokeWeight(.3)
    stroke(40, 235, 180)
    for i, b in enumerate(branches):
        line(b[0].x , b[0].y, b[1].x , b[1].y)
    popStyle()
    
    
    #Draw Medial-Axis
    pushStyle()
    stroke(255, 70, 122)
    strokeWeight(1.8)
    for b in bones:
        line(b[0].x , b[0].y, b[1].x , b[1].y)
    popStyle()
    
    
    noLoop()


             
                    

class Skeleton(object):
    def __init__(self, verts, tol = 0):
        self.vertices = verts
        self.branches = []
        self.bones = []
        self.rays = []
        self.tol = tol
        self.run()
        
        
    def run(self):
        
        '''Runs the algorithm.'''
        
        self.initialize
        self.reduce
                
        
    @property
    def initialize(self):
        
        '''Starts by computing the rays originating from the corners of the polygon. 
           Each ray object also stores the coordinates of its neighboring polygon edges (previous and next).'''
        
        lv = len(self.vertices)
        
        for i in xrange(lv):
            pv = self.vertices[(i-1)%lv]  #previous vertex
            cv = self.vertices[i] #current vertex
            nv = self.vertices[(i+1)%lv] #next vertex
                    
            b = self.bisector(pv, cv, cv, nv)
            ept = b.setMag(100000).add(cv) #end-point of the ray (far from start-point)
            self.rays.append((cv, ept, (pv, cv), (cv, nv))) #start-point, end-point, previous edge, next edge 
            
       
       
    @property            
    def reduce(self):
        
        '''Finds consecutive rays intersecting each other before either intersects its other neighbor. 
           Replaces these pairs of rays by a single one originating from their intersection point and 
           aiming towards the average direction of its neighboring polygon edges.'''
        
        while len(self.rays) > 2:
        
            intersections = [] # Nearest intersection points from parent
            lr = len(self.rays)
        
            for i in xrange(lr):
                pr = self.rays[(i-1)%lr] #previous ray 
                cr = self.rays[i] #current ray
                nr = self.rays[(i+1)%lr] #next ray
                
                mind, X = 100000, None #min-distance, intersection point
                
                for r in (pr, nr):
                    x = self.intersect(cr[0], cr[1], r[0], r[1])
                    if x:
                        d = x.dist(cr[0])
                        if d < mind:
                            mind = d
                            X = x

                intersections.append((X, mind))
            
            candidates = []
            for id, ((X1, d1), (X2, d2)) in enumerate(zip(intersections, intersections[1:]+intersections[:1])): 
                if (X1 and X2) != None and X1 == X2:
                    dsum = d1 + d2 #sum of distances from selected rays (s-pt) to their intersection point
                    candidates.append((id, X1, dsum))
                    
            if not candidates:
                
                #making sure that NO consecutive rays intersect each other before either intersects its other neighbor
                if len(set(intersections)) == len(intersections):
                    for id, ((X1, d1), (X2, d2)) in enumerate(zip(intersections, intersections[1:]+intersections[:1])):
                        if (X1 and X2) != None:
                            
                            #sending to candidates anyway
                            dsum = d1 + d2 #sum of distances from selected rays (s-pt) to their intersection point
                            candidates.append((id, X1, dsum))
        
            srays = sorted(candidates, key=lambda t: t[2]) #rays sorted by distance from their s-pt to their intersection pt
            selectedrays =  [r for r in srays if r[2] == srays[0][2]] #in cases when several rays are equally distant from their intersection pt
            
            if selectedrays:
                offset = 0
                
                for ray in selectedrays:
                    id, X, _ = ray
                    r1 = self.rays[id-offset]
                    r2 = self.rays[(id+1-offset)%lr]
                    
                    #stores bones (segments from parent rays to intersection point 'X')
                    for r in (r1[0], r2[0]):
                        if r not in self.vertices:
                            d1 = X.dist(r1[2][1])
                            d2 = X.dist(r2[3][0])
                            if (d1 + d2) / 2 > self.tol:
                                self.bones.append((r, X))
                            else:
                                self.branches.append((r, X))
                        else:
                            self.branches.append((r, X))
                        

                    #compute direction of the new ray
                    b = self.bisector(r1[2][0], r1[2][1], r2[3][0], r2[3][1]) #s-pt & e-pt of 1st edge, s-pt & e-pt of 2nd edge
                    ept = X.copy().add(b.setMag(100000))
                    
                    #check if ray points in the right direction
                    si = self.vertices.index(r2[3][0]) #s-pt index
                    ei = self.vertices.index(r1[2][1]) #e-pt index
                    slice = self.vertices[si:ei] if ei > si else self.vertices[si:] + self.vertices[:ei]
                    
                    if not self.contain([X] + slice, X.copy().add(b.setMag(1))): 
                        ept.mult(-1)
    
                    #delete parents rays from array list and insert their 'child' instead  
                    self.rays[id-offset] = (X, ept, r1[2], r2[3])
                    del self.rays[(id+1-offset)%lr]
                    offset += 1
                    lr = len(self.rays)
                
            else: 
                print "Error: no rays have been found for reduction. A shared intersection point is probably missing."
                break
                   
                                
        #connect start-points of the last 2 rays    
        self.bones.append((self.rays[0][0], self.rays[1][0]))
                        
         
         
    def bisector(self, p1, p2, p3, p4):
        
        '''Returns a PVector.
           Computes the bisector of a corner between edge p1-p2 and edge p3-p4.'''
        
        dir1 = (p2 - p1).normalize() #direction between 1st point and 2nd point of edge 1
        dir2 = (p4 - p3).normalize() #direction between 1st point and 2nd point of edge 2

        dsum = dir1 + dir2

        return PVector(dsum.y, -dsum.x).normalize() 
        
        
        
    def intersect(self, p1, p2, p3, p4):
        
        '''Returns a Boolean.
           Checks if 2 lines are intersecting. Option: returns location of intersection point'''
        
        uA = ((p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y)) 
        uB = ((p2.x-p1.x)*(p1.y-p3.y) - (p2.y-p1.y)*(p1.x-p3.x)) / ((p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y))
        
        if uA >= 0 and uA <= 1 and uB >= 0 and uB <= 1:
                                            
            secX = p1.x + (uA * (p2.x-p1.x))
            secY = p1.y + (uA * (p2.y-p1.y))
            
            return PVector(round(secX), round(secY))
        
        
        
    def contain(self, verts, pt):
    
        '''Returns a boolean.
           Determine whether a point lies inside a shape/polygon or not.'''
        
        N = len(verts)
        isInside = 0
        
        for i in xrange(N):
            v1 = verts[(i+1)%N]
            v2 = verts[i]
            
            if (v2.y > pt.y) != (v1.y > pt.y):
                if pt.x < (v1.x - v2.x) * (pt.y - v2.y) / (v1.y - v2.y) + v2.x:
                    isInside = not isInside
                    
        return isInside

Please note that this version is made simple for clarity and doesn’t address all the special cases (polygon with holes, with a large number of nearly-parallel edges, those that exhibit symmetrical conditions…) that a more robust library would probably handle with ease.
Yet this is probably one of the few implementations that makes possible to select the bones above a chosen threshold distance ( tol in the code) from the polygon’s vertices. These inner bones are the very segments that form the medial axis.

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

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

2 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

1 Like

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

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

Now do it like that for the airplane :wink:.

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

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.

image image


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 = [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

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.

2 Likes