Write some examples / documentation on edge cases, tolerances, etc.

[Huite]

E.g. something like:

"""
Edge cases and numeric robustness
=================================

Edge cases and numeric robustness are an especially difficult problems in
computational geometry. This examples illustrates a few aspects.
"""
# %%
import matplotlib.pyplot as plt
import numba_celltree
from numba_celltree import demo
from numba_celltree.constants import Point

import numpy as np

import numba_celltree.geometry_utils

# %%
# We'll set up a basic example of a mesh with three quad faces.

nodes = np.array(
    [ 
        [0.0, 0.0], # 0
        [2.0, 0.0], # 1
        [2.0, 2.0], # 2
        [0.0, 2.0], # 3
        [4.0, 0.0], # 4
        [4.0, 4.0], # 5
        [0.0, 4.0], # 6
    ]
)
faces = np.array(
    [
        [0, 1, 2, 3],
        [1, 4, 5, 2],
        [3, 2, 5, 6],
    ]
)

fig, ax = plt.subplots()
demo.plot_edges(*nodes.T, demo.edges(faces, -1), ax=ax)

# %%
# We can the first single face and test a few points.

polygon = nodes[faces[0]]

points = np.array(
    [
        [-0.5, -0.5],
        [0.0, 0.0],
        [0.5, 0.0],
        [1.0, 1.0],
        [1.0, 2.0],
        [2.0, 1.0],
        [2.0, 2.0],
    ]
)
inside = [
    numba_celltree.geometry_utils.point_in_polygon(Point(*xy), polygon) for xy in points
]

fig, ax = plt.subplots()
demo.plot_edges(*nodes.T, demo.edges(faces, -1), ax=ax)
ax.scatter(*points.T, c=inside)

# %%
# All but a single point (1.0, 1.0) are edges cases. Some are identified as
# inside the polygon, some outside. For a single polygon, this is a reasonable
# answer as the edge is neither inside or out. However, for a mesh, the
# situation is different: a point at (2.0, 2.0) is not part of any face, but it
# is certainly within the mesh.
#
# This second example demonstrates edge cases and floating point roundoff:

nodes2 = np.array([
    [0.0, 0.0],
    [3.0, 0.0],
    [1.0, 1.0],
    [0.0, 2.0],
    [3.0, 2.0],
]
)
faces2 = np.array([
    [0, 1, 2],
    [0, 2, 3],
    [2, 4, 3],
])
points2 = np.array([
    [0.0, 0.0],
    [0.05, 0.05],
    [0.15, 0.15],
    [0.35, 0.35],
    [0.55, 0.55],
    [0.75, 0.75],
    [0.95, 0.95],
    [1.0, 1.0],
]) 
triangle = nodes2[faces2[0]]

inside = [
    numba_celltree.geometry_utils.point_in_polygon(Point(*xy), triangle) for xy in points2
]

fig, ax = plt.subplots()
demo.plot_edges(*nodes2.T, demo.edges(faces2, -1), ax=ax)
ax.scatter(*points2.T, c=inside)

# %%
# The first point (which is a vertex) is identified as inside. The next two
# points fall outside; but the other points on the diagonal all fall inside.
# Then, the last point which is also a vertex, falls outside.
#
# An alternative scheme is to include points on the edge as falling inside.
# Unfortunately, due to floating point roundoff, we cannot rely on an exact
# condition to identify points on the edge. We rely on a tolerance instead to
# points that are very close to an edge.

# The celltree uses 64-bit floats to describe coordinates. The precision of
# 64-bit floats is not infinite, and so numbers are effectively rounded. This
# can lead to inconsistencies. A widely used solution is robust predicates (by
# Shewchuk), but this has a cost in terms of code complexity and performance.

inside = [
    numba_celltree.geometry_utils.point_in_polygon_or_on_edge(Point(*xy), polygon) for xy in points
]

fig, ax = plt.subplots()
demo.plot_edges(*nodes.T, demo.edges(faces, -1), ax=ax)
ax.scatter(*points.T, c=inside)

# %%
# A potential problem with this approach is that interior edges are shared by
# two faces and a point may fall into two faces at once.

points = np.array(
    [
        [0.0-1e-10, 0.0-1e-10],
        [2.0+1e-10, 0.0],
        [2.0, 2.0+1e-10],
    ]
)
inside = [
    numba_celltree.geometry_utils.point_in_polygon_or_on_edge(Point(*xy), polygon) for xy in points
]
print(inside)
tree = numba_celltree.CellTree2d(nodes, faces, -1)
face_indices = tree.locate_points(points)

fig, ax = plt.subplots()
demo.plot_edges(*nodes.T, demo.edges(faces, -1), ax=ax)
ax.scatter(*points.T, c=inside)

# %%
# Strictly speaking, all points fall outside of the first face (at a distance
# of at least 1e-10), but they are identified as falling inside as they fall
# within the tolerance.
#
# Numba celltree was created primarily with geopatial applications in mind,
# with meter as the typical unit. In that context, feature coordinates are not
# known at even mm resolutions and such false positives are acceptable.
#
# Please open an issue if this is problematic for your use case.