# Intersection of level sets

**URL:** <https://discourse.julialang.org/t/intersection-of-level-sets/114913>\
**Category:** Performance\
**Tags:** question, performance\
**Created:** [May 29, 2024, 12:59pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913 "2024-05-29T12:59:03Z")\
**Posts on this page:** 1\
**Showing post:** 4

<div class="post-metadata">

**Author:** ![cmarcotte](https://avatars.discourse-cdn.com/v4/letter/c/a3d4f5/32.png) [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)\
**Post date:** [May 30, 2024, 12:35pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/4 "2024-05-30T12:35:41Z")

</div>

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](https://gist.github.com/cmarcotte/5eedff00c8a9e21eeacfcb3d9ec3d2a2) 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 `Segment`s – 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](https://discourse.julialang.org/t/finding-crossing-points-of-two-vectors-of-vertices/15786/7) thread, which appears extremely relevant.
  - `GeometryTypes` is archived in favor of `GeometryBasics`, and the documentation is unclear if or how intersection tests for `LineString`s 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`](https://github.com/StellaOrg/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.

---

_[View the full topic](https://discourse.julialang.org/t/intersection-of-level-sets/114913)._