Algorithm to close a polygon

Question

I have a part of perimeter of a polygon and need to close it.Please refer this image

As I can see there is only one unique way to close the polygon without dividing the polygon and without the edges intersecting.

And the closing edges would be b->c,d->e,f->g,h->a

Is there any algo to achieve this?

I can think of only one brute force method, try every possible combination and check if it forms a closed polygon(Any good algos to check if it is closed polygon?)

Is there any better way or a known algorithm?

Note: The vertices should be connected by single straight lines only and polygon is not necessarily convex

Also, You can safely assume that these segments always form a polygon because I get these line segments from a polygon and Im trying to recreate the polygon

Answer 1

I think that in "well-behaved" (small gaps, not too irregular shape, etc.) cases, one might get away with following approach. The idea is to assume that the solution (particular permutation of the input line segments which are then assumed to be connected with straight lines) minimizes the length of the resulting MultiLineString defining the boundary of the polygon of interest.

To tackle this problem, the implementation below uses the 2-opt heuristic to the traveling salesman problem. It proceeds in following steps:

1. the set of vertices is defined as the union of the endpoints of all input line segments
2. it tries to connect these points in order to minimize the total length of the resulting MultiLineString under the constraint that the points belonging to the same input line segment are always connected (i.e., the 2-opt algorithm is allowed to split only edges connecting different line segments - this is handled by the extra if condition in the main double for-loop).

The result is then:

import logging
import random
import sys
from shapely.geometry import LineString, Polygon
from shapely.ops import polygonize, linemerge

#prevent shapely from showing an error message on is_valid tests  
logger = logging.getLogger()
logger.setLevel(logging.ERROR)

#input lines (LineStrings)
lines = [
    [(3.15,3.94), (4.06,3.91), (4.27,3.49)],
    [(0.84,2.99), (0.97,3.67), (1.68,3.91), (2.68,3.92)],
    [(4.46,3.23), (5.12,2.97), (4.60,2.00)],
    [(4.13,1.44), (4.41,0.68), (1.14,1.99)]
]
random.shuffle(lines)

N, pnts = 0, []
pnt2line = {}
for line_id, line in enumerate(lines):
    #for each line, save its endpoints and remember
    #to which line each point belongs
    for pnt in [line[0], line[-1]]:
        pnt2line[N] = line_id
        pnts.append(pnt)
        N += 1

#as initial guess, try to connect these points sequentially
route = [i for i in range(0, N)]

def nrm_idx(N, idx):
    return (N + idx) % N

def get_polygon(route):
    #for given route, attempt to construct the resulting polygon
    segments = []
    m = len(route)
    for idx in range(0, m):
        i, j = route[idx], route[nrm_idx(m, idx+1)]
        if pnt2line[i] == pnt2line[j]:
            #if two consecutive points belong to the same line, consider this line
            segments.append(lines[pnt2line[i]])
        else:
            #otherwise, connect these points with a straight line
            segments.append([pnts[i], pnts[j]])

    return Polygon(linemerge(segments))

def get_weight(route):
    P = get_polygon(route)
    return P.length if P.is_valid else sys.maxsize

def edge_is_fixed(pnt_i, pnt_j):
    #check if an edge specified by point pnt_i/pnt_j can be dissected or not
    #in the latter case, the points belong to the same line/line segment
    return (pnt2line[pnt_i] == pnt2line[pnt_j])

def opt_swap(route, i, k):
    #perform 2-opt swap
    return route[0:i] + route[i:k+1][::-1] + route[k+1:]

flag = True
while flag:
    flag = False
    best_weight = get_weight(route)

    for i in range(0, N-1):
        for k in range(i+1, N):

            if edge_is_fixed(route[nrm_idx(N, i-1)], route[i]) or edge_is_fixed(route[k], route[nrm_idx(N, k+1)]):
                continue

            new_route = opt_swap(route, i, k)
            weight = get_weight(new_route)

            if weight < best_weight:
                route = new_route[:]
                best_weight = weight
                flag = True

P = get_polygon(route)
for x, y in P.exterior.coords:
    print(x, y)

For your input (approximated), the result is then indeed: