# generating randomly oriented non-intersecting cylinders

I could not find anything relevant with "googling" (only this)

I can create a set of randomly oriented cylinders (code found here)

cylinders =
Table[{RandomReal[{-100, 100}, {2, 3}], RandomReal[5]}, {100}];
Graphics3D[{EdgeForm[None],
Directive[Opacity@RandomReal[{.4, .9}], Hue[RandomReal[]]],
Cylinder[First@#, Last@#]} & /@ cylinders, Boxed -> False,
ImageSize -> 800]


How can we keep the cylinders that they do not intersect each other?

• Do all the cylinders have the same radius? Nov 27, 2015 at 16:03
• @J.M. No, unfortunately. The radius is controlled by RandomReal[5] in Table[{RandomReal[{-100, 100}, {2, 3}], RandomReal[5]}, {100}]. Nov 27, 2015 at 16:09
• You can adapt the rejection method used here to this case; all you need is a way to find the distance between two line segments, such as the routine in this answer. Nov 27, 2015 at 16:41
• Ok. I will try to do on my own and post an answer. Thanks for the advice. Nov 27, 2015 at 16:46
• @J.M. If I understand correctly, I will use the routine in order to generate no intersecting line segments which will be axes of cylinders? Dec 1, 2015 at 9:38

This answer replaces an earlier one that deleted only some of the intersecting cylinders. (My thanks to paw for pointing this out.) It also is much faster.

The square of the distance between points at p1 and p2 is

p1.p1 + p2.p2 - 2 p1.p2


and a cylinder axis can be parameterized by pi + dp t, where pi is one end of the axis, dp is the vector from pi to the other end of the axis, and 0 < t < 1. Thus, the parameterized square of the distance between points on the axes of two cylinders is

((p1.p1 + p2.p2 - 2 p1.p2) /. {p1 -> p1i + dp1 t1, p2 -> p2i + dp2 t2}) // tf // Expand
(* p1i.p1i - 2 p1i.p2i + p1i.(dp1 t1) - 2 p1i.(dp2 t2) + p2i.p2i +
p2i.(dp2 t2) + (dp1 t1).p1i - 2 (dp1 t1).p2i + (dp1 t1).(dp1 t1) -
2 (dp1 t1).(dp2 t2) + (dp2 t2).p2i + (dp2 t2).(dp2 t2) *)


where

tf[e_] := e /. Dot[z1_ + z2_, z3_ + z4_] ->
Dot[z1, z3] + Dot[z2, z3] + Dot[z1, z4] + Dot[z2, z4]


is used to distribute Dot over Plus. (FunctionExpand does not do so.) Next, the scalars t1 and t2 need to be factored from Dot.

Map[# /. z_ :> (z /. t1 -> 1) t1^Count[z, t1, 4] &,
Map[# /. z_ :> (z /. t2 -> 1) t2^Count[z, t2, 4] &, %]]
(* t1^2 dp1.dp1 - 2 t1 t2 dp1.dp2 + t1 dp1.p1i - 2 t1 dp1.p2i + t2^2 dp2.dp2 +
t2 dp2.p2i + t1 p1i.dp1 - 2 t2 p1i.dp2 + p1i.p1i - 2 p1i.p2i + t2 p2i.dp2 + p2i.p2i *)


Finally, t1 and t2 that minimize this quantity are obtained

FullSimplify[Solve[{D[%, t1] == 0, D[%, t2] == 0}, {t1, t2}],
TransformationFunctions -> {Automatic, tf1}]
(* {{t1 -> (dp2.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp2 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2),
t2 -> (dp1.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp1 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2)}} *)


where

tf1[e_] := e /. Dot[z1_, z2_] -> Dot[z2, z1]


A function determining the distance between two cylinders can be written, based on this result.

int[cyl1_, cyl2_] :=
Module[{p1i = cyl1[[1, 1]], dp1 = cyl1[[1, 2]] - cyl1[[1, 1]], r1 = cyl1[[2]],
p2i = cyl2[[1, 1]], dp2 = cyl2[[1, 2]] - cyl2[[1, 1]], r2 = cyl2[[2]],
loc, t1, t2}, loc = {t1 -> (dp2.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp2 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2),
t2 -> (dp1.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp1 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2)};
(-r2/Norm[dp2] < (t1 /. loc) < 1 + r2/Norm[dp2]) &&
(-r1/Norm[dp1] < (t2 /. loc) < 1 + r1/Norm[dp1]) &&
((t1^2 dp1.dp1 - 2 t1 t2 dp1.dp2 + t1 dp1.p1i - 2 t1 dp1.p2i + t2^2 dp2.dp2 +
t2 dp2.p2i + t1 p1i.dp1 - 2 t2 p1i.dp2 + p1i.p1i - 2 p1i.p2i + t2 p2i.dp2 +
p2i.p2i) /. loc) < (r1 + r2)^2]


The conditions (-r2/Norm[dp2] < (t1 /. loc) < 1 + r2/Norm[dp2]) and similarly for t2 handle the ends of the cylinders. With this distance function completed, the rest of the computation is straightforward.

cylinders = Table[{RandomReal[{-100, 100}, {2, 3}], RandomReal[5]}, {100}];
nint = ParallelTable[Or @@ Table[int[cylinders[[i]], cylinders[[j]]], {j, i + 1, 100}],
{i, 100}];
c = Cases[{cylinders, nint} // Transpose, {z_, False} -> z, {1}];
Graphics3D[{EdgeForm[None], Directive[Opacity@RandomReal[{.4, .9}], Hue[RandomReal[]]],
Cylinder[First@#, Last@#]} & /@ c, Boxed -> False, ImageSize -> 800]


In this particular run, only 45 of the cylinders remain. The entire calculation takes only seconds.

• I'm afraid your code thinks that the cylinders Cylinder[{{-1, 0, 0.1}, {0, 0, 0}}, 1] and Cylinder[{{1, 0, 0}, {2, 0, 0.1}}, 1] are intersecting. (I had to tilt them slightly otherwise it gave me a 1/0 error.)
– user484
Nov 28, 2015 at 3:12
• @Rahul You are correct on both counts. Collinear axes cause problems, and the intersection criterion at the ends errs on the side of caution. There probably is a simply way to handle the former. The latter is deliberate in order to avoid extremely complicated expressions for handling cylinder ends. I prefer to discard a few too many cylinders. Nov 28, 2015 at 3:23

Not exactly a solution and certainly not sufficient at all.

The following function can check whether two cylinders intersect or not.

   IntCyl[cyl1_, cyl2_] :=
RegionQ[DiscretizeRegion[RegionIntersection[cyl1, cyl2]] // Quiet]


For example

c1 = Cylinder[{{1, 1, 1}, {3, 5, 1}}];
c2 = Cylinder[{{-4, 5, 0}, {2, 5, 3}}];
c3 = Cylinder[{{-2, 0, 0}, {2, 0, 0}}];
c4 = Cylinder[{{0, 0, -2}, {0, 0, 2}}];
GraphicsRow@{Graphics3D[{c1, c2}], Graphics3D[{c3, c4}]}


IntCyl[c1, c2] // Timing
IntCyl[c3, c4] // Timing


{0.461932, False}
{2.11221, True}

An extension to more cylinders is not attractive at all...Even with 10 cylinders Mathematica needs a lot of time.

cylinders =
Table[{RandomReal[{-100, 100}, {2, 3}], RandomReal[5]}, {10}];
cylinderslist = Cylinder[First@#, Last@#] & /@ cylinders;
Timing[nint =
Table[Or @@
Table[IntCyl[cylinderslist[[i]], cylinderslist[[j]]], {j, i + 1,
10}], {i, 10}];]


{142.887, Null}