pythonmathgeometryunionpolygon

unknown bug in finding points inside a polygon and calculating the union of several polygons


I hope you are well. I am writing a code that finds the union of several polygons. I know that libraries do this but I need to do it from scratch. When I ran my code I got this output that was obviously wrong because of the left-bottom part.

enter image description here

I really don't know what could be the reason for this problem.

When I searched in my code and wrote some test codes I found another bug that I don't know if is related to the bug of the union part or not but anyway, this is a part of that union function. the problem is in the contains_point function. It works great when a point is inside or outside of the polygon.

enter image description here

but not when the point is on the edges of the polygon.

enter image description here

I first thought maybe it was because of floating point numbers that some points stand inside and rest outside but after, I realized that there is no red point on the top and left-top segments of the polygon. so it couldn't be for floating point behavior but I don't know the reason for this bug.

please help me to find it. is the problem with my algorithm or implementation or something else?

union part of the code

def find_intersect(a1: Point, a2: Point, b1: Point, b2: Point) -> dict | None:
    t_top = (b2.x - b1.x) * (a1.y - b1.y) - (b2.y - b1.y) * (a1.x - b1.x)
    u_top = (b1.y - a1.y) * (a1.x - a2.x) - (b1.x - a1.x) * (a1.y - a2.y)
    bottom = (b2.y - b1.y) * (a2.x - a1.x) - (b2.x - b1.x) * (a2.y - a1.y)
    if bottom != 0:
        t = t_top / bottom
        u = u_top / bottom
        if 0 <= t <= 1 and 0 <= u <= 1:
            return {
                "x": lerp(a1.x, a2.x, t),
                "y": lerp(a1.y, a2.y, t),
                "offset": t,
            }
    return None


class Polygon:
    def __init__(self, points: list) -> None:
        self.points = points
        self.segments = []
        for i, _ in enumerate(self.points):
            self.segments.append(
                Segment(self.points[i], self.points[(i + 1) % len(self.points)])
            )

    @staticmethod
    def union(polygons: list) -> list:
        i = 0
        while i < len(polygons) - 1:
            j = i + 1
            while j < len(polygons):
                polygons[i].break_segments(polygons[j])
                j += 1
            i += 1
        kept_segments = []
        for i, polygon_i in enumerate(polygons):
            for segment in polygon_i.segments:
                keep = True
                for j, polygon_j in enumerate(polygons):
                    if i != j:
                        if polygon_j.contains_segment(segment):
                            keep = False
                            break
                if keep:
                    kept_segments.append(segment)
        return kept_segments

    def break_segments(self, other: Self) -> None:
        segments_1 = self.segments
        segments_2 = other.segments
        for i, _ in enumerate(segments_1):
            for j, _ in enumerate(segments_2):
                intersection = find_intersect(
                    segments_1[i].start,
                    segments_1[i].end,
                    segments_2[j].start,
                    segments_2[j].end,
                )
                if (
                    intersection
                    and intersection["offset"] != 0
                    and intersection["offset"] != 1
                ):
                    point = Point(intersection["x"], intersection["y"])
                    temporary = segments_1[i].end
                    segments_1[i].end = point
                    segments_1.insert(i + 1, Segment(point, temporary))
                    temporary = segments_2[j].end
                    segments_2[j].end = point
                    segments_2.insert(j + 1, Segment(point, temporary))

    def contains_segment(self, segment: Segment) -> bool:
        return self.contains_point(segment.midpoint())

contains_point function of the polygon class

def contains_point(self, point: Point) -> bool:
    px = point.x
    py = point.y
    intersections = 0
    for segment in self.segments:
        x1 = segment.start.x
        y1 = segment.start.y
        x2 = segment.end.x
        y2 = segment.end.y
        if (py < y1) != (py < y2):
            if ((px < x1) != (px < x2)) and (
                (x2 - x1) * (py - y1) == (px - x1) * (y2 - y1)
            ):
                return True
            if px < (x2 - x1) * (py - y1) / (y2 - y1) + x1:
                intersections += 1
    if intersections % 2 == 1:
        return True
    return False

I want to calculate the union of polygons (picture 1) but correctly.


Solution

  • I think the issue with calculating the union lies within the contains_point function, and your problem in that function is specifically with comparing floating point numbers. As you said, some points stand on the right side of the line and others are on the left and barely or never you can find a point exactly on the line. Therefore, the first if statement checking for the same slope between the line and the new line from the start to the new point is unnecessary. However, for the second if statement, you should include a tiny threshold to account for points that are almost on the lines.

    You may need something like this code:

    def contains_point(self, point: Point, threshold: float = 0.01) -> bool:
        px = point.x
        py = point.y
        intersections = 0
        for segment in self.segments:
            x1 = segment.start.x
            y1 = segment.start.y
            x2 = segment.end.x
            y2 = segment.end.y
            if (py < y1) != (py < y2):
                if px - threshold < (x2 - x1) * (py - y1) / (y2 - y1) + x1:
                    intersections += 1
        if intersections % 2 == 1:
            return True
        return False
    

    The code has a problem.

    It is theoretically possible that a point was on the right of one of the right lines (or on the right of one of the left lines) of the polygon but so close to it that when you subtract the threshold from the x coordinate of the point it goes inside (or outside) the polygon but it happens so rarely. I don't have any idea how to solve this new problem.