Union Polygon Algorithm

Question

This is my attempt at creating a Union Algorithm for an arbitrary set of any number of arbitrary simple polygons. It works for both convex and concave polygons.

public static List<List<Vector2>> UnionPolygon(List<List<Vector2>> uvPolygons)
{
    List<List<Vector2>> newPolygons = new List<List<Vector2>>();

    #region Fix Winding Order
    for (int i = 0; i < uvPolygons.Count; i++)
    {
        float sum = new float();
        for (int j = 0; j < uvPolygons[i].Count; j++)
        {
            Vector3 a = uvPolygons[i][j];
            Vector3 b = uvPolygons[i][(j + 1) % uvPolygons[i].Count];
            sum += (b.x - a.x) * (b.y + a.y);
        }
        if (Mathf.Sign(sum) < 0)
            uvPolygons[i].Reverse();
    }
    #endregion

    #region Find Outside Point
    Vector2 pointOutside = new Vector2();
    for (int i = 0; i < uvPolygons.Count; i++)
    {
        for (int j = 0; j < uvPolygons[i].Count; j++)
        {
            if (uvPolygons[i][j].x > pointOutside.x)
            {
                pointOutside = uvPolygons[i][j];
            }
        }
    }
    pointOutside += new Vector2(1f, 0);
    #endregion

    while (uvPolygons.Count > 0)
    {
        List<Vector2> newPolygon = new List<Vector2>();

        int P = 0; //currentPolygon
        int I = 0; //currentIndex

        #region Setup First Point
        //Make sure starting point is not in any polygons
        for (int p = 0; p < uvPolygons.Count; p++)
        {
            if (p != P)
            {
                if (!CheckIndex(I, 0, uvPolygons[0].Count))
                    break;
                if (IsPointInPolygon(uvPolygons[0][I], uvPolygons[p], pointOutside))
                {
                    p = 0;
                    I++;
                    continue;
                }
            }
        }
        #endregion

        HashSet<int> hashedPolygons = new HashSet<int>();

        while (true)
        {
            hashedPolygons.Add(P);

            if (newPolygon.Contains(uvPolygons[P][I % uvPolygons[P].Count]))
                break;
               
            Vector2 a = uvPolygons[P][I % uvPolygons[P].Count];
            Vector2 b = uvPolygons[P][(I + 1) % uvPolygons[P].Count];

            newPolygon.Add(a);

            bool intersected = false;
            Vector3 intersection = new Vector3();
            Vector3 closestIntersection = new Vector3();
            float closestDistance = new float();

            int tp = new int();
            int ti = new int();
            for (int p = 0; p < uvPolygons.Count; p++)
            {
                if (p != P)
                {
                    for (int i = 0; i < uvPolygons[p].Count; i++)
                    {
                        Vector2 x = uvPolygons[p][i];
                        Vector2 y = uvPolygons[p][(i + 1) % uvPolygons[p].Count];

                        if (AreLinesIntersecting(a, b, x, y, false))
                        {
                            if (LineLineIntersection(a, (b - a).normalized, x, (y - x).normalized, out intersection))
                            {
                                if (newPolygon.Contains(intersection))
                                    continue;

                                if (!intersected)
                                {
                                    closestIntersection = intersection;
                                    closestDistance = Vector3.Distance(a, intersection);

                                    tp = p;
                                    ti = i;
                                    intersected = true;
                                }
                                else if (Vector3.Distance(a, intersection) < closestDistance)
                                {
                                    closestIntersection = intersection;
                                    closestDistance = Vector3.Distance(a, intersection);
                                    tp = p;
                                    ti = i;
                                }
                            }
                        }
                    }
                }
            }
            if (intersected)
            {
                newPolygon.Add(closestIntersection);
                P = tp;
                I = ti;          
            }
            I++;
        }

        if (hashedPolygons.Count == 1)
        {
            newPolygons.Add(new List<Vector2>(newPolygon));
            uvPolygons.RemoveAt(0);
            continue;
        }

        for (int i = uvPolygons.Count - 1; i >= 0; i--)
            if (hashedPolygons.Contains(i))
                uvPolygons.RemoveAt(i);
        if (uvPolygons.Count > 0)
            uvPolygons.Add(new List<Vector2>(newPolygon));
    }
    return newPolygons;
}

public static bool CheckIndex(int i, int min, int max)
{
    if (i >= min && i < max)
        return true;
    else
        return false;
}
public static bool AreLinesIntersecting(Vector2 p1x, Vector2 p1y, Vector2 p2x, Vector2 p2y, bool shouldIncludeEndPoints)
{
    bool isIntersecting = false;

    float denominator = (p2y.y - p2x.y) * (p1y.x - p1x.x) - (p2y.x - p2x.x) * (p1y.y - p1x.y);

    //Make sure the denominator is > 0, if not the lines are parallel
    if (denominator != 0f)
    {
        float u_a = ((p2y.x - p2x.x) * (p1x.y - p2x.y) - (p2y.y - p2x.y) * (p1x.x - p2x.x)) / denominator;
        float u_b = ((p1y.x - p1x.x) * (p1x.y - p2x.y) - (p1y.y - p1x.y) * (p1x.x - p2x.x)) / denominator;

        //Are the line segments intersecting if the end points are the same
        if (shouldIncludeEndPoints)
        {
            //Is intersecting if u_a and u_b are between 0 and 1 or exactly 0 or 1
            if (u_a >= 0f && u_a <= 1f && u_b >= 0f && u_b <= 1f)
            {
                isIntersecting = true;
            }
        }
        else
        {
            //Is intersecting if u_a and u_b are between 0 and 1
            if (u_a > 0f && u_a < 1f && u_b > 0f && u_b < 1f)
            {
                isIntersecting = true;
            }
        }

    }

    return isIntersecting;
}
public static bool LineLineIntersection(Vector3 linePoint1, Vector3 lineVec1, Vector3 linePoint2, Vector3 lineVec2, out Vector3 intersection)
{
    Vector3 lineVec3 = linePoint2 - linePoint1;
    Vector3 crossVec1and2 = Vector3.Cross(lineVec1, lineVec2);
    Vector3 crossVec3and2 = Vector3.Cross(lineVec3, lineVec2);

    float planarFactor = Vector3.Dot(lineVec3, crossVec1and2);

    //is coplanar, and not parallel
    if (Mathf.Abs(planarFactor) < 0.0001f
            && crossVec1and2.sqrMagnitude > 0.0001f)
    {
        float s = Vector3.Dot(crossVec3and2, crossVec1and2)
                / crossVec1and2.sqrMagnitude;
        intersection = linePoint1 + (lineVec1 * s);
        return true;
    }
    else
    {
        intersection = Vector3.zero;
        return false;
    }
}
public static bool IsPointInPolygon(Vector2 point, List<Vector2> polygonPoints, Vector2 pointOutside)
{
    int numIntersections = 0;

    for (int j = 0; j < polygonPoints.Count; j++)
    {

        Vector2 uv1 = polygonPoints[j];
        Vector2 uv2 = polygonPoints[(j + 1) % polygonPoints.Count];

        if (AreLinesIntersecting(point, pointOutside, uv1, uv2, true))
            numIntersections++;

    }

    if (numIntersections == 0 || numIntersections % 2 == 0)
        return false;
    else
        return true;
}

Briefly on how this algorithm works is it follows a polygon through its winding order until it detects an intersection. Upon intersection, it jumps to the intersected polygon. Repeat until it has made its way around to the beginning.

Setup requirements for this to work correctly:

• The Starting Index Must Be Outside All Other Polygons
• All Polygons Must Share Winding Order
• All Polygons Must Be Simple (i.e. not open or self intersecting)

Finally, if two polygons share an intersecting edge, the algorithm will "jump over" the second polygon. (as seen in figure 1) So it keeps track of the number of starting polygons and will repeat until all are accounted for.

It has a superficial resemblance to the Weiler Atherton algorithm

Answer 1

Regions

Well, regions can help group some things together, but most of the times having regions indicates the code in question has some problems.
Regions grouping methods inside a class normally indicates that the method has too many responsibilities which is a design flaw.
Regions inside a method is just a no go, because they indicate that the method itself is doing too much and tends to be a god-method.
#region Fix Winding Order just screams for a method which is doing just what the regions-text says, which applies for #region Find Outside Point and #region Setup First Point as well.

---

float sum = new float();  

Why don't you use float sum = 0f; ? that would be consistent. You have int P = 0; as well instead of int P = new int(); Beeing consistent in the way the code is written makes it much easier to read and understand.

---

public static bool CheckIndex(int i, int min, int max)
{
    if (i >= min && i < max)
        return true;
    else
        return false;
}  

A construct, like if a condition is true return true otherwise return false, can and should be shortened to just return the codition like

public static bool CheckIndex(int i, int min, int max)
{
    return i >= min && i < max;
}  

In addition, omitting braces althought they might be optional can lead to hidden and therefor hard to find bugs.

---

If AreLinesIntersecting() is behaving as it should, it could be greatly enhanced by

• returning early for denominator ==0f, well at the second glance that only could be true. If it is true then the comment //Make sure the denominator is > 0, if not the lines are parallel is a lie.
• combining the if and the else` branche into one statement

like so (taking into account that the code returns the same results as the former method):

public static bool AreLinesIntersecting(Vector2 p1x, Vector2 p1y, Vector2 p2x, Vector2 p2y, bool shouldIncludeEndPoints)
{
    float denominator = (p2y.y - p2x.y) * (p1y.x - p1x.x) - (p2y.x - p2x.x) * (p1y.y - p1x.y);

    //Make sure the denominator is > 0, if not the lines are parallel
    if (denominator == 0f) { return false; }

    float u_a = ((p2y.x - p2x.x) * (p1x.y - p2x.y) - (p2y.y - p2x.y) * (p1x.x - p2x.x)) / denominator;
    float u_b = ((p1y.x - p1x.x) * (p1x.y - p2x.y) - (p1y.y - p1x.y) * (p1x.x - p2x.x)) / denominator;

    return u_a > 0f && u_a < 1f && u_b > 0f && u_b < 1f;
}

Answer 2

Style-wise, your code could use some improvements.

First, you use #regions as comments of sort, which is not their intended use. You could just use comments instead, or do without them completely.

You could also use more descriptive variable names in some instances. For example, in the following lines:

        int P = 0; //currentPolygon
        int I = 0; //currentIndex

why not use currentPolygon and currentIndex as variable names? It would make reading the code much easier down the line.

Some other variable names are downright confusing. Consider the signature of this function:

public static bool AreLinesIntersecting(Vector2 p1x, Vector2 p1y, Vector2 p2x, Vector2 p2y, bool shouldIncludeEndPoints)

To me, it would make sense for p1x and p1y to be names for the x- and y-coordinates of point 1, but this doesn't match the Vector2 data type. In fact, it looks like they are the start and end points of line 1, so line1Start and line1End would be much more readable names.

You could also improve documentations by adding documentation comments to your methods. This would be quite valuable since you have some conditions on the polygons you pass to the function (simple polygons all with the same winding order) without validation before executing the algorithm.

---

Logic-wise, my main issue with your code is your mixed use of Vector2 and Vector3. Since your code solves a purely planar problem, Vector3 are no use here. In the first section fixing winding order, it looks like your Vector3s are a typo, which went unnoticed because Unity allows implicit casting between Vector2 and Vector3, but you don't use the third dimension at all here.

You also use Vector3s later for the LineLineIntersection method. Again, since your problem is planar, there should be no use for the third dimension. It looks like you reused this algorithm from a 3D application: it starts with a coplanarity check, when in fact, since you're working in 2D, you are guaranteed to work with coplanar vectors.

Cross products of 2D vectors are a float, you can then work with that instead of relying on the Vector3.Cross method, which is wasteful and confusing here.

And in fact, your AreLinesIntersecting method almost goes all the way in finding the intersection between the line, I believe all that is left to do is to assign p1x + u_a * (p2x - p1x) to an out parameter and you can get rid of the LineLineIntersection completely.