Deltares / numba_celltree

Celltree data structure for searching for points, lines, boxes, and cells (convex polygons) in a two dimensional unstructured mesh.
https://deltares.github.io/numba_celltree/
MIT License
11 stars 3 forks source link

Faster point in convex polygon check #8

Open Huite opened 1 year ago

Huite commented 1 year ago

Checking for separating planes instead of relying on Jordan ray theorem, something like this:

@nb.njit(inline="always")
def point_in_convex_polygon(p: Point, poly: FloatArray) -> bool:
    """Half-plane check."""
    positive = False 
    negative = False 
    length = len(poly)
    v0 = as_point(poly[-1])
    for i in range (length):
        v1 = as_point(poly[i])
        U = to_vector(p, v0)
        W = to_vector(v0, v1)
        if cross_product(U, W) > 0:
            positive = True
        else:
            negative = True

        if positive and negative:
            return False

        v0 = v1

    return True

@nb.njit(inline="always")
def point_in_convex_polygon_or_on_edge(p: Point, poly: FloatArray) -> bool:
    """Half-plane check."""
    positive = False 
    negative = False 
    length = len(poly)
    v0 = as_point(poly[-1])
    for i in range(length):
        v1 = as_point(poly[i])
        U = to_vector(p, v0)
        V = to_vector(p, v1)
        W = to_vector(v0, v1)

        twice_area = abs(cross_product(U, V))
        if twice_area < TOLERANCE_ON_EDGE:
            # Project the point onto W.
            W = to_vector(v0, v1)
            if W.x != 0.0:
                t = (p.x - v0.x) / W.x
            elif W.y != 0.0:
                t = (p.y - v0.y) / W.y
            else:
                # The vector has a length of zero. Skip.
                continue
            if 0 <= t <= 1:
                # Now it's on the edge.
                return True

        if cross_product(U, W) > 0:
            positive = True
        else:
            negative = True

        if positive and negative:
            return False

        v0 = v1
        U = V

    return True

This seems to be approximately two times faster than the current check. However, it's obviously less robust (cannot deal with concave polygons). Annoyingly (and expectedly), it reports different results for edge cases such as a point located on a vertex. Given only moderate speed-up and less robustness, it's probably best to stick to the current implementation (also given status quo bias...).