Area of self-intersecting polygon

Question

Calculating the area of a simple irregular polygon is trivial. However, consider the self-intersecting polygon ABCDEF shown at left below:

                   

If we use the linked-to formula above traversing the points in polygon order, we get an area of 0. (The 'clockwise' area cancels out the 'counter-clockwise' area.)

However, if we sort the points radially around a center and calculate the area, we get the incorrect area of the polygon ABEDCF at right above.

How can I best find the visible area of a self-intersecting polygon? (If the answer requires creating phantom points for each intersection, please provide details for how to best find the intersections and how then to traverse them in correct order.)

This question arose when investigating edge cases for my solution to this question.

Defining the Area

I define the 'area' as the amount of pixels visible when filling the polygon using either the "nonzero" or "evenodd" rules. I will accept an answer for either of these, though both would be better. Note that I explicitly do not define the area for self-overlapping to count the overlapping area twice.

Answer 1

You can try Bentley–Ottmann with the following pseudo code from this page

The Bentley-Ottmann Algorithm

The input for the Bentley-Ottmann algorithm is a collection OMEGA={Li} of line segments Li, and its output will be a set LAMBDA={Ij} of intersection points. This algorithm is referred to as a "sweep line algorithm" because it's operation can be visualized as having another line, a "sweep line" SL, sweeping over the collection OMEGA and collecting information as it passes over the individual segments Li. The information collected for each position of SL is basically an ordered list of all segments in OMEGA that are currently being intersected by SL. The data structure maintaining this information is often also called the "sweep line". This class structure also detects and outputs intersections as it discovers them. The process by which it discovers intersections is the heart of the algorithm and its efficiency.

To implement the sweep logic, we must first linearly order the segments of OMEGA to determine the sequence in which SL encounters them. That is, we need to order the endpoints {Ei0,Ei1}i=1,n of all the segments Li so we can detect when SL starts and stops intersecting each segment of OMEGA. Traditionally, the endpoints are ordered by increasing x first and then increasing y-coordinate values, but any linear order will do (some authors prefer decreasing y first and then increasing x). With the traditional ordering, the sweep line is vertical and moves from left to right as it encounters each segment, as shown in the diagram:

Pic-sweepline

At any point in the algorithm, the sweep line SL intersects only those segments with one endpoint to the left of (or on) it and the other endpoint to the right of it. The SL data structure keeps a dynamic list of these segments by: (1) adding a segment when its leftmost endpoint is encountered, and (2) deleting a segment when its rightmost endpoint is encountered. Further, the SL orders the list of segments with an "above-below" relation. So, to add or delete a segment, its position in the list must be determined, which can be done by a worst-case O(log-n) binary search of the current segments in the list. In addition, besides adding or deleting segments, there is another event that changes the list structure; namely, whenever two segments intersect, then their positions in the ordered list must be swapped. Given the two segments, which must be neighbors in the list, this swap is an O(log-n) operation.

To organize all this, the algorithm maintains an ordered "event queue" EQ whose elements cause a change in the SL segment list. Initially, EQ is set to the sweep-ordered list of all segment endpoints. But as intersections between segments are found, then they are also added to EQ in the same sweep-order as used for the endpoints One must test, though, to avoid inserting duplicate intersections onto the event queue. The example in the above diagram shows how this can happen. At event 2, segments S1 and S2 cause intersection I12 to be computed and put on the queue. Then, at event 3, segment S3 comes between and separates S1 and S2. Next, at event 4, S1 and S3 swap places on the sweep line, and S1 is brought next to S2 again causing I12 to be computed again. But, there can only be one event for each intersection, and I12 cannot be put on the queue twice. So, when an intersection is being put on the queue, we must find its potential x-sorted location in the queue, and check that it is not already there. Since there is at most one intersect point for any two segments, labeling an intersection with identifiers for the segments is sufficient to uniquely identify it. As a result of all this, the maximum size of the event queue = 2n+k.le.2n+n2, and any insertion or deletion can be done with a O(log(2n+n2)) = O(log-n) binary search.

But, what does all this have to do with efficiently finding the complete set of segment intersections? Well, as segments are sequentially added to the SL segment list, their possible intersections with other eligible segments are determined. When a valid intersection is found, then it is inserted into the event queue. Further, when an intersection-event on EQ is processed during the sweep, then it causes a re-ordering of the SL list, and the intersection is also added to the output list LAMBDA. In the end, when all events have been processed, LAMBDA will contain the complete ordered set of all intersections.

However, there is one critical detail, the heart of the algorithm, that we still need to describe; namely, how does one compute a valid intersection? Clearly, two segments can only intersect if they occur simultaneously on the sweep-line at some time. But this by itself is not enough to make the algorithm efficient. The important observation is that two intersecting segments must be immediate above-below neighbors on the sweep-line. Thus, there are only a few restricted cases for which possible intersections need to be computed:

When a segment is added to the SL list, determine if it intersects with its above and below neighbors.

When a segment is deleted from the SL list, its previous above and below neighbors are brought together as new neighbors. So, their possible intersection needs to be determined.

At an intersection event, two segments switch positions in the SL list, and their intersection with their new neighbors (one for each) must be determined. This means that for the processing of any one event (endpoint or intersection) of EQ, there are at most two intersection determinations that need to be made.

One detail remains, namely the time needed to add, find, swap, and remove segments from the SL structure. To do this, the SL can be implemented as a balanced binary tree (such as an AVL, a 2-3, or a red-black tree) which guarantees that these operations will take at most O(log-n) time since n is the maximum size of the SL list. Thus, each of the (2n+k) events has at worst O(log-n) processing to do. Adding up the initial sort and the event processing, the efficiency of the algorithm is: O(nlog-n)+O((2n+k)log-n)=O((n+k)log-n).

Pseudo-Code: Bentley-Ottmann Algorithm

Putting all of this together, the top-level logic for an implementation of the Bentley-Ottmann algorithm is given by the following pseudo-code:

    Initialize event queue EQ = all segment endpoints;
    Sort EQ by increasing x and y;
    Initialize sweep line SL to be empty;
    Initialize output intersection list IL to be empty;

    While (EQ is nonempty) {
        Let E = the next event from EQ;
        If (E is a left endpoint) {
            Let segE = E's segment;
            Add segE to SL;
            Let segA = the segment Above segE in SL;
            Let segB = the segment Below segE in SL;
            If (I = Intersect( segE with segA) exists) 
                Insert I into EQ;
            If (I = Intersect( segE with segB) exists) 
                Insert I into EQ;
        }
        Else If (E is a right endpoint) {
            Let segE = E's segment;
            Let segA = the segment Above segE in SL;
            Let segB = the segment Below segE in SL;
            Delete segE from SL;
            If (I = Intersect( segA with segB) exists) 
                If (I is not in EQ already) 
                    Insert I into EQ;
        }
        Else {  // E is an intersection event
            Add E’s intersect point to the output list IL;
            Let segE1 above segE2 be E's intersecting segments in SL;
            Swap their positions so that segE2 is now above segE1;
            Let segA = the segment above segE2 in SL;
            Let segB = the segment below segE1 in SL;
            If (I = Intersect(segE2 with segA) exists)
                If (I is not in EQ already) 
                    Insert I into EQ;
            If (I = Intersect(segE1 with segB) exists)
                If (I is not in EQ already) 
                    Insert I into EQ;
        }
        remove E from EQ;
    }
    return IL;
}

This routine outputs the complete ordered list of all intersection points.

Answer 2

This is from top of my mind, I'd assume there is no hole in your polygon, with hole it will be more complicated and you should first remove holes from your poly:

1. First find convex hull of your polygon, for this you need to find convex hull of your polygon vertices. And compute convex hull area.

2. After that, find all intersection of your polygon.

3. You should subtract extra polygons which doesn't belong to your original polygon from convex hull to find your polygon area, name them badpoly. badpolys always have at least one border on convex hull, such that this border does not belong to your original polygon, name them badborder, by iterating over convex hull you can find all badborders, but for finding other borders of badpoly, next connected border to given badborder which has smallest angle relative to badborder is one of a borders of your badpoly, you can continue this to find all borders of your badpoly and then calulate its area, also by repeating this way you can calculate area of all badpolys.