4

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?

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
Dimitris
  • 4,794
  • 22
  • 50

2 Answers2

8

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]

enter image description here

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

Community
  • 1
bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
  • 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.) –  Nov 28 '15 at 03: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. – bbgodfrey Nov 28 '15 at 03:23
1

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}]}

enter image description here

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}

Dimitris
  • 4,794
  • 22
  • 50