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

@micycle Thanks for your reply & advice, I’m looking forward to checking out your library :slight_smile:

I’ve revisited this topic recently after unearthing a straight-skeleton implementation that is ~50x faster than camp-skeleton (however, it can produce malformed outputs). It will be integrated into the next version of PGS.

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

2 Likes

I Converted the Java Medial Axis from Straight Skeleton to Unity3d C#

using System.Collections.Generic;
using System.Linq;
using UnityEngine;

/**
 * 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
{
    private float tol;
    public List<Vector2> vertices;
    public List<Ray> rays;
    public List<Bone> bones;
    public List<Bone> branches;


    public SolubSkeleton(float tol, List<Vector2> vertices)
    {
        this.tol = tol;
        this.vertices = vertices;
        rays = new List<Ray>();
        bones = new List<Bone>();
        branches = new List<Bone>();
        Run();
    }
    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.Count;

        for (int i = 0; i < lv; i++)
        {
            Vector2 pv = vertices[(i - 1 + lv) % lv];  // previous vertex
            Vector2 cv = vertices[i]; // current vertex
            Vector2 nv = vertices[(i + 1) % lv]; // next vertex

            Vector2 b = Bisector(pv, cv, cv, nv);
            Vector2 ept = cv + b.normalized * 100000; // 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); // start-point, end-point, previous edge, next edge
        }
    }

    private void Reduce()
    {
        /*
     * Check for ray-ray collision
    */
        while (rays.Count > 2)
        {
            List<Intersection> intersections = new List<Intersection>();
            int lr = rays.Count();

            for (int i = 0; i < lr; i++)
            {
                Ray pRay = rays[(i - 1 + lr) % lr]; // previous ray
                Ray cr = rays[i]; // current ray
                Ray nRay = rays[(i + 1) % lr]; // next ray

                float minD = 100000; // min-distance
                Vector2 a = Vector2.zero; // intersection point (optional)

                // Check current ray for intersection with prev and next
                Vector2 x = Intersect(cr, pRay); // check with prev
                if (x != Vector2.zero) // rays intersected
                {
                    float d = (x - (Vector2)cr.startPoint).magnitude;
                    if (d < minD)
                    {
                        minD = d;
                        a = x;
                    }
                }

                x = Intersect(cr, nRay); // check with next
                if (x != Vector2.zero) // rays intersected
                {
                    float d = (x - (Vector2)cr.startPoint).magnitude;
                    if (d < minD)
                    {
                        minD = d;
                        a = x;
                    }
                }

                intersections.Add(new Intersection(a, minD));
            }

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

            for (int i = 0; i < intersections.Count; i++)
            {
                Intersection i1 = intersections[i];
                Intersection i2 = intersections[(i + 1) % intersections.Count]; // wrap around to 0 for last pair

                if (i1.X() != Vector2.zero && i2.X() != Vector2.zero) // if X not null and equivalent
                {
                    if (i1.X() == 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.Count == 0)
            {
                // making sure that NO consecutive rays intersect each other before either
                // intersects its other neighbor
                if (new HashSet<Intersection>(intersections).Count == intersections.Count)
                {
                    // similar loop to above
                    for (int id = 0; id < intersections.Count; id++)
                    {
                        // iterate over consecutive pairs of Iterations
                        Intersection i1 = intersections[id];
                        Intersection i2 = intersections[(id + 1) % intersections.Count]; // wrap around to 0
                        if (i1.X() != Vector2.zero && i2.X() != Vector2.zero)
                        {
                            // 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
            List<Candidate> srays = new List<Candidate>(candidates);
            srays.Sort();

            List<Candidate> selectedRays = new List<Candidate>();
            foreach (Candidate r in srays)
            {
                if (r.minD == srays[0].minD)
                {
                    selectedRays.Add(r);
                }
            }

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

                foreach (var ray in selectedRays)
                {
                    int id = ray.id;
                    Vector2 X = ray.intersectionPoint;
                    Ray r1 = rays[id - offset];
                    Ray r2 = rays[(id + 1 - offset) % lr];

                    // stores bones (segments from parent rays to intersection point 'X')
                    foreach (Ray r in new Ray[] { r1, r2 })
                    {
                        if (!vertices.Contains(r.startPoint))
                        {
                            float d1 = Vector2.Distance(X, r1.prevEdge.p2);
                            float d2 = Vector2.Distance(X, 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
                    Vector2 b = Bisector(r1, r2);
                    Vector2 ept = X + b.normalized * 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

                    List<Vector2> slice = new List<Vector2>();
                    slice.Add(X); // moved from '[X] + slice' to here
                    if (ei > si)
                    {
                        for (int i = si; i < ei; i++)
                        {
                            slice.Add(vertices[i]); // vertices[si:ei]
                        }
                    }
                    else
                    {
                        for (int i = si; i < vertices.Count; i++)
                        {
                            slice.Add(vertices[i]); // vertices[si:]
                        }
                        for (int i = 0; i < ei; i++)
                        {
                            slice.Add(vertices[i]); // self.vertices[:ei]
                        }
                    }

                    if (!Contains(slice, X + b.normalized))
                    {
                        ept *= -1;
                    }
                    // delete parents rays from array list and insert their 'child' instead
                    rays[id - offset] = new Ray(X, ept, r1.prevEdge, r2.nextEdge);
                    rays.RemoveAt((id + 1 - offset) % lr);
                    offset++;
                    lr = rays.Count;

                }

            }else {
             Debug.LogError(
                    "Error: no rays have been found for reduction. A shared intersection point is probably missing.");
                break;
            }
        }
        bones.Add(new Bone(rays[0].startPoint, rays[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 Vector2 Bisector(Vector2 p1, Vector2 p2, Vector2 p3, Vector2 p4)
    {
        Vector2 dir1 = (p2 - p1).normalized; // direction between 1st point and 2nd point of edge 1
        Vector2 dir2 = (p4 - p3).normalized; // direction between 1st point and 2nd point of edge 2

        Vector2 dsum = dir1 + dir2;

        return new Vector2(dsum.y, -dsum.x).normalized;
    }

    private static Vector2 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 Vector2 Intersect(Vector2 p1, Vector2 p2, Vector2 p3, Vector2 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 Vector2(Mathf.Round(secX), Mathf.Round(secY));
        }

        return default;
    }

    private static Vector2 Intersect(Ray r1, Ray r2)
    {
        return Intersect(r1.startPoint, r1.endPoint, r2.startPoint, r2.endPoint);
    }

    private static bool Contains(List<Vector2> verts, Vector2 pt)
    {
        int N = verts.Count;
        bool isInside = false;

        for (int i = 0; i < N; i++)
        {
            Vector2 v1 = verts[(i + 1) % N];
            Vector2 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 = !isInside;
                }
            }
        }

        return isInside;
    }

    /**
	 * Simple container.
	 *
	 * Rays originate from the corners of the polygon, moving inwards.
	 */
    public class Ray
    {
        public Vector2 startPoint; // the vertex ray originated from
        public Vector2 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(Vector2 startPoint, Vector2 endPoint, Edge e1, Edge e2)
        {
            this.startPoint = startPoint;
            this.endPoint = endPoint;
            this.prevEdge = e1;
            this.nextEdge = e2;
        }
    }

    /**
	 * Simple container
	 */
    public class Edge
    {
        public Vector2 p1; // from
        public Vector2 p2; // to

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

    /**
	 * Simple container
	 */
    class Intersection
    {
        public float minD;
        public Vector2 intersectionPoint;

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

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

    public class Candidate : System.IComparable<Candidate>
    {
        public int id;
        public Vector2 intersectionPoint;
        public float minD;

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

        public int CompareTo(Candidate other)
        {
            if (minD > other.minD)
            {
                return 1;
            }
            else if (minD < other.minD)
            {
                return -1;
            }
            else
            {
                return 0;
            }
        }
    }

    public class Bone
    {
        public Vector2 sp1; // startPoint of ray 1
        public Vector2 sp2; // startPoint of ray 2

        public Bone(Vector2 sp1, Vector2 sp2)
        {
            this.sp1 = sp1;
            this.sp2 = sp2;
        }
    }
}

Sample

using System.Collections.Generic;
using UnityEngine;

namespace DefaultNamespace
{   [ExecuteInEditMode]
    public class Skeleton : MonoBehaviour
    {
        public List<Vector2> verts = new List<Vector2>
        {
            new Vector2(901.12f, 300),
            new Vector2(853.48f, 279.76f),
            new Vector2(768.04f, 270.04f),
            new Vector2(701.56f, 261.04f),
            new Vector2(684.4f, 253.72f),
            new Vector2(684.4f, 250.48f),
            new Vector2(701.32f, 234.76f),
            new Vector2(701.2f, 231.28f),
            new Vector2(685.84f, 223.48f),
            new Vector2(631.12f, 220.6f),
            new Vector2(570.52f, 218.32f),
            new Vector2(534.88f, 200.68f),
            new Vector2(434.2f, 50.56f),
            new Vector2(362.68f, 50.56f),
            new Vector2(320.32f, 218.44f),
            new Vector2(306.16f, 213.52f),
            new Vector2(258.28f, 140.68f),
            new Vector2(191.44f, 140.68f),
            new Vector2(164.8f, 253.12f),
            new Vector2(233.68f, 262.96f),
            new Vector2(245.92f, 265.12f),
            new Vector2(252.64f, 268.72f),
            new Vector2(278.2f, 268.12f),
            new Vector2(265.0f, 276.04f),
            new Vector2(244.24f, 284.68f),
            new Vector2(244.24f, 315.32f),
            new Vector2(265.0f, 323.96f),
            new Vector2(278.2f, 331.88f),
            new Vector2(252.64f, 331.28f),
            new Vector2(245.92f, 334.88f),
            new Vector2(233.68f, 337.04f),
            new Vector2(164.8f, 346.88f),
            new Vector2(191.44f, 459.32f),
            new Vector2(258.28f, 459.32f),
            new Vector2(306.16f, 386.47998f),
            new Vector2(320.32f, 381.56f),
            new Vector2(362.68f, 549.44f),
            new Vector2(434.2f, 549.44f),
            new Vector2(534.88f, 399.32f),
            new Vector2(570.52f, 381.68f),
            new Vector2(631.12f, 379.4f),
            new Vector2(685.84f, 376.52002f),
            new Vector2(701.2f, 368.72f),
            new Vector2(701.32f, 365.24f),
            new Vector2(684.4f, 349.52002f),
            new Vector2(684.4f, 346.28f),
            new Vector2(701.56f, 338.96f),
            new Vector2(768.04f, 329.96f),
            new Vector2(853.48f, 320.24f)
        };

        public List<Vector2> vertices;
        public List<Ray> rays;
        public List<SolubSkeleton.Bone> bones;
        public List<SolubSkeleton.Bone> branches;

       public float tol = 20;

        void Update()
        {
            SolubSkeleton skeleton = new SolubSkeleton(tol, verts);
            branches = skeleton.branches;
            bones = skeleton.bones;
        }

        private void OnDrawGizmos()
        {
            Gizmos.color = Color.yellow;
            
            for (int i = 0; i < verts.Count; i++)
            {
                Vector2 p1 = verts[i];
                Vector2 p2 = verts[(i + 1) % verts.Count];
                Gizmos.DrawLine(new Vector3(p1.x, 0, p1.y),new Vector3(p2.x, 0, p2.y) );
            }

            Gizmos.color = Color.red;

            foreach (Vector2 p in verts)
            {
                Gizmos.DrawSphere(new Vector3(p.x, 0, p.y), 2f); // Adjust the radius as needed
            }


            // Draw Branches
            Gizmos.color = new Color32(214, 235, 222, 255);
            foreach (SolubSkeleton.Bone b in branches)
            {
                Gizmos.DrawLine(new Vector3(b.sp1.x, 0, b.sp1.y),new Vector3(b.sp2.x, 0, b.sp2.y) ); 
            }

            // Draw Medial-Axis
            Gizmos.color = new Color32(103, 255, 101, 255);
            foreach (SolubSkeleton.Bone b in bones)
            {
                Gizmos.DrawLine(new Vector3(b.sp1.x, 0, b.sp1.y),new Vector3(b.sp2.x, 0, b.sp2.y) );
            }
        }
    }
}