Creating random configurations of spherocylinders or cylinders

About the setting:

We have a 3D simulation box with side $l$ and our catesian coordinate system is set with its origin at the centre of the box. We have a number $N$ of spherocylinders of aspect ratio $L/D=20,$ where $D$ is our basic unit of length, namely, the base diameter of the cylinder (or the spherical caps) and $L$ the length of the cylinder (without the caps). We can take an example of a box with $l=100.$

Question:

The aim is to generate random configurations for the $N$ spherocylinders, such that none of them are in overlap, so not even touching, nor are they in contact with the walls of the box. This is useful for studying packing properties of such anisotropic objects, in contrast to the well-known sphere packing studies. I am curious to find out whether in Mathematica, there are existing libraries that can help one get a head start for implementing a routine that would generate valid configurations. A configuration is simply a list with the coordinates of the centres of the particles and their respective orientation coordinates, so 6 numbers per particle.

• The shape you want is called CapsuleShape[] in Mathematica. – J. M. is in limbo Nov 8 '17 at 17:31
• @J.M. Exactly, thank you. – user929304 Nov 8 '17 at 17:36
• – Jason B. Nov 8 '17 at 18:29

3 Answers

Okay, I could not resist to give it another try. This is somewhat related to my research anyway. I use a new answer since this is going to be quite a different approach and there might still something to learn from the old approach.

This time, we write a heavily optimized CompiledFunction, called CapsuleIntersectingQ that (hopefully correctly) detects wether two capsules with given centerlines and given radii intersect. We do that by discussing a convex quadratic function f (see below) in two variables that describes quadratic distances of points on the infinite lines through the centerlines of the capsules. Actually, we need the minimizers of this function on the unit square, so we need to discuss several cases (9 cases to be precise, depending on wether the minimum lies in the interior, on an edge, or in a corner of the unit square). The generic case can solved by applying Newton's method symbolically (only one step is needed, since f is quadratic). I am not 100% sure wether I handle the boundary cases correctly, so please report if you find any intersecting capsules where they should not be or any other oddities.

(* a CompiledFunction that is supposed to test pairs of capsules for intersection *)
CapsuleIntersectingQ =
Quiet[Block[{pp, qq, f, Df, DDf, x, y, S, T, S0, S1, T0, T1, s, t, p1, p2, p3, p4, p5, p6, q1, q2, q3, q4, q5, q6},
pp = {{p1, p2, p3}, {p4, p5, p6}};
qq = {{q1, q2, q3}, {q4, q5, q6}};
f = {x, y} \[Function] Evaluate[Total[({(1 - x), x}.pp - {(1 - y), y}.qq)^2]];
Df = {x, y} \[Function] Evaluate[D[f[x, y], {{x, y}, 1}]];
DDf = {x, y} \[Function] Evaluate[D[f[x, y], {{x, y}, 2}]];
{S, T} = LinearSolve[DDf[0, 0], -Df[0, 0]];
S0 = -Df[0, 0][[1]]/DDf[0, 0][[1, 1]];
S1 = -Df[0, 1][[1]]/DDf[0, 1][[1, 1]];
T0 = -Df[0, 0][[2]]/DDf[0, 0][[2, 2]];
T1 = -Df[1, 0][[2]]/DDf[1, 0][[2, 2]];
With[{
sint = N@S, tint = N@T,
F00 = N@f[0, 0], F01 = N@f[0, 1], F10 = N@f[1, 0],
F11 = N@f[1, 1], Fst = N@f[S, T],
s0 = N@S0, Fs0 = N@f[S0, 0],
s1 = N@S1, Fs1 = N@f[S1, 1],
t0 = N@T0, F0t = N@f[0, T0],
t1 = N@T1, F1t = N@f[1, T1]
},
Compile[{{p, _Real, 2}, {r1, _Real}, {q, _Real, 2}, {r2, _Real}},
Block[{s, t, c, p1, p2, p3, p4, p5, p6, q1, q2, q3, q4, q5, q6},
p1 = CompileGetElement[p, 1, 1];
p2 = CompileGetElement[p, 1, 2];
p3 = CompileGetElement[p, 1, 3];
p4 = CompileGetElement[p, 2, 1];
p5 = CompileGetElement[p, 2, 2];
p6 = CompileGetElement[p, 2, 3];
q1 = CompileGetElement[q, 1, 1];
q2 = CompileGetElement[q, 1, 2];
q3 = CompileGetElement[q, 1, 3];
q4 = CompileGetElement[q, 2, 1];
q5 = CompileGetElement[q, 2, 2];
q6 = CompileGetElement[q, 2, 3];
c = Min[{F00, F01, F10, F11}];
If[c <= (r1 + r2)^2,
True,
s = sint; t = tint;
If[0. <= s <= 1. && 0. <= t <= 1.,
c = Fst,
If[0. <= s0 <= 1., c = Min[c, Fs0]];
If[0. <= s1 <= 1., c = Min[c, Fs1]];
If[0. <= t0 <= 1., c = Min[c, F0t]];
If[0. <= t1 <= 1., c = Min[c, F1t]];
];
c <= (r1 + r2)^2
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]
]
]];


For completeness, here is our method to produce random rotation matrices from random vectors.

getRotationMatrix =
Compile[{{u, _Real, 1}},
Block[{uu, r, cosr, sinr}, uu = u[[1]]^2 + u[[2]]^2 + u[[3]]^2;
r = Sqrt[uu];
cosr = Cos[r];
sinr = Sin[r];
{{(u[[1]]^2 + cosr u[[2]]^2 + cosr u[[3]]^2)/
uu, (u[[1]] u[[2]] - cosr u[[1]] u[[2]] - u[[3]] r sinr)/
uu, (u[[1]] u[[3]] - cosr u[[1]] u[[3]] + u[[2]] r sinr)/
uu}, {(u[[1]] u[[2]] - cosr u[[1]] u[[2]] + u[[3]] r sinr)/
uu, (cosr u[[1]]^2 + u[[2]]^2 + cosr u[[3]]^2)/
uu, (u[[2]] u[[3]] - cosr u[[2]] u[[3]] - u[[1]] r sinr)/
uu}, {(u[[1]] u[[3]] - cosr u[[1]] u[[3]] - u[[2]] r sinr)/
uu, (u[[2]] u[[3]] - cosr u[[2]] u[[3]] + u[[1]] r sinr)/
uu, (cosr u[[1]]^2 + cosr u[[2]]^2 + u[[3]]^2)/uu}}],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];


And here is the main program. Note that I remove everything related to RegionDistance.

(*number of capsules to stuff into the box*)
n = 10000;
(*maximal number of trials to stuff a random capsule into the box*)
maxiters = 1000;
(*edge length of the box*)
l = 100.;
(*length of capsules*)
L = 15;
(*radius of capsules*)
r = 3.;
(*prototypical centerline of capsules; new capsules are generated by rotating and translating it randomly*)
p0 = DeveloperToPackedArray[{-L/2 {1., 0., 0.}, L/2 {1., 0., 0.}}];
(*a bag for rule them all*)
ptbag = InternalBag[{}];
(*initialize the distance function so that the first capsule will be feasible*)
distToCenterlines = Function[x, 2 l, Listable];
(*a counter for the capsules*)
ncapsules = 0;
(*the box (for graphical purposes only*)
cube = Cuboid[-l/2 {1., 1., 1.}, l/2 {1., 1., 1.}];
(*the box where the centerlines have to be contained in*)
safetycube = Cuboid[-(l/2 - r) {1., 1., 1.}, (l/2 - r) {1., 1., 1.}];
(*generate a function that test if a point is contained in safetycube*)
insafetycubeQ = RegionMember[safetycube];
(*this is the amount of random data to generate at once*)
m = 1000;
(*c is a counter keeping track of the amount of random data (matlist and shiftlist below) that has been used already;initialize the counter so that a refresh will be kicked off immediately in the main loop*)
c = m;
(*use a normal distribution to generate random rotation vectors*)
distro = MultinormalDistribution[{0., 0., 0.}, DiagonalMatrix[{1., 1., 1.}]];
Dynamic[{j, iter}]
(*finally,the main loop*)
Do[iter = 0;
(*Boolean that is supposed to be equal to False if the new capsule is feasible*)
b = True;
(*generate random capsules until you find a feasible one*)
While[b,
++iter;
If[iter >= maxiters, Break[]];
(*check if the randomness reservoir is not empty*)
If[c >= m,(*if empty,refill the reservoir*)
matlist = getRotationMatrix[RandomVariate[distro, m]];
shiftlist = RandomReal[{r - l/2, l/2 - r}, {m, 3}];
c = 0;
];
c++;
(*generate the random capsule*)
p = p0.matlist[[c]] + ConstantArray[shiftlist[[c]], 2];
(*perform the cheap test against the box first*)
b = ! And @@ (insafetycubeQ[p]);
If[! b && ncapsules > 0,
centerlines = ArrayReshape[InternalBagPart[ptbag, All], {ncapsules, 2, 3}];
b = Or @@ CapsuleIntersectingQ[centerlines, r, p, r]
];
];
If[iter >= maxiters,
Print["I tried my best but with ", maxiters, " trials, I was able to stuff only ", ncapsules, " capsules into the box."];
Break[]
];
(*yes,we have one more capsule...*)
ncapsules += 1;
(*stuff the centerline of new capsule into the bag*)
InternalStuffBag[ptbag, Flatten[p], 6];
, {j, 1, n}]


This time, I let it run until nothing seems to fit in any more. Here is the result:

Graphics3D[{
FaceForm[{Darker@Blue, Opacity[.1]}], Specularity[White, 30],
cube, safetycube,
FaceForm[{Darker@Orange, Opacity[1.0]}], Specularity[White, 30],
Map[CapsuleShape[#, r] &, centerlines]
},
Lighting -> "Neutral", PlotRange -> Transpose[List @@ cube],
Boxed -> False, SphericalRegion -> True
]


• Feel free to experiment with it. You merely have to change the part where p is generated from p0. You can also use other random distributions of orthogonal matrices. Look for MatrixPropertyDistribution in the documentation and see how JasonB applied it. You can of course also alter the distribution for the shift vectors. I am satisfied for know. (If there are no bugs in the code.) Well, maybe I will think about applying a tree-structure to get the complexity down from $O(n^2)$ to $O(n \log(n))$... – Henrik Schumacher Nov 9 '17 at 14:34
• Alright, indeed it s a very fun problem. Thanks again for all your patience and dedication, much much appreciated. I guess a naive way of fixing the orientations would be to fix the value of c in p=p0.matlist[[c]]+... – user929304 Nov 9 '17 at 14:43
• Yeah, or simply replacing matlist[[c]] by a rotation matrix of your choice... – Henrik Schumacher Nov 9 '17 at 14:44
• Henrik, thank you again for your patience and all the kind and prompt replies, really really appreciated. I have learnt a lot from our discussions and your suggestions. – user929304 Nov 10 '17 at 13:17
• You're welcome! – Henrik Schumacher Nov 10 '17 at 13:34

Any approach to this will depend on efficiently deciding if two capsules overlap. You could do this using the built-in RegionDisjoint

capsulesOverlapRegion[cap1 : CapsuleShape[{p1a_List, p1b_List}, d1_],
cap2 : CapsuleShape[{p2a_List, p2b_List}, d2_]] := !
RegionDisjoint[
BoundaryDiscretizeRegion[cap1],
BoundaryDiscretizeRegion[cap2]
]


but this isn't terribly fast.

SeedRandom[42];
{c1, c2} = {CapsuleShape[RandomReal[1, {2, 3}], .25],
CapsuleShape[RandomReal[1, {2, 3}], .25]};


capsulesOverlapRegion[c1, c2] // RepeatedTiming
(* {0.033, True} *)


Instead, I will use the compiled function described in this post, which measures the minimum distance between two line segments. If this distance is greater than the sum of the diameters, the capsules do not intersect. Here I'll drop the CapsuleShape wrapper, and represent a given capsule as a list of {point1, point2, diameter}

distSegToSegComp = << "https://git.io/vDJMK";
capsulesOverlap[CapsuleShape[{p1a_List, p1b_List}, d1_],
CapsuleShape[{p2a_List, p2b_List}, d2_]] := (d1 + d2) >
distSegToSegComp[{p1a, p1b}, {p2a, p2b}]

capsulesOverlap[c1, c2] // RepeatedTiming
(* {9.*10^-6, True} *)


Now you can use a loop, generating a random capsule at each step and rejecting it if it overlaps with another capsule.

With a speedy test, we can just use a While loop to get the required number of shapes. Here is an example that generates 300 pseudorandomly oriented capsules, with length 0.3 and diameter 0.05, all within a box extending from {-1,-1,-1} to {1,1,1}.

length = 0.3;
diameter = 0.05;
nCapsuleGoal = 300;
maxAttempts = 3000;

boxLowerPoint = {-1, -1, -1};
boxUpperPoint = {1, 1, 1};
box = BoundaryDiscretizeRegion[
Cuboid[diameter + boxLowerPoint, boxUpperPoint - diameter]];
pointInBox =
And @@ Thread[(diameter + boxLowerPoint) < # < (-diameter +
boxUpperPoint)] &;

distSegToSegComp = << "https://git.io/vDJMK";
capsulesOverlap[{p1a_List, p1b_List, d1_}, {p2a_List, p2b_List,
d2_}] := (d1 + d2) > distSegToSegComp[{p1a, p1b}, {p2a, p2b}]

capsuleInBox[{p1_List, p2_List, d_}] :=
pointInBox[p1] && pointInBox[p2];

\[ScriptCapitalD] =
MatrixPropertyDistribution[r.{0, 0, length},
r \[Distributed] CircularRealMatrixDistribution[3]];

capsules = {};
nCapsules = 0;
nAttempts = 0;
randomCapsule[] :=
With[{start = RandomPoint[box]}, {start,
start + RandomVariate[\[ScriptCapitalD]], diameter}];

While[nCapsules < nCapsuleGoal && nAttempts < maxAttempts,
capsule = randomCapsule[];
If[capsuleInBox[capsule] &&
AllTrue[capsules, ! capsulesOverlap[capsule, #] &],
AppendTo[capsules, capsule];
nCapsules++;];
nAttempts++;];


visualize the results with

Graphics3D[{CapForm[Round], Tube[{#1, #2}, #3] & @@@ capsules, Blue,
Opacity[0.1], Cuboid[boxLowerPoint, boxUpperPoint]}, Boxed -> False]


• Ah, iteresting. I did not know about MatrixPropertyDistribution until now. – Henrik Schumacher Nov 8 '17 at 19:47
• Unfortunately, I cannot test your code for I am running on version 11.0.1 and RegionWithin was new in 11.1 . It's time to get the new version, though... – Henrik Schumacher Nov 8 '17 at 19:50
• @user929304 - I removed reference to RegionWithin, and changed the return to give data instead of CapsuleShape objects. After running the code, the symbol capsules is a list of the form {{x1, y1, z1}, {x2, y2, z2}, diameter}. It should be straightforward for you to change that to suit your needs. – Jason B. Nov 8 '17 at 22:47
• I don't follow your question about normalization though. I just used randomly chosen parameters for the length of the capsules, and the size of the box. You can set these however you like. – Jason B. Nov 8 '17 at 22:47
• @user929304 - that makes sense, physically you can reach a case where there aren't enough cavities big enough to find randomly. I don't know what methods people use to get around that. – Jason B. Nov 9 '17 at 15:05

Disclaimer

In the meantime, user929304 found out that my collision test is incorrect. I posted a new answer with a rather new approach.

Strategy

I use a Bag to store the lines between the two sphere centers of each existing capsule. When generating a new capsule I use

1. RegionMember to check wether a new capsule fits into the box and
2. RegionDistance of a MeshRegion containing all centerlines of existing cpasules to check for collisions between the new capsule and the existent capsules.

It is advantageous to pack all existing capsules into a single MeshRegion and to use its RegionDistanceFunction since internally, Mathematica uses a k-d-tree (or a similar data structure): This avoids computing the distance from a new capsule to each preexisting capsule.

Here is the code along with several comments:

(*number of capsules to stuff into the box*)
n = 100;
(*maximal number of trials to stuff a random capsule into the box*)
maxiters = 1000;
(*edge length of the box*)
l = 10.;
(*length of capsules*)
L = 0.5;
(*radius of capsules*)
r = 0.125;
(*prototypical centerline of capsules; new capsules are generated by rotating and translating it randomly*)
p0 = DeveloperToPackedArray[{-L/2 {1., 0., 0.}, L/2 {1., 0., 0.}}];
(*a bag for rule them all*)
ptbag = InternalBag[{}];
(*generating edge combinatorics in advance safes much time later*)
edges = Transpose[{2 Range[n] - 1, 2 Range[n]}];
(*initialize the distance function so that the first capsule will be feasible*)
distToCenterlines = Function[x, 2 l, Listable];
(*a counter for the capsules*)
ncapsules = 0;
(*the box (for graphical purposes only*)
cube = Cuboid[-l/2 {1., 1., 1.}, l/2 {1., 1., 1.}];
(*the box where the centerlines have to be contained in*)
safetycube = Cuboid[-(l/2 - r) {1., 1., 1.}, (l/2 - r) {1., 1., 1.}];
(*generate a function that test if a point is contained in safetycube*)
insafetycubeQ = RegionMember[safetycube];
(*the following is a 1000 times faster substitute for RotationMatrix; sorry for the mess*)
getRotationMatrix = Compile[{{u, _Real, 1}},
Block[{uu, r, cosr, sinr},
uu = u[[1]]^2 + u[[2]]^2 + u[[3]]^2;
r = Sqrt[uu];
cosr = Cos[r];
sinr = Sin[r];
{{
(u[[1]]^2 + cosr u[[2]]^2 + cosr u[[3]]^2)/uu,
(u[[1]] u[[2]] - cosr u[[1]] u[[2]] - u[[3]] r sinr)/
uu, (u[[1]] u[[3]] - cosr u[[1]] u[[3]] + u[[2]] r sinr)/uu
}, {
(u[[1]] u[[2]] - cosr u[[1]] u[[2]] + u[[3]] r sinr)/uu,
(cosr u[[1]]^2 + u[[2]]^2 + cosr u[[3]]^2)/uu,
(u[[2]] u[[3]] - cosr u[[2]] u[[3]] - u[[1]] r sinr)/uu
}, {
(u[[1]] u[[3]] - cosr u[[1]] u[[3]] - u[[2]] r sinr)/
uu, (u[[2]] u[[3]] - cosr u[[2]] u[[3]] + u[[1]] r sinr)/uu,
(cosr u[[1]]^2 + cosr u[[2]]^2 + u[[3]]^2)/uu
}}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];

(*this is the amount of random data to generate at once*)
m = 1000;
(* c is a counter keeping track of the amount of random data (matlist and shiftlist below) that has been used already; initialize the counter so that a refresh will be kicked off immediately in the main loop*)
c = m;
(*use a normal distribution to generate random rotation vectors*)
distro = MultinormalDistribution[{0., 0., 0.},DiagonalMatrix[{1., 1., 1.}]];
Dynamic[{j, iter}]
(*finally, the main loop*)
Do[iter = 0;
(*a Boolean that is supposed to be equal to False if the new capsule is feasible*)
b = True;
(*generate random capsules until you find a feasible one*)
While[b,
++iter;
If[iter >= maxiters, Break[]];
(*check if the randomness reservoir is not empty*)
If[c >= m,
(*if empty, refill the reservoir*)
matlist = getRotationMatrix[RandomVariate[distro, m]];
shiftlist = RandomReal[{r - l/2, l/2 - r}, {m, 3}];
c = 0;
];
c++;
(*generate the random capsule*)
p = p0.matlist[[c]] + ConstantArray[shiftlist[[c]], 2];
(*perform the cheap test against the box first*)
b = ! And @@ (insafetycubeQ[p]);
If[! b, b = Min[distToCenterlines[p]] <= 2. r];
];
(*when the while loop exits, p represents a feasible capsule or we used maxiters trials*)
If[iter >= maxiters,
Print["I tried my best but I could only stuff ", ncapsules,
" into the box."];
Break[]
];
(*yes,we have one more capsule...*)
ncapsules += 1;
(*stuff the centerline end of new capsule into the bag*)
InternalStuffBag[ptbag, Flatten[p], 6];
(*recompute distance function to existing centerlines*)
distToCenterlines = RegionDistance[
MeshRegion[
Partition[InternalBagPart[ptbag, All], 3],
Line[edges[[1 ;; ncapsules]]]
]
];
, {j, 1, n}]


Extract the information as desired from by the OP:

capsulecenterlines = ArrayReshape[InternalBagPart[ptbag, All], {ncapsules, 2, 3}];
capsulecenters =  0.5 (capsulecenterlines[[All, 2]] + capsulecenterlines[[All, 1]]);
capsuledirections = (capsulecenterlines[[All, 2]] - capsulecenterlines[[All, 1]])/L;


These can be exported, e.g., to CSV files...

Export[FileNameJoin[{$HomeDirectory, "capsulecenterdata.csv"}], capsulecenters] Export[FileNameJoin[{$HomeDirectory, "capsuledirectiondata.csv"}], capsuledirections]


... and also reimported (with check for correctness):

Import[FileNameJoin[{$HomeDirectory, "capsulecenterdata.csv"}], "Data"] == capsulecenters Import[FileNameJoin[{$HomeDirectory, "capsuledirectiondata.csv"}], "Data"] == capsuledirections

(* True *)
(* True *)


And here is the result in graphical form:

Graphics3D[{
FaceForm[{Darker@Blue, Opacity[.1]}], Specularity[White, 30],
cube, safetycube,
FaceForm[{Darker@Orange, Opacity[1.0]}], Specularity[White, 30],
Map[CapsuleShape[#, r] &, capsulecenterlines]
},
Lighting -> "Neutral",
PlotRange -> Transpose[List @@ cube],
Boxed -> False,
SphericalRegion -> True
]


Suggestions for improvement:

I think it would be possible to speed things up if one could get a grip directly on the k-d-Tree. In the current implementation, the tree is rebuilt completely after a new capsule is added. But there has to be cheaper a method to create the new k-d-Tree by modifying the old one.

I am not 100% sure if the generated rotation matrices a uniformly distributed over rotation group $\operatorname{SO}(3)$. On the other hand, the generated distribution is rotation invariant, isn't it? So we sample from the normalized Haar measure, aren't we?

Edit

In the meantime, I made several improvements:

• The greatest effect came from defining the edge combinatorics right in the beginning so that a lot of messy copying, packing and unpacking of arrays can be skipped.

• Generate random data in larger chunks since this is a bit faster. That has also some effect.

• Tried to use an InternalBag for collecting the endpoints of the center lines. The effect however is almost not measurable.

• Replaced RotationMatrix with a much quicker compiled function.

• Very nice, and fast as well – Jason B. Nov 8 '17 at 19:42
• @JasonB. The idea was to use the KDTree in RegionsDistance such that a new capsule need not to be tested against all existing capsules... – Henrik Schumacher Nov 8 '17 at 19:44
• Time is always short ;o) I wrote some comments into the code and also extracted the center coordinates and orientation vectors for you. – Henrik Schumacher Nov 8 '17 at 22:44
• You may increase m but it is not a must. Large m may improve the performance marginally. In the meantime I also fixed the issue with the first "dummy" capsule. I needed it because Mathematica does not like empty MeshRegions. – Henrik Schumacher Nov 9 '17 at 10:05
• Huh, you found overlaps? That's bad. Precision issues are almost surely not the reason. Must be a coding error. I have just edited the presented code several times. Could you please test it again and tell me if this problem persists? – Henrik Schumacher Nov 9 '17 at 10:13