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

intersect_edges does not appear to be entirely robust #26

Closed Huite closed 1 week ago

Huite commented 1 week ago

This is fine:

import numpy as np
import numba_celltree

triangles = np.fromfile("c:/tmp/triangles.bin", dtype=np.int64).reshape((-1, 3))
vertices = np.fromfile("c:/tmp/vertices.bin", dtype=np.float64).reshape((-1, 2))

tree = numba_celltree.CellTree2d(vertices, triangles, -1)

# %%

xmin = vertices[:, 0].min()
ymin = vertices[:, 1].min()
xmax = vertices[:, 0].max()
ymax = vertices[:, 1].max()
x = np.linspace(xmin, xmax, 1000)
y = np.linspace(ymin, ymax, 1000)
yy, xx = [a.ravel() for a in np.meshgrid(y, x, indexing="ij")]
points = np.column_stack((xx, yy))
edges = points.reshape((-1, 2, 2))

# %%

tree.intersect_edges(edges)
# %%

But testing a grid against its own edges either enters a infinite loop or something else:

edges = numba_celltree.demo.edges(triangles, -1)
edge_coords = vertices[edges]
tree.intersect_edges(edge_coords)

Needs testing with a smaller example.

Huite commented 1 week ago

This is due to dx, dy not getting updated: https://github.com/Deltares/numba_celltree/blob/374a7a4a7606952d2487b0d96ec5c6603cfd286f/numba_celltree/algorithms/cohen_sutherland.py#L74

It can be triggered by this:

Image

Where b is equal to (xmin, ymax).

Liang-Barsky is also not consistent here, since it returns True for intersection although the resulting edge is a single point.

Huite commented 1 week ago

Fixed by f47aa91c49af61657c8d8bb5d8aa60556e251518