4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Do all the cylinders have the same radius? $\endgroup$ – J. M. is away Nov 27 '15 at 16:03
  • $\begingroup$ @J.M. No, unfortunately. The radius is controlled by RandomReal[5] in Table[{RandomReal[{-100, 100}, {2, 3}], RandomReal[5]}, {100}]. $\endgroup$ – Dimitris Nov 27 '15 at 16:09
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M. is away Nov 27 '15 at 16:41
  • $\begingroup$ Ok. I will try to do on my own and post an answer. Thanks for the advice. $\endgroup$ – Dimitris Nov 27 '15 at 16:46
  • $\begingroup$ @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? $\endgroup$ – Dimitris Dec 1 '15 at 9:38
8
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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.) $\endgroup$ – Rahul Nov 28 '15 at 3:12
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Nov 28 '15 at 3:23
1
$\begingroup$

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}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.