Sort algorithm to create a polygon from points with only right angles

Question

Given a set of (x, y) coordinates in some random order, can they be sorted so that a polygonal path can be drawn with only 90^o internal or external angles.

It is know that such a path exists, but it is not know what order the edge points of the polygon need to be connected.

The closest solutions readily findable in SO are:

• Algorithm to create a polygon from points

• Algorithm to create polygon

Both of these use polar coordination to order the points, and will produce a star like polygon, for which only some of the corners are 90^o angles.

[NOTE This is a reposting of a deleted question: Sort algorithm to create a polygon from points with only right angle. I had developed a solution and went to post it only to find that the question had been deleted. I am reposting it here because others may find it useful.]

Answer 1

I created the deleted question, and I found a similar solution myself.

After trying my solution with many different polygon formats; in some cases, it didn't work. So I tried yours, same result.

Try these points:

points = [[-72, 0], [53, -36], [-72, 17], [-36, 35], [53, 53], [71, 0], [-54, 17], [53, 0], [-54, 0], [-54, 35], [-36, 53], [-54, -36], [71, 17], [53, 17]]

<svg xmlns="http://www.w3.org/2000/svg" width="380" height="200.0" viewBox="-26.0 -40.0 75.0 97.0"><path stroke="black" stroke-width="7" fill="none" d="M -72.0 0.0 L -72.0 17.0 L -54.0 17.0 L -54.0 0.0 L 53.0 0.0 L 53.0 -36.0 L -54.0 -36.0 L -54.0 35.0 L -36.0 35.0 L -36.0 53.0 L 53.0 53.0 L 53.0 17.0 L 71.0 17.0 L 71.0 0.0 Z"></path><path stroke="red" stroke-width="2" fill="none" d="M -72 0 L -72 0 L -54 0 L -54 -36 L 53 -36 L 53 0 L 71 0 L 71 17 L 53 17 L 53 53 L -36 53 L -36 35 L -54 35 L -54 17 L -72 17 Z"></path><circle fill="blue" cx="-5.571428571428571" cy="12.285714285714286" r="7"></circle></svg>

The polar sort is giving good results here, but not your algorithm.

You take the closest point (when there are multiple on the same axis) and this is wrong in some cases.

Here's the workaround that I found, based on your rightSort function.

First, before adding a new possible value to the sorted array, I'm checking if it does cross over a past sorted point. In this case it will not add it.

Second, I added an hashset that save what path is leading to a dead end because it can get stuck in an infinite loop.

Here's a rust solution, assuming the points are integers, and x/y aligned:

use glam::IVec2;
use std::collections::HashSet;
use anyhow::{anyhow, Result};

fn right_sort(mut points: Vec<IVec2>) -> Result<Vec<IVec2>> {
    let mut sorted = vec![points[0].clone()];
    let next_point = closest_point(sorted.last().unwrap(), &points);
    let mut x = (sorted.last().unwrap().x - next_point.x) == 0;
    let mut dead_ends = HashSet::new();

    let mut popped = vec![];
    points.remove(0);
    while points.len() > 0 {
        let mut ndxes = vec![];
        if x {
            for ndx in 0..points.len() {
                if (points[ndx].x - sorted.last().unwrap().x).abs() == 0
                    && !check_points_on_line(sorted.last().unwrap(), &points[ndx], &sorted)
                    && !is_dead_end(&sorted, &points[ndx], &dead_ends) {
                    ndxes.push(ndx);
                }
            }
            if ndxes.len() == 0 {
                let hashcode = get_path_hashcode(&sorted);
                dead_ends.insert(hashcode);
                popped.push(sorted.pop().unwrap());
                if sorted.len() == 0 {
                    return Err(anyhow!("No solution found!"));
                }
                x = false;
            }
            else {
                let mut closest_dist = (points[ndxes[0]].y - sorted.last().unwrap().y).abs();
                let mut ndx_closest = ndxes[0];
                for ndx in &ndxes[1..] {
                    if (points[*ndx].y - sorted.last().unwrap().y).abs() < closest_dist {
                        closest_dist = (points[*ndx].y - sorted.last().unwrap().y).abs();
                        ndx_closest = *ndx;
                    }
                }
                sorted.push(points.remove(ndx_closest));
                x = false;
                if popped.len() > 0 {
                    points.append(&mut popped);
                    popped = vec![];
                }
            }
        }
        else {
            for ndx in 0..points.len() {
                if (points[ndx].y - sorted.last().unwrap().y).abs() == 0
                    && !check_points_on_line(sorted.last().unwrap(), &points[ndx], &sorted)
                    && !is_dead_end(&sorted, &points[ndx], &dead_ends) {
                    ndxes.push(ndx);
                }
            }
            if ndxes.len() == 0 {
                let hashcode = get_path_hashcode(&sorted);
                dead_ends.insert(hashcode);
                popped.push(sorted.pop().unwrap());
                if sorted.len() == 0 {
                    return Err(anyhow!("No solution found!"));
                }
                x = true;
            }
            else {
                let mut closest_dist = (points[ndxes[0]].x - sorted.last().unwrap().x).abs();
                let mut ndx_closest = ndxes[0];
                for ndx in &ndxes[1..] {
                    if (points[*ndx].x - sorted.last().unwrap().x).abs() < closest_dist {
                        closest_dist = (points[*ndx].x - sorted.last().unwrap().x).abs();
                        ndx_closest = *ndx;
                    }
                }
                sorted.push(points.remove(ndx_closest));
                x = true;
                if popped.len() > 0 {
                    points.append(&mut popped);
                    popped = vec![];
                }
            }
        }
    }
    if popped.len() > 0 {
        sorted.append(&mut popped);
    }

    Ok(sorted)
}

fn closest_point(focus: &IVec2, points: &Vec<IVec2>) -> IVec2 {
    let mut closest_point;
    if focus.x != points[0].x && focus.y != points[0].y {
        closest_point = &points[0];
    }
    else {
        closest_point = &points[1];
    }
    let mut closest_dist = (focus.x - closest_point.x).pow(2) + (focus.y - closest_point.y).pow(2);
    for point in points {
        if focus.x != point.x && focus.y != point.y {
            let dist = (focus.x - point.x).pow(2) + (focus.y - point.y).pow(2);
            if dist < closest_dist {
                closest_dist = dist;
                closest_point = point;
            }
        }
    }
    closest_point.clone()
}

fn check_points_on_line(p1: &IVec2, p2: &IVec2, sorted_points: &Vec<IVec2>) -> bool {
    for point in sorted_points {
        if point.eq(p1) {
            continue ;
        }
        if p1.x == p2.x && point.x == p1.x {
            if (p1.y < point.y && p2.y > point.y)
                || (p2.y < point.y && p1.y > point.y) {
                return true;
            }
        }
        else if p1.y == p2.y && point.y == p1.y {
            if (p1.x < point.x && p2.x > point.x)
                || (p2.x < point.x && p1.x > point.x) {
                return true;
            }
        }
    }
    false
}

fn is_dead_end(sorted: &Vec<IVec2>, focus: &IVec2, dead_ends: &HashSet<String>) -> bool {
    let mut hashcode = get_path_hashcode(sorted);
    hashcode.push_str(&format!("{}{}", focus.x, focus.y));
    dead_ends.contains(&hashcode)
}

fn get_path_hashcode(path: &Vec<IVec2>) -> String {
    let mut hashcode = String::new();
    for point in path {
        hashcode.push_str(&format!("{}{}", point.x, point.y));
    }
    hashcode
}

Answer 2

To sort randomly ordered (x, y) rectalinear points into rectalinear polygonal order:

1. find the center of the points
2. find the remotest point
3. find the nearest point to the remotest point
4. find the angle between the remotest point and nearest remote point and the x/y axis (probably could be any two "nearest" points but remotest nearest points should reduce the likelihood of any ambiguity)
5. rotate all points to be x-y axis aligned
6. pick any point as a start point as the first stepping point
7. find the nearest point as the next point
8. if the stepping point and the next point are x-axis aligned look for the next nearest y-axis aligned point
9. if the stepping point and the next point are y-axis aligned look for the next nearest x-axis aligned point
10. if there is no next axis aligned point, back track one point at a time, temporarily removing the back tracked points from available next points until another next back tracked axis aligned point is found and then add the back tracked points back to the available next points (back tracking is necessary because it is possible to get into an enclave with no path out, but that is not a valid polygon)
11. make the next point the stepping point
12. alternate between x and y axis aligned next nearest points
13. repeat from 10 until all points are used
14. rotate points back to the their original alignment

The code below is a rough implementation in python. It will produce a number of SVG files for comparison.

points = [(156.40357183517773, 23.19316100057085), (83.97002318399646, 188.27914171909507), (518.4511031561435, 60.897074118366035), (799.3826769425817, 214.44658030407507), (304.1247347870089, -2.8540656494687013), (593.7387174567936, 199.93582818685073), (773.3354502925422, 66.72541735224387), (625.6142873407109, 92.7726440022834), (428.65273673826925, 127.50227953566946), (379.41234908765887, 136.184688419016), (446.0175545049623, 225.98305483689026), (448.871620154431, 530.1077896238992), (509.768694272797, 11.65668646775564), (373.58400585378104, 391.06903555541453), (602.4211263401401, 249.17621583746111), (182.45079848521726, 170.91432395240204), (616.9318784573643, 43.53225635167299), (165.08598071852424, 72.43354865118125), (312.80714367035546, 46.3863220011417), (225.86284290194985, 417.1162622054541), (399.63123250382057, 538.7901985072457), (66.60520541730344, 89.79836641787429)]

def genSVG(points):
    path = "M " + str(points[0][0]) + " " + str(points[0][1]) + " "
    minX = points[0][0]
    minY = points[0][1]
    maxX = minX
    maxY = minY
    for point in points[1:]:
        path += "L " + str(point[0]) + " " + str(point[1]) + " "
        if point[0] < minX:
            minX = point[0]
        elif point[0] > maxX:
            maxX = point[0]
        if point[1] < minY:
            minY = point[1]
        elif point[1] > maxY:
            maxY = point[1]
    path += "Z"
    path = '<path fill="grey" d="' + path + '"/>'

    viewbox = ' viewbox="' + str(minX-1) + ' ' + str(minY-1) + ' ' + str(maxX+1) + ' ' + str(maxY+1) + '"'
    
    width = ' width="' + str((maxX - minX + 2)) + '"'
    height = ' height="' + str((maxY - minY + 2)) + '"'

    return '<svg ' + 'xmlns="http://www.w3.org/2000/svg"' + width + height + viewbox + '>' + path + '</svg>'

def genSVGover(points, overs, center):
    path = "M " + str(points[0][0]) + " " + str(points[0][1]) + " "
    minX = points[0][0]
    minY = points[0][1]
    maxX = minX
    maxY = minY
    for point in points[1:]:
        path += "L " + str(point[0]) + " " + str(point[1]) + " "
        if point[0] < minX:
            minX = point[0]
        elif point[0] > maxX:
            maxX = point[0]
        if point[1] < minY:
            minY = point[1]
        elif point[1] > maxY:
            maxY = point[1]
    path += "Z"
    path = '<path stroke="black" stroke-width="7" fill="none" d="' + path + '"/>'

    viewbox = ' viewbox="' + str(minX-4) + ' ' + str(minY-4) + ' ' + str(maxX+4) + ' ' + str(maxY+4) + '"'
    
    width = ' width="' + str((maxX - minX + 8)) + '"'
    height = ' height="' + str((maxY - minY + 8)) + '"'

    over = "M " + str(overs[0][0]) + " " + str(overs[0][1]) + " "
    for point in overs:
        over += "L " + str(point[0]) + " " + str(point[1]) + " "
    over += "Z"
    over = '<path stroke="red" stroke-width="2" fill="none" d="' + over + '"/>'
        
    return '<svg ' + 'xmlns="http://www.w3.org/2000/svg"' + width + height + viewbox + '>' + path + over + '<circle fill="blue" cx="' + str(center[0]) + '" cy="' + str(center[1]) + '" r="7" />' + '</svg>'

import math
def rotate(points, theta):
    rotated = []
    cosTheta = math.cos(theta)
    sinTheta = math.sin(theta)
    for point in points:
        rotated.append(( cosTheta * point[0] + sinTheta * point[1], -sinTheta * point[0] + cosTheta * point[1] ))
    return rotated

def closest(focus, points):
    if ( points[0] != focus ):
        closestPoint = points[0]
    else:
        closestPoint = points[1]
    closestDist = ( focus[0] - closestPoint[0] )**2 + ( focus[1] - closestPoint[1] )**2
    for point in points:
        if point != focus:
            dist = ( focus[0] - point[0] )**2 + ( focus[1] - point[1] )**2
            if dist < closestDist:
                closestDist = dist
                closestPoint = point
    return closestPoint

def rotangle(points):
    focus = remotest(points)
    closestPoint = closest(focus, points)
    if abs(focus[0] - closestPoint[0]) < tolerance or abs(focus[1] - closestPoint[1]) < tolerance:
        return 0
    else:
        return math.atan2(focus[1] - closestPoint[1], focus[0] - closestPoint[0])

tolerance = 0.000000000001
def rightSort(points):
    sorted = [ points[0] ]
    nextPoint = closest(sorted[-1], points)
    x = abs( sorted[-1][0] - nextPoint[0]) < tolerance
    popped = []
    del points[0]
    while len(points) > 0:
        ndxes = []
        if x:
            for ndx in range(len(points)):
                if abs(points[ndx][0] - sorted[-1][0]) < tolerance:
                    ndxes.append(ndx)
            if len(ndxes) == 0:
                popped.append(sorted.pop())
                x = False
            else:
                closestDist = abs(points[ndxes[0]][1] - sorted[-1][1])
                ndxClosest = ndxes[0]
                for ndx in ndxes[1:]:
                    if abs(points[ndx][1] - sorted[-1][1]) < closestDist:
                        ndxClosest = ndx
                sorted.append(points[ndxClosest])
                del points[ndxClosest]
                x = False
                if popped:
                    points += popped
                    popped = []
        else:
            for ndx in range(len(points)):
                if abs(points[ndx][1] - sorted[-1][1]) < tolerance:
                    ndxes.append(ndx)
            if len(ndxes) == 0:
                popped.append(sorted.pop())
                x = True
            else:
                closestDist = abs(points[ndxes[0]][0] - sorted[-1][0])
                ndxClosest = ndxes[0]
                for ndx in ndxes[1:]:
                    if abs(points[ndx][0] - sorted[-1][0]) < closestDist:
                        ndxClosest = ndx
                sorted.append(points[ndxClosest])
                del points[ndxClosest]
                x = True
                if popped:
                    points += popped
                    popped = []
    if popped:
        sorted += popped
    return sorted

def center(points):
    return ( sum(point[0] for point in points) / len(points),
                sum(point[1] for point in points) / len(points) )

def remotest(points):
    centerPoint = center(points)
    print( "center", centerPoint )
    remotestPoint = points[0]
    remotestDist = ( centerPoint[0] - remotestPoint[0] )**2 + ( centerPoint[1] - remotestPoint[1] )**2
    for point in points[1:]:
        dist = ( centerPoint[0] - point[0] )**2 + ( centerPoint[1] - point[1] )**2
        if dist > remotestDist:
            remotestDist = dist
            remotestPoint = point
    print( "remotest", remotestPoint )
    return remotestPoint

def squaredPolar(point, centerPoint):
    return ( math.atan2(point[1] - centerPoint[1], point[0] - centerPoint[0]),
            ( point[0] - centerPoint[0] )**2 + ( point[1] - centerPoint[1] )**2 )

def polarSort(points):
    centerPoint = center(points)
    presorted = []
    for point in points:
        presorted.append(( squaredPolar(point, centerPoint), point ))
    presorted.sort()
    sorted = []
    for point in presorted:
        sorted.append(point[1])
    return sorted
    
htmlFile = open("polygon.html", "w")
htmlFile.write("<html><body>")
htmlFile.write(genSVG(points))
htmlFile.write("</body></html>")
htmlFile.close()

angle = rotangle(points)
print( "angle", angle * 180 / math.pi )

htmlFile = open("rightgon.html", "w")
htmlFile.write("<html><body>")
htmlFile.write(genSVGover(rotate(rightSort(rotate(points, angle)), -angle), polarSort(points), center(points)))
htmlFile.write("</body></html>")
htmlFile.close()

htmlFile = open("polargon.html", "w")
htmlFile.write("<html><body>")
htmlFile.write(genSVG(polarSort(points)))
htmlFile.write("</body></html>")
htmlFile.close()

The image below is an unsorted points "polygon".

<svg xmlns="http://www.w3.org/2000/svg" width="734.7774715252783" height="543.6442641567144" viewbox="65.60520541730344 -3.8540656494687013 800.3826769425817 539.7901985072457"><path fill="grey" d="M 156.40357183517773 23.19316100057085 L 83.97002318399646 188.27914171909507 L 518.4511031561435 60.897074118366035 L 799.3826769425817 214.44658030407507 L 304.1247347870089 -2.8540656494687013 L 593.7387174567936 199.93582818685073 L 773.3354502925422 66.72541735224387 L 625.6142873407109 92.7726440022834 L 428.65273673826925 127.50227953566946 L 379.41234908765887 136.184688419016 L 446.0175545049623 225.98305483689026 L 448.871620154431 530.1077896238992 L 509.768694272797 11.65668646775564 L 373.58400585378104 391.06903555541453 L 602.4211263401401 249.17621583746111 L 182.45079848521726 170.91432395240204 L 616.9318784573643 43.53225635167299 L 165.08598071852424 72.43354865118125 L 312.80714367035546 46.3863220011417 L 225.86284290194985 417.1162622054541 L 399.63123250382057 538.7901985072457 L 66.60520541730344 89.79836641787429 Z"/></svg>

The image below is the rendering of one output file. It shows:

• blue dot is the center of the (x, y) coordinates
• red polygon is the polar sorted polygon
• black polygon is the right angle sorted polygon

<svg xmlns="http://www.w3.org/2000/svg" width="740.7774715252784" height="549.6442641567145" viewbox="62.60520541730345 -6.854065649468694 803.3826769425818 542.7901985072458"><path stroke="black" stroke-width="7" fill="none" d="M 156.40357183517776 23.19316100057085 L 165.08598071852424 72.43354865118125 L 66.60520541730345 89.7983664178743 L 83.97002318399647 188.2791417190951 L 182.4507984852173 170.91432395240207 L 225.86284290194988 417.1162622054542 L 373.5840058537811 391.0690355554146 L 399.63123250382057 538.7901985072458 L 448.87162015443107 530.1077896238993 L 379.41234908765887 136.184688419016 L 428.65273673826937 127.50227953566947 L 446.01755450496233 225.9830548368903 L 593.7387174567937 199.93582818685076 L 602.4211263401402 249.17621583746114 L 799.3826769425818 214.44658030407507 L 773.3354502925423 66.72541735224388 L 625.614287340711 92.7726440022834 L 616.9318784573644 43.532256351673 L 518.4511031561435 60.89707411836606 L 509.76869427279706 11.656686467755648 L 312.8071436703555 46.3863220011417 L 304.1247347870089 -2.8540656494686942 Z"/><path stroke="red" stroke-width="2" fill="none" d="M 182.45079848521726 170.91432395240204 L 182.45079848521726 170.91432395240204 L 66.60520541730344 89.79836641787429 L 165.08598071852424 72.43354865118125 L 156.40357183517773 23.19316100057085 L 379.41234908765887 136.184688419016 L 312.80714367035546 46.3863220011417 L 304.1247347870089 -2.8540656494687013 L 428.65273673826925 127.50227953566946 L 509.768694272797 11.65668646775564 L 518.4511031561435 60.897074118366035 L 616.9318784573643 43.53225635167299 L 625.6142873407109 92.7726440022834 L 773.3354502925422 66.72541735224387 L 799.3826769425817 214.44658030407507 L 593.7387174567936 199.93582818685073 L 602.4211263401401 249.17621583746111 L 446.0175545049623 225.98305483689026 L 448.871620154431 530.1077896238992 L 399.63123250382057 538.7901985072457 L 373.58400585378104 391.06903555541453 L 225.86284290194985 417.1162622054541 L 83.97002318399646 188.27914171909507 Z"/><circle fill="blue" cx="409.6874424591604" cy="177.00212769986794" r="7" /></svg>