Intersection of level sets

General Usage Performance
,
1

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.

1 Like
2

I’d build spatial trees for each contour so you can quickly find which segments might intersect and avoid most intersection tests.

3

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?

4

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

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.

1 Like
6

Thank you very much. I’m not familiar with the direct array approach – is this demonstrated in Meshes.jl or elsewhere?

7

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.

2 Likes
8

As a note for those who stumble on this thread in the future:
I recommend using the 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):

test
test1920Ă—1440 535 KB

1 Like
9

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!

10

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 their volume(product_field, algorithm = :iso, isorange = 2*eps(eltype(product_field)), isovalue = 0.0) method or the contour method on that page.

1 Like
11

Thank you! I think your suggestion is likely to be very helpful! I’ll have a try! :handshake: :handshake:

12

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. :woozy_face:

13

Yes, perhaps that is better.

1 Like