I'm looking for an algorithm that returns all combinations of coordinates ordering that creates polygons with no lines crossing. All coordinates have to be present in the ordering.
For example here are 2 possible ordering of 5 dots:
Any idea how to solve this? Thanks in advance!
I am unaware of an existing algorithm to do this, so I will attempt to find some lower bounds for one, starting with the suggestions in the comments:
n! possible permutations of n points(n-1)!/2, as MBo suggests.a, b, c, d are consecutive in the hull, a, c, b, d will result in a crossing. This greatly reduces the amount of permutations to check. You can build the convex hull in O(n log n) or better, so it is definitely worth it if the hull has many points: a hull with h points allows you to focus on (h-1)!/2 less candidates. On the flip side, if h is a triangle, you gain nothing.Let us build a backtracking algorithm based on the observation that, given a valid closed polygon, and a point p not part of this polygon, all segments e1, e2 such that e1, p and p, e2 cross no polygon segments result in a new valid closed polygon that includes p. We can avoid double-counting (n! vs (n-1)!/2) by choosing a starting triangle.
function find_polygons(points) {
let current = polygon_of(points[0], points[1], points[2]);
let remaining = copy_of(points).removeAll(current);
let answers = new Set();
let goal = points.length();
find_polygons_recursive(goal, current, remaining, answers);
return answers;
}
function find_polygons_recursive(goal, current, remaining, answers) {
if (current.length() == goal) {
answers.add_if_absent(current);
return;
}
for (point p in remaining) {
let without_p = copy_of(remaining).remove(p);
for (segment s in current.segments()) {
if ( ! current.intersects(segment(s1, p)) and
! current.intersects(segment(s2, p))) {
next = copy_of(current).remove(s)
.add(segment(s1,p));
.add(segment(s2,p));
find_polygons_recursive(goal, next, without_p, answers);
}
}
}
}
At the cost of n^4 space and time, you can pre-compute a matrix of n^2xn^2 booleans indicating which segments intersect with which other segments. This should greatly improve the speed of intersection tests.
Here is the C# implementation of @tucxi algorithm described as the accepted answer if someone need it.
Start by calling the function ComputAllPossibleOrdering(coords);
The convex Hull algorithm described here as been used to find the starting points: https://gist.github.com/dLopreiato/7fd142d0b9728518552188794b8a750c (the improvement was huge, depending on the polygon, around 50x faster than starting with a triangle)
Also the intersection algorithm here: https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
The number of possible ordering has been limited to 200. You can easily remove it.
public class UtilsGeometry
{
private List<Tuple<int, int>> crossing_segments = new List<Tuple<int, int>>();
private List<Segment> all_segments = new List<Segment>();
public class Segment
{
public Coords Pt1;
public Coords Pt2;
public int Id;
public Segment(int id, Coords pt1, Coords pt2)
{
Pt1 = pt1;
Pt2 = pt2;
Id = id;
}
}
public class Polygon
{
public List<Segment> Segments = new List<Segment>();
public Polygon(List<Segment> segments)
{
Segments = segments;
}
}
// Given three collinear points p, q, r, the function checks if
// point q lies on line segment 'pr'
private Boolean onSegment(Coords p, Coords q, Coords r)
{
if (q.Lon.Value <= Math.Max(p.Lon.Value, r.Lon.Value) && q.Lon.Value >= Math.Min(p.Lon.Value, r.Lon.Value) &&
q.Lat.Value <= Math.Max(p.Lat.Value, r.Lat.Value) && q.Lat.Value >= Math.Min(p.Lat.Value, r.Lat.Value))
return true;
return false;
}
// To find orientation of ordered triplet (p, q, r).
// The function returns following values
// 0 --> p, q and r are collinear
// 1 --> Clockwise
// 2 --> Counterclockwise
private int orientation(Coords p, Coords q, Coords r)
{
// See https://www.geeksforgeeks.org/orientation-3-ordered-points/
// for details of below formula.
double val = (q.Lat.Value - p.Lat.Value) * (r.Lon.Value - q.Lon.Value) -
(q.Lon.Value - p.Lon.Value) * (r.Lat.Value - q.Lat.Value);
if (val == 0) return 0; // collinear
return (val > 0) ? 1 : 2; // clock or counterclock wise
}
// The main function that returns true if line segment 'p1q1'
// and 'p2q2' intersect.
private Boolean doIntersect(Coords p1, Coords q1, Coords p2, Coords q2)
{
// Find the four orientations needed for general and
// special cases
int o1 = orientation(p1, q1, p2);
int o2 = orientation(p1, q1, q2);
int o3 = orientation(p2, q2, p1);
int o4 = orientation(p2, q2, q1);
// General case
if (o1 != o2 && o3 != o4)
return true;
// Special Cases
// p1, q1 and p2 are collinear and p2 lies on segment p1q1
if (o1 == 0 && onSegment(p1, p2, q1)) return true;
// p1, q1 and q2 are collinear and q2 lies on segment p1q1
if (o2 == 0 && onSegment(p1, q2, q1)) return true;
// p2, q2 and p1 are collinear and p1 lies on segment p2q2
if (o3 == 0 && onSegment(p2, p1, q2)) return true;
// p2, q2 and q1 are collinear and q1 lies on segment p2q2
if (o4 == 0 && onSegment(p2, q1, q2)) return true;
return false; // Doesn't fall in any of the above cases
}
private List<Segment> ComputeAllSegments(List<Coords> coords)
{
List<Segment> segments = new List<Segment>();
int index = 1;
int id_segment = 0;
foreach (Coords coord in coords)
{
for (var i = index; i < coords.Count; i++)
{
segments.Add(new Segment(id_segment, coord, coords[i]));
id_segment++;
}
index++;
}
return segments;
}
private List<Tuple<int, int>> ComputeCrossingSegments(List<Segment> segments)
{
List<Tuple<int, int>> crossing_segments = new List<Tuple<int, int>>();
int index = 1;
int cnt = 0;
foreach (Segment s in segments)
{
for (var i = index; i < segments.Count; i++)
{
if (s.Pt1 == segments[i].Pt1 || s.Pt1 == segments[i].Pt2 || s.Pt2 == segments[i].Pt1 || s.Pt2 == segments[i].Pt2)
{
//Don't consider segments sharing same point
continue;
}
else if (doIntersect(s.Pt1, s.Pt2, segments[i].Pt1, segments[i].Pt2))
{
crossing_segments.Add(new Tuple<int, int>(s.Id, segments[i].Id));
}
}
index++;
}
return crossing_segments;
}
public List<Polygon> ComputeAllPossiblePolygons(List<Coords> coords)
{
List<Polygon> polygons = new List<Polygon>();
if (coords.Count <= 3)
{
return polygons;
}
//Let's compute all crossin segments
all_segments = ComputeAllSegments(coords);
crossing_segments = ComputeCrossingSegments(all_segments);
List<Coords> convex = ConvexHull.ComputeConvexHull(coords);
List<Segment> polygon_segments = new List<Segment>();
for(int i = 0; i < convex.Count-1; i++)
{
polygon_segments.Add(FindSegment(convex[i], convex[i + 1]));
}
polygon_segments.Add(FindSegment(convex.Last(), convex.First()));
Polygon current = new Polygon(polygon_segments);
List<Coords> remaining = new List<Coords>();
foreach (Coords coord in coords)
{
if (!convex.Exists(t=>t.Lon.Value == coord.Lon.Value && t.Lat.Value == coord.Lat.Value))
{
remaining.Add(coord);
}
}
int goal = coords.Count;
FindPolygonsRecursive(goal, current, remaining, ref polygons);
return polygons;
}
public List<List<Coords>> ComputAllPossibleOrdering(List<Coords> coords)
{
List<int> existing_ordering = new List<int>();
List<List<Coords>> all_ordering = new List<List<Coords>>();
List<Polygon> polygons = ComputeAllPossiblePolygons(coords);
foreach (Polygon polygon in polygons)
{
List<Coords> ordering = new List<Coords>();
//Start from any segment
Segment s = polygon.Segments.First();
ordering.Add(s.Pt1);
ordering.Add(s.Pt2);
polygon.Segments.Remove(s);
Coords last_coords = s.Pt2;
for (int i = 0; i < coords.Count - 1; i++)
{
if (polygon.Segments.Exists(t => t.Pt1 == last_coords))
{
Segment seg = polygon.Segments.First(t => t.Pt1 == last_coords);
if(seg.Pt2 != s.Pt1)
{
ordering.Add(seg.Pt2);
}
last_coords = seg.Pt2;
polygon.Segments.Remove(seg);
}
else if (polygon.Segments.Exists(t => t.Pt2 == last_coords))
{
Segment seg = polygon.Segments.First(t => t.Pt2 == last_coords);
if (seg.Pt1 != s.Pt1)
{
ordering.Add(seg.Pt1);
}
last_coords = seg.Pt1;
polygon.Segments.Remove(seg);
}
}
int nbre = 0;
for (int i = 0; i < ordering.Count; i++)
{
nbre += i * ordering[i].Id;
}
if (!existing_ordering.Contains(nbre))
{
existing_ordering.Add(nbre);
all_ordering.Add(ordering);
}
}
return all_ordering;
}
private void FindPolygonsRecursive(int goal, Polygon current, List<Coords> remaining, ref List<Polygon> polygons)
{
if (current.Segments.Count == goal)
{
if (!polygons.Contains(current))
{
polygons.Add(current);
}
return;
}
foreach (Coords p in remaining)
{
List<Coords> without_p = new List<Coords>();
foreach (Coords coord in remaining)
{
if (coord != p) { without_p.Add(coord); }
}
foreach (Segment s in current.Segments)
{
Segment new_segment1 = FindSegment(s.Pt1, p);
Segment new_segment2 = FindSegment(s.Pt2, p);
if (!Intersect(new_segment1, current) && !Intersect(new_segment2, current))
{
List<Segment> next_segments = new List<Segment>();
foreach(Segment previous_segment in current.Segments)
{
if (previous_segment != s)
{
next_segments.Add(previous_segment);
}
}
Polygon next = new Polygon(next_segments);
next.Segments.Add(new_segment1);
next.Segments.Add(new_segment2);
if (polygons.Count > 200)
{
return;
}
FindPolygonsRecursive(goal, next, without_p, ref polygons);
}
}
}
}
private Segment FindSegment(Coords pt1, Coords pt2)
{
if (all_segments.Exists(t => t.Pt1 == pt1 && t.Pt2 == pt2))
{
return all_segments.First(t => t.Pt1 == pt1 && t.Pt2 == pt2);
}
if (all_segments.Exists(t => t.Pt1 == pt2 && t.Pt2 == pt1))
{
return all_segments.First(t => t.Pt1 == pt2 && t.Pt2 == pt1);
}
return null;
}
private bool Intersect(Segment s, Polygon polygon)
{
foreach (Segment segment in polygon.Segments)
{
int first = (s.Id > segment.Id) ? segment.Id: s.Id;
int second = (s.Id > segment.Id) ? s.Id : segment.Id;
if (crossing_segments.Contains(new Tuple<int, int>(first, second)))
{
return true;
}
}
return false;
}
}
public static class ConvexHull
{
public static List<Coords> ComputeConvexHull(List<Coords> points, bool sortInPlace = false)
{
if (!sortInPlace)
points = new List<Coords>(points);
points.Sort((a, b) =>
a.Lon.Value == b.Lon.Value ? a.Lat.Value.CompareTo(b.Lat.Value) : (a.Lon.Value > b.Lon.Value ? 1 : -1));
// Importantly, DList provides O(1) insertion at beginning and end
CircularList<Coords> hull = new CircularList<Coords>();
int L = 0, U = 0; // size of lower and upper hulls
// Builds a hull such that the output polygon starts at the leftmost Vector2.
for (int i = points.Count - 1; i >= 0; i--)
{
Coords p = points[i], p1;
// build lower hull (at end of output list)
while (L >= 2 && (p1 = hull.Last).Sub(hull[hull.Count - 2]).Cross(p.Sub(p1)) >= 0)
{
hull.PopLast();
L--;
}
hull.PushLast(p);
L++;
// build upper hull (at beginning of output list)
while (U >= 2 && (p1 = hull.First).Sub(hull[1]).Cross(p.Sub(p1)) <= 0)
{
hull.PopFirst();
U--;
}
if (U != 0) // when U=0, share the Vector2 added above
hull.PushFirst(p);
U++;
}
hull.PopLast();
return hull;
}
private static Coords Sub(this Coords a, Coords b)
{
Coords coord = new Coords();
coord.Lon = a.Lon.Value - b.Lon.Value;
coord.Lat = a.Lat.Value - b.Lat.Value;
return coord;
}
private static float Cross(this Coords a, Coords b)
{
return (float)(a.Lon.Value * b.Lat.Value - a.Lat.Value * b.Lon.Value);
}
private class CircularList<T> : List<T>
{
public T Last
{
get
{
return this[this.Count - 1];
}
set
{
this[this.Count - 1] = value;
}
}
public T First
{
get
{
return this[0];
}
set
{
this[0] = value;
}
}
public void PushLast(T obj)
{
this.Add(obj);
}
public T PopLast()
{
T retVal = this[this.Count - 1];
this.RemoveAt(this.Count - 1);
return retVal;
}
public void PushFirst(T obj)
{
this.Insert(0, obj);
}
public T PopFirst()
{
T retVal = this[0];
this.RemoveAt(0);
return retVal;
}
}
}