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 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. This is seemingly robust, but stunningly slow.

A complete gist is available here. 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:

```
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.