# 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:** 13
**Page:** 1
### Author:  [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)
#### Post date: [May 29, 2024, 12:59pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/1 "2024-05-29T12:59:03Z")
Hello!
I am computing the intersections of two level sets computed from some input arrays. The idea is that, having retrieved the [contours from Makie](https://discourse.julialang.org/t/extracting-vector-of-contour-lines-from-makie-contour/104866) we can compute the intersections of the curves by iterating over sequential pairs of vertices. The approach I took many years ago was to form a 4x4 linear system to compute the intersection point (and lerp parameters at the intersection between the two linear segments) [like this code](https://uk.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections). This is seemingly robust, but stunningly slow.
A complete gist is available [here](https://gist.github.com/cmarcotte/5eedff00c8a9e21eeacfcb3d9ec3d2a2). The slow portion of the code is the intersection computation when the input field is large (or the length of the contours are long). The relevant function is:
```julia
function computeIntersects!(intersects::Vector{NTuple{2,T}}, cs0, cs1) where T <: Number
A = zeros(T,4,4)
A[[1,2],3] .= T(-1.0)
A[[3,4],4] .= T(-1.0)
B = zeros(T,4)
X = zeros(T,4)
# iterate over levels
for c0 in cs0.contours, c1 in cs1.contours # O(c^2)
# iterate over isolated lines
for l0 in c0.lines
r0 = Rectangle( mapreduce(p->p[1], min, l0.vertices), mapreduce(p->p[1], max, l0.vertices), mapreduce(p->p[2], min, l0.vertices), mapreduce(p->p[2], max, l0.vertices))
for l1 in c1.lines # O(l^2)
r1 = Rectangle( mapreduce(p->p[1], min, l1.vertices), mapreduce(p->p[1], max, l1.vertices), mapreduce(p->p[2], min, l1.vertices), mapreduce(p->p[2], max, l1.vertices))
# iterate over vertices on each line only if the bounding boxes for l0 and l1 intersect:
if rectIntersection(r0, r1) # O(v^2)
# compute intersections
for (p00, p01) in PairIterator(l0.vertices)
A[1,1] = (p01[1]-p00[1])
A[3,1] = (p01[2]-p00[2])
B[1] = -p00[1]
B[3] = -p00[2]
for (q00, q01) in PairIterator(l1.vertices)
A[2,2] = (q01[1]-q00[1])
A[4,2] = (q01[2]-q00[2])
B[2] = -q00[1]
B[4] = -q00[2]
try
X .= A\B
if 0 <= X[1] && X[1] < 1.0 && 0 <= X[2] && X[2] < 1
push!(intersects, (X[3], X[4]))
end
catch
X .= T(NaN)
end
end
end
end
end
end
end
return nothing
end
```
This takes about 30ms on my machine for the (very simple) example in the gist. Random arrays are a much more significant computation (clocking in around 1.5s for an input array which is 1/64 the size). My question is if there is something about the implementation (especially in `computeIntersects!`) which is obviously not a good idea.
---
### Author:  [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)
#### Post date: [May 29, 2024, 2:41pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/2 "2024-05-29T14:41:36Z")
I’d build spatial trees for each contour so you can quickly find which segments might intersect and avoid most intersection tests.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 29, 2024, 5:57pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/3 "2024-05-29T17:57:04Z")
I believe that the whole idea of implicit modeling (level sets) is to avoid such bottlenecks with topological operations such as intersection. If you have level set representations of geometries, you should be able to compute intersections with bit-wise instructions?
---
### Author:  [@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")
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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 30, 2024, 12:45pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/5 "2024-05-30T12:45:53Z")
It is not usual to extract contours to do geometric processing. It is more efficient (usually) to operate on the arrays of values directly (level sets of the form f(x) \< L).
If you really need to extract the segments from the contours, with a explicit discretization of vertices, then search for the “sweep-line algorithm” for intersection between two sets of segments. It has much better computational complexity than the double loop with quadratic complexity implemented above.
This algorithm will be added to Meshes.jl at some point. If you want to contribute, a `Meshes.intersection` method, that is very welcome.
---
### Author:  [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)
#### Post date: [May 30, 2024, 1:31pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/6 "2024-05-30T13:31:05Z")
Thank you very much. I’m not familiar with the direct array approach – is this demonstrated in `Meshes.jl` or elsewhere?
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 30, 2024, 1:53pm UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/7 "2024-05-30T13:53:31Z")
The array approach should be something simple that relies on the f(x) \< Lf and g(x) \< Lg implicit representations. If you want to find all x for which f(x) = L (e.g. contour) and g(x) = L, you can filter all pixels of the two arrays to match the criteria. After that, you can extract the vertices of the filtered pixels. Makes sense?
Perhaps you can gain more insight by reading the implicit modeling literature. There is also ImplicitEquations.jl to experiment with the concepts.
---
### Author:  [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)
#### Post date: [June 3, 2024, 9:38am UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/8 "2024-06-03T09:38:36Z")
As a note for those who stumble on this thread in the future:
_I recommend using the [ImplicitBVH.jl](https://github.com/StellaOrg/ImplicitBVH.jl) package to accelerate the generation of a small set of actual intersection tests_.
Even using a `BSphere{T}` for each segment of each contour for two input `(512,512)` random fields (i.e., overestimating the intersection density significantly in an already extreme test case), determining the potentially intersecting elements is eminently tractable. In an example, this approach reduces the number of intersection tests to be performed from `261587 * 260824 = 68228167688` to `2022690` in under `365.776 ms` (with ~80% GC time).
Even then sweeping the initial segments (since the `BSphere` loses information about the segment points) is fast enough. I’ve updated the gist again, and attach a Haring-esque `(128,128)` example (intersections in green):
---
### Author:  [@WuSiren](https://discourse.julialang.org/u/WuSiren)
#### Post date: [June 5, 2024, 2:49am UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/9 "2024-06-05T02:49:46Z")
Hello! Thank you for your generous sharing! Now I’m encountered with a similar problem where I need to visualize the intersection of the following sets: f\_1(x)\leq c\_1, f\_2(x)\leq c\_2, \cdots, f\_n(x)\leq c\_n where x\in\mathbb{R}^3 (or simply visualize the surface of the intersection). Will your solution apply to this problem?
Thanks!
---
### Author:  [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)
#### Post date: [June 5, 2024, 10:35am UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/10 "2024-06-05T10:35:43Z")
Hello, I am not entirely convinced the approach detailed in the gist above will apply as well for your problem, simply because you have inequalities, rather than equalities (as in my problem). For my use, the discretized level sets are always co-dimension 1 with the original space, and thus their intersections are (with high confidence) co-dimension 2 with the original space.
In your problem, you would expect the intersection of these inequalities to be \Omega\_i \subseteq \mathbb{R}^3, and so \partial\Omega\_i is a collection of 3D points defining a set of surfaces with (approximate) dimension 2.
For the visualization in particular, you may wish to ask a question of the Makie developers in a new topic; I would expect that as a baseline you may wish to compute \prod\_{i=1}^{n} (f\_i(x)-c\_i) as a single field over x \in \Omega = \cup\_{i=1}^n \Omega\_i, and then [call](https://docs.makie.org/stable/reference/plots/volume) their `volume(product_field, algorithm = :iso, isorange = 2*eps(eltype(product_field)), isovalue = 0.0)` method or the contour method on that page.
---
### Author:  [@WuSiren](https://discourse.julialang.org/u/WuSiren)
#### Post date: [June 5, 2024, 10:58am UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/11 "2024-06-05T10:58:49Z")
Thank you! I think your suggestion is likely to be very helpful! I’ll have a try! 🤝 🤝
---
### Author:  [@WuSiren](https://discourse.julialang.org/u/WuSiren)
#### Post date: [June 5, 2024, 11:07am UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/12 "2024-06-05T11:07:33Z")
> [@cmarcotte](#):
>
> I would expect that as a baseline you may wish to compute \prod\_{i=1}^{n} (f\_i(x)-c\_i) as a single field over x \in \Omega = \cup\_{i=1}^n \Omega\_i
But why is it \prod\_{i=1}^{n} (f\_i(x)-c\_i)? ~~Maybe what you mean is \sum\_i|f\_i(x) - c\_i|=0?~~ Don’t quite understand. 🥴
---
### Author:  [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)
#### Post date: [June 6, 2024, 9:26am UTC](https://discourse.julialang.org/t/intersection-of-level-sets/114913/13 "2024-06-06T09:26:57Z")
Yes, perhaps that is better.