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

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

}

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

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

// 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)
{
}
}

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)
{
}
else
{
}
}
else
{
}
}

// 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++)
{
}
}
else
{
for (int i = si; i < vertices.Count; i++)
{
}
for (int i = 0; i < ei; i++)
{
}
}

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

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