I need an algorithm to calculate the outer boundary of a polygon

Question

I have a polygon that can have an arbitrary number of self-intersections defined as an array of boundary points (two consecutive points form an edge in the polygon boundary). I want to calculate just the external boundary of the polygon (i.e. ignoring all the holes generated by self-intersections). (Not the convex hull, but I want the outer envelope of the polygon (all its edges should lie on the edges of the polygon boundary))

Here is some code below to generate a random polygon

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

typedef struct Point2D
{
    double x;
    double y;
} Point2D;

int main()
{
    srand( time( NULL ) );
    const unsigned int POLY_SIZE = 20;
    //Polygon, two consecutive array indices form a line in the boundary of the
    // polygon. last point assumed connects to the first point
    Point2D *pPoly = malloc( sizeof( Point2D ) * POLY_SIZE );

    //Generate the Polygon
    const double polyRangeMax = 100.0;
    const double polyRangeMin = -100.0;
    const double polyPointRange = (polyRangeMax - polyRangeMin); 
    const double polyRangeScale = polyPointRange / RAND_MAX;
    for( unsigned int polyPointIdx = 0; polyPointIdx < POLY_SIZE; ++polyPointIdx )
    {
        pPoly[polyPointIdx].x = polyRangeMin + ( rand() * polyRangeScale );
        pPoly[polyPointIdx].y = polyRangeMin + ( rand() * polyRangeScale );
    }

    for( unsigned int polyPointIdx = 0; polyPointIdx < POLY_SIZE; ++polyPointIdx )
    {
        printf("<%f,%f>",pPoly[polyPointIdx].x,pPoly[polyPointIdx].y);
        if( polyPointIdx == POLY_SIZE-1 )
        {
            printf("\n");
        }
        else
        {
            printf(" ");
        }
    }


    //Calculate Outer boundary of the polygon (it itself is a polygon), takes at most same amount of memory as pPoly above;
    unsigned int boundaryPolyPointCount = 0; //TODO calculate this alongside the poly
    Point2D *pOuterBoundaryPoly = malloc( sizeof( Point2D ) * POLY_SIZE );


    //TODO



    for( unsigned int polyPointIdx = 0; polyPointIdx < boundaryPolyPointCount; ++polyPointIdx )
    {
        printf("<%f,%f>",pOuterBoundaryPoly[polyPointIdx].x,pOuterBoundaryPoly[polyPointIdx].y);
        if( polyPointIdx == POLY_SIZE-1 )
        {
            printf("\n");
        }
        else
        {
            printf(" ");
        }
    }

    free( pOuterBoundaryPoly );
    free( pPoly );
    return 0;
}

Every answer I see on Stack Overflow refers to JavaScript or Python libraries and never the actual algorithms themselves. If there is already an answer, let me know and I'll close this question, but I couldn't find anything.

Answer 1

Bentley-Ottmann, and a sketch of an algorithm using it, were mentioned in the comments. That's the most efficient way of calculating the boundary. But if your polygon has less than, say, 1000 vertices and self-intersections, you can use a simpler algorithm and still get realtime performance.

First, identify an extreme vertex of the polygon (say, the one with the lowest X coordinate, with ties broken by the Y coordinate). Of the two edges coming out from that vertex, pick the one which is CW of the other (or equivalently has a lower y/x slope). That's your starting point and your starting edge, which are known to be in the boundary.

Now, loop over all the other edges in the polygon. For each one, check whether it intersects the current edge at a point "forward" (not necessarily >x or >y) of the starting point. Of the edges that intersect, track which one intersects the current edge first (that is, the intersection point is closest to the start point of the edge). Ignore edges which intersect at the start point. Remember, there's always at least one valid intersecting edge, the one which follows the current one in the original polygon. Whichever edge is closest, is the next segment in the boundary, and the intersection is the next starting point of that segment.

Eventually, the next starting point will be the original vertex. Then you know you're done.

This algorithm needs to be tweaked a bit if you expect to have points on line segments, or multiple edges intersecting in the same place. To handle these, only accept an intersection between the middle of one edge and the end of another edge if that other edge continues to the left of the current edge, and if multiple intersections are "nearest" then pick the one that results in the greatest right (clockwise) turn.