pythonalgorithmpolygonconcave

Ordering points that roughly follow the contour of a line to form a polygon

I have a list of points (x,y). A line is drawn somewhere 'inside' these points (within the area they form). I have been trying different algorithms for ordering these points to create a polygon, with these constraints:

  1. The polygon is not self-intersecting.
  2. The resulting polygon is not necessarily convex; a convex hull will not suffice
  3. The polygon needs to go through all of the points. I have tried the alphashape concave hull algorithm, but I was not able to get it to work (it ignored some points, and I wasn't sure what to tune the alpha value to)
  4. The resulting polygon will not intersect with the drawn line. I think this is relevant to remove a lot of the ambiguity for concave polygons

Some clarification for how the points are ordered: they roughly follow the profile of a contoured line, which is unknown and different from the drawn line. The drawn line will be intrinsically quite similar to this unspecified line however.

I have also tried the very simplistic angle sort algorithm, which unsurprisingly does not work. However, running the algorithm at differnt points along the line gave more promising results, but I am not sure how to use the data.

Example of problem

Above is a screenshot of the following test program I cobbled together:

import pygame
import numpy as np
from scipy.interpolate import splprep, splev
import sys
import random

# Pygame setup
pygame.init()
screen = pygame.display.set_mode((1000, 1000))
pygame.display.set_caption("Smoothed Curve with Draggable Points")
clock = pygame.time.Clock()

# Colors
WHITE = (255, 255, 255)
BLUE  = (0, 0, 255)
RED   = (255, 0, 0)
GREEN = (0, 200, 0)
BLACK = (0, 0, 0)

drawing = False
points = []
smoothed = []
normals = []

# Dot settings
DOT_RADIUS = 6
DRAG_RADIUS = 10

# Dot class for draggable points
class Dot:
    def __init__(self, pos):
        self.x, self.y = pos
        self.dragging = False

    def pos(self):
        return (int(self.x), int(self.y))

    def is_mouse_over(self, mx, my):
        return (self.x - mx)**2 + (self.y - my)**2 <= DRAG_RADIUS**2

    def start_drag(self):
        self.dragging = True

    def stop_drag(self):
        self.dragging = False

    def update_pos(self, mx, my):
        if self.dragging:
            self.x, self.y = mx, my

# crude test arrangement of dots.
random_dots = [Dot(tup) for tup in [(159, 459),(133, 193),(286, 481),(241, 345),(411, 404),(280, 349),(352, 471),(395, 361),(85, 390),(203, 321),(41, 281),(58, 348),(175, 275),(75, 185),(385, 443),(44, 219),(148, 229),(215, 477),(338, 339),(122, 430)]]

def downsample(points, step=2):
    return points[::step]

def smooth_curve(points, smoothness=150):
    if len(points) < 4:
        return [], []
    
    points = downsample(points, step=2)
    x, y = zip(*points)
    try:
        tck, u = splprep([x, y], s=100.0, k=3)
        u_new = np.linspace(0, 1, smoothness)

        # Evaluate curve points
        x_vals, y_vals = splev(u_new, tck)

        # Derivatives: dx/du, dy/du
        dx_du, dy_du = splev(u_new, tck, der=1)

        gradients = np.array(dy_du) / np.array(dx_du)  # dy/dx

        return list(zip(x_vals, y_vals)), list(zip(dx_du, dy_du))
    except:
        return [], []

def draw_normals(screen, smoothed, derivatives, spacing=10, length=30):
    for i in range(0, len(smoothed), spacing):
        x, y = smoothed[i]
        dx, dy = derivatives[i]

        # Normal vector is (-dy, dx)
        nx, ny = -dy, dx
        norm = np.hypot(nx, ny)
        if norm == 0:
            continue
        nx /= norm
        ny /= norm

        x1 = x - nx * length / 2
        y1 = y - ny * length / 2
        x2 = x + nx * length / 2
        y2 = y + ny * length / 2

        pygame.draw.line(screen, GREEN, (x1, y1), (x2, y2), 2)

def generate_random_dots(n=20, w=800, h=600):
    return [Dot((random.randint(0, w), random.randint(0, h))) for _ in range(n)]

# Main loop
dragged_dot = None
running = True
while running:
    screen.fill(WHITE)
    mx, my = pygame.mouse.get_pos()

    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

        elif event.type == pygame.MOUSEBUTTONDOWN:
            drawing = True
            points = []
            smoothed = []
            normals = []

            # Check for draggable dot
            for dot in random_dots:
                if dot.is_mouse_over(mx, my):
                    dragged_dot = dot
                    dot.start_drag()
                    drawing = False
                    break

        elif event.type == pygame.MOUSEBUTTONUP:
            drawing = False
            if dragged_dot:
                dragged_dot.stop_drag()
                dragged_dot = None
            else:
                smoothed, normals = smooth_curve(points)

        elif event.type == pygame.MOUSEMOTION:
            if dragged_dot:
                dragged_dot.update_pos(mx, my)
            elif drawing:
                points.append((mx, my))

        elif event.type == pygame.KEYDOWN:
            if event.key == pygame.K_p:
                random_dots = generate_random_dots()
            if event.key == pygame.K_l:
                for i in random_dots: print(f"({i.x}, {i.y})")

    # Draw random draggable dots
    for dot in random_dots:
        pygame.draw.circle(screen, BLACK, dot.pos(), DOT_RADIUS)

    # Draw raw stroke
    if len(points) > 1:
        pygame.draw.lines(screen, BLUE, False, points, 1)

    # Draw smoothed curve
    if len(smoothed) > 1:
        pygame.draw.lines(screen, RED, False, smoothed, 3)

    # Draw normals
    if smoothed and normals:
        draw_normals(screen, smoothed, normals)

    pygame.display.flip()
    clock.tick(60)

pygame.quit()
sys.exit()

What algorithm could I use to accomplish this?

EDIT:

Some details about the diagram. The black dots are the input points to form the polygon. The red is the drawn line, which is a smoothed version of the blue one. The green lines show some equidistant test points on the drawn line. The required output order is the order in which the points have to be joined up to form the polygon, i.e. A,B,C,D for a rectangle. The input order of the points is not ordered in any meaningful way.

Solution

  • For completeness I show a geometric implementation that is written in vectorised Numpy. Your code is in Python so I also use Python. For a few reasons, this approach benefits from decimating the interior line down to a low point count.

    import numpy as np
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection


def polysort(
    dots: np.ndarray, line: np.ndarray,
) -> tuple[
    np.ndarray,  # order index
    np.ndarray,  # line segment centroids
    np.ndarray,  # line segment centroids for each dot
]:
    # Between every line point and the next one, get segment centroids
    bounds = np.lib.stride_tricks.sliding_window_view(line, window_shape=2, axis=0)
    centroids = bounds.mean(axis=-1)

    # Get basis vectors for each segment
    basis = np.einsum(
        # i: point index
        # j: x,y
        # k: prev bound, next bound
        # m: u, v
        # n: parallel/normal
        'ijk,k,jmn->imn',
        bounds,
        # Difference between current and next bound, approximately normalised
        (-1e-4, +1e-4),
        # Transform to u, v space
        (
            ((1,  0),
             (0, -1)),
            ((0,  1),
             (1,  0)),
        ),
    )

    # Deltas between each dot and each line point,
    # the index of the closest centroid,
    # and the closest centroid
    diffs = dots[:, np.newaxis] - centroids[np.newaxis]
    idx = np.vecdot(diffs, diffs, axis=-1).argmin(axis=-1)
    chosen_centroids = centroids[idx]

    # Projection of each dot into its segment's uv space
    u, v = np.einsum('ik,ijk->ji', dots - chosen_centroids, basis[idx])
    chirality = v > 0

    # Sort by chirality, then segment index, then segment projection
    keys = np.stack((u, idx, chirality))
    keys[:2, ~chirality] *= -1
    return np.lexsort(keys), centroids, chosen_centroids


def plot(
    dots: np.ndarray, line: np.ndarray, order: np.ndarray,
    centroids: np.ndarray, chosen_centroids: np.ndarray,
) -> plt.Figure:
    fig, ax = plt.subplots()
    ax.scatter(*dots.T)
    ax.scatter(*centroids.T)
    ax.plot(*line.T)

    ax.add_collection(LineCollection(
        np.stack((chosen_centroids, dots), axis=1)[order],
        edgecolors='blue', alpha=0.2,
    ))

    for o, (dx, dy) in enumerate(dots[order]):
        ax.text(dx+5, dy+5, str(o))

    return fig


def demo() -> None:
    # Sample data
    dots = np.array((
        (159, 459), (133, 193), (286, 481), (241, 345),
        (411, 404), (280, 349), (352, 471), (395, 361),
        ( 85, 390), (203, 321), ( 41, 281), ( 58, 348),
        (175, 275), ( 75, 185), (385, 443), ( 44, 219),
        (148, 229), (215, 477), (338, 339), (122, 430),
    ))
    line = np.array((
        ( 94, 209),
        (117, 329),
        (200, 400),
        (287, 419),
        (375, 387),
    ))

    order, centroids, chosen_centroids = polysort(dots=dots, line=line)

    plot(
        dots=dots, line=line, order=order,
        centroids=centroids, chosen_centroids=chosen_centroids,
    )
    plt.show()


if __name__ == '__main__':
    demo()

    diagram