"""
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.
E.g. something like: