Intersection of level sets

I have managed to speed up the calculation significantly (~50x) by switching to a mutable struct for the Rectangle and updating an axis-aligned bounding-box at both the line and segment level, and only proceeding to the linear solve in the case of the latter intersection returning true. The gist is updated (now with plotting!). The timing is now about 625us for the original simple example (two contours from size (512,512) fields, 32ms for random input of size (64,64), and 9s for random input of size (256,256) yielding 24107 intersections. The previously 4h calculation is now down to 6m.

I can’t quite parse your meaning, @juliohm – do you mean an SDF-based approach? Or should I be able to call something like mutualIntersections = intersect(cs0, cs1) and get nice dispatch to something sub-quadratic?

There are three issues, as I understand it:

  1. Determining whether two segments intersect based on their coordinates, which is simple enough (I think Meshes.jl does this really well on Segments – certainly more robustly than the 4x4 linear solve); and
  2. Improving the intersection testing efficiency for two sets of line segments so it’s not quadratic complexity. The new gist implementation is still quadratic but the prefactor is now tractably small.
  3. Improving the intersection testing efficiency for multiple sets of chains, so it’s not similarly quadratic. This new approach is still quadratic, but should have early exits for my use-case.

Some further thoughts for future people who stumble on this thread:

  • I also came across this thread, which appears extremely relevant.
    • GeometryTypes is archived in favor of GeometryBasics, and the documentation is unclear if or how intersection tests for LineStrings are done – from a quick look, it doesn’t seem to be implemented.
    • CollisionDetection appears to solve precisely this problem as @Jeff_Emanuel suggests (octrees), though there is some issue if the input data is e.g. Float32 and it appears to be unmaintained.
  • The Bounding Volume Hierarchy structure seems like the correct approach, for which there is a very nice package ImplicitBVH.jl; it does not appear to include support for 2D elements (a BRect{T} to their BBox{T}), but perhaps could support the GeometryBasics types as leafs, or might be useful for embedding the 2D contours into 3D, and it explicitly handles dual-bvh intersections efficiently.
  • Meshes.jl has implemented all the edge-cases for Segment intersections, and it appears Meshes.jl has a method for intersection(::chain, ::chain) but it’s unclear whether this is building a tree for efficient testing, or falling back to a quadratic method, as here.