# 3D intersection algorithm for cylinders

The problem

The input is a list of $$N$$ cylinders in 3D space, and the output should be a list of $$M \leq N(N-1)/2$$ pairs of cylinders that intersect. ($$M$$ depends on the input data, obviously.)

If it matters, the cylinders are very thin (with a diameter less than 1% of the length for all cylinders), and a solution for "rounded cylinders" would work for me (it probably simplifies the geometry calculations). A "rounded cylinder" is a cylinder with half-spheres at the extremities; formally, for a starting point $$S$$, an endpoint $$E$$ and a radius $$r$$, the rounded cylinder $$(S, E, r)$$ is defined as the set $$\{P \mid \exists Q \in [S,E], \|PQ\| \leq r\}$$.

The obvious solution

It is easy enough to do in $$O(N^2)$$ time and $$O(max(M, N))$$ space: the pseudocode of my current implementation is (for rounded cylinders):

Ncyl = length(cylinder_list)
output = {}
for i = 1, 2, ... Ncyl:
for j = i+1, i+2, ... Ncyl:
(S1, E1, r1) = cylinder_list[i]
(S2, E2, r2) = cylinder_list[j]
find P∈[S1, E1], Q∈[S2, E2] such that ||PQ|| is minimal  # this is the costly line, says the profiler
if ||PQ|| < r1 + r2:
return output


Better performance?

Any algorithm will have a (time and space) worst-case in $$O(N^2)$$ (at least) because the output list can itself be of that length. However, the above algorithm guarantees $$O(N^2)$$ time even on "friendly" data because it tests all possible intersections.

In my use case, the cylinders are fairly spread apart in space (the longest cylinder is less than one tenth of the diameter of the whole set of cylinders). Furthermore, they occupy a small fraction of space and $$M\sim N$$ (for values of $$N$$ up to 2000 or so, above that it times out). This suggests to me that there could be an improvement by some "sweeping plane" algorithm similar to Bentley-Ottmann. However, I did not find a straightforward way to do Bentley-Ottmann in 3D (in 2D after sweeping you end up ordering points on a line which is easy enough, but in 3D there is no obvious ordering for a plane).

• O(n²) complexity can be be avoided with data structures that eliminate interaction between objects that are spread apart. k-d trees for example. – d3m4nz3 May 12 '20 at 14:48
• Minimal distance between two segments in 3D can be found quickly - please see: math.stackexchange.com/questions/1873911/… – HEKTO May 17 '20 at 20:42
• @d3m4nz3 k-d trees gave me the pointer I needed. Should I write an answer myself or do you want to go ahead? – JMU May 18 '20 at 13:38
• HEKTO: yes, I know. My implementation of segment distance is basically the same as in your link; maybe it is not optimal but it certainly is within a factor 2 of the optimal; so I am pretty sure I need to attack the complexity (via the number of times the segment distance is computed), not the constant. – JMU May 18 '20 at 13:43
• If you know how to calculate the shortest distance between segments, then I think the trivial algorithm will work fine. Anyway, you have only 2000 segments - these hierarchical data structures and sweeping make sense for millions of segments – HEKTO May 18 '20 at 15:15

If you could accept an approximate solution for this problem, then the following approach would help you. Your cylinders are thin, so if two of them intersect, then their axial segments almost touch each other. If you project such a pair of almost touching segments on each of coordinate planes you'll see that their projections almost always intersect.

So - project your set of segments on three coordinate planes and apply the Bentley-Ottmann (or similar) algorithm three times to find all the intersections. Select pairs of segments, which intersect on all three planes, calculate distance between them and verify the cylinders intersection.

This method can give you false negatives - for example in the case, when two cylinders are collinear with small gap between them. However, this method won't give you false positives, because you explicitly verify each suspicious pair of cylinders.

• It is tricky to see in which configurations that algorithm is an improvement. The problem is that you might have a lot of false negatives so that the plane intersection steps are O(N^2) (or worse), even if there are very few intersections in volume. For instance, in a box 1x1x1, with cylinders of random orientation, length 1 and very low diameter (say 10^-100), there are very few intersections in volume, but (probably?) O(N^2) intersection in projection planes. On the other hand, if cylinders have a strong preferential orientation, one projection plane might prune the search space enough. – JMU May 24 '20 at 13:17
• (edit) The idea can be improved slightly to eliminate false negatives: after each projection and before calculating in-plane intersections, lengthen each segment by the cylinder diameter on either end (this might create more false positives, but very few when the aspect ratio of the cylinders is >100). – JMU May 24 '20 at 13:22
• A plausible variant of this idea is to pick a random projection to 2D (i.e., pick a random pair of vectors, then project onto the plane defined by those two vectors, and solve the 2D problem). – D.W. May 24 '20 at 17:39

@d3m4nz3 pointed in a comment to the use of k-d trees. That is the solution for the "intersection of spheres" problem, and I suspect it is (asymptotically) optimal when the cylinders are thick. The pseudocode would be something like:

max_sphere_rad = max( ||SE|| + r for (S, E, r) in cylinder_list)
tree = (build kd tree of cylinders based on their center location)

output = {}
loop for (Si, Ei, ri) in cylinder_list
potential_intersections = query tree for cylinders whose center is at most at distance (||Si Ei||/2 + r + max_sphere_rad) of (Si+Ei)/2
loop for (Sj, Ej, rj) in potential_intersections
if (Si, Ei, ri) and (Sj, Ej, rj) intersect:

Building the kd tree is $$O(N log N)$$. The question is whether it really speeds up the subsequent range search for potential_intersections. That in turn depends on the exact way the kd tree was built vs. the distribution of cylinder center points etc., but the rough idea is that if the output is of size $$K$$, the query has complexity $$O(K + log N)$$ (for a reasonably balanced kd tree). In particular:
• if $$K\leq O(log N)$$ (very few potential intersection points per cylinder), the whole algorithm runs in $$O(N log N)$$
• If $$K=O(N)$$ (large portions of the dataset interact), the query will always be at least $$O(N)$$, making the whole algorithm run in $$O(N^2)$$
For a rough analysis of the impact of the cylinder distribution on the algorithm complexity, assume $$N$$ cylinders with constant length $$L$$ (and radii well below $$L$$), that are uniformly-distributed in a volume $$V$$. A typical cylinder will have $$O(N L^3/V)$$ neighbors in potential_intersections. Based on the above, such neighbors can be extracted in close to $$O(N L^3/V)$$ time and the whole thing runs then in $$O(N^2 L^3/V)$$. At constant $$L^3/V$$ (i.e. add more cylinders in a constant-volume space), the algorithm is quadratic in the number of cylinders; at constant $$N/V$$ (i.e. increase volume size and add cylinders at iso-density), it is linear.