6
$\begingroup$

For a more complex collision detection approach, I need an efficient (aka fast) algorithm to prune the elements by using boundary boxes. A classical algorithm is the so called sweep and prune (SAP) algorithm (as outlined in this Web-Article). However so far I did not succeed in vectorize the SAP algorithm, so I went for an sequential implementation. Furthermore I did not suceed in taking advantage of the quoted temporal coherence by exploiting an insertionsort approach, since Mathematica sort is not open to external tweaking of the sort algorithm as far as I know. The sequential code of the SAP algorithm for one axis looks as follows.

SweepAndPruneAxis[objaxes_] := 
 Module[{axisList, activeList, collcand, j, i},
   axisList = Ordering[objaxes, All, First[#1] < First[#2] &];
   activeList = {First[axisList]};
   collcand = {};
   For[j = 2, j <= Length[axisList], j++,
     i = 1;
     While[i > 0 && i <= Length[activeList] && i <= Length[axisList],
       If[objaxes[[axisList[[j]], 1]] <= objaxes[[activeList[[i]], 2]],
         collcand = Append[collcand, {activeList[[i]], axisList[[j]]}];
         i++;
         ,
         activeList = Delete[activeList, i];
       ];
     ];
     activeList = Append[activeList, axisList[[j]]];
   ];
   Return[collcand];
 ];

Here the objaxes contains {xmin,xmax} pairs of coordinates which get analyzed with respect of overlap and a list of indices of the elements in objaxes which overlap are returned.

The following code snippet demonstrates the capabilities of the SAP algorithm.

 SeedRandom[35];
 objects = Sort /@ RandomReal[{0, 15}, {10, 2}];
 Graphics[MapIndexed[{Rectangle[{#1[[1]], -First[#2]}, {#1[[2]], -First[#2] - 0.5}], 
 Text[First[#2], {0.1, -First[#2] - 0.25}]} &, objects]]
 Sort[Sort /@ SweepAndPruneAxis[objects]]

which yields the output

SAP Algorithm output

and the sorted list of collision pairs

{{1, 3}, {1, 5}, {1, 6}, {1, 7}, {1, 8}, {1, 9}, {2, 6}, {3, 5}, {3, 7}, {3, 8}, {3, 9}, {4, 6}, {5, 7}, {5, 9}, {7, 8}, {7, 9}, {7, 10}}

This approach can be extended to e.g. 3 dimensions by simply taking the intersection of the 3 resulting SAP results for the 3 coordinate axes using

 Intersection[Sequence @@ (SweepAndPruneAxis /@ 
   (Transpose[Transpose[objboxes, {1, 3, 2}]][[1 ;; dim]])),
   SameTest -> (Sort[#1] == Sort[#2] &)
 ]

where here now objboxes is a list of 3 min/max pairs for the respective axis.

My main issue is now that an "almost" brute force implementation which is using on-board commands of Mathematica is computationally more efficient in getting the same result:

 Extract[Subsets[Range[Length[objboxes]], {2}], 
   Transpose[Position[Apply[IntervalIntersection, 
     Map[Interval, Subsets[Transpose[Transpose[objboxes, {1, 3, 2}]][[#]], {2}] & /@
       Range[dim], {3}], {2}]], {Repeated[Interval[_List], {dim}]}]]

with objboxes as described above and dim as the dimension of the objboxes sets.

This can be demonstrated by defining a "bracket" function as

FindCollisionCandidates[objboxes_, dim_: 3, sap_: True] := If[sap,
  Intersection[
   Sequence @@ (SweepAndPruneAxis /@ (Transpose[
         Transpose[objboxes, {1, 3, 2}]][[1 ;; dim]])), 
   SameTest -> (Sort[#1] == Sort[#2] &)]
  ,
  Extract[Subsets[Range[Length[objboxes]], {2}], 
   Position[
    Transpose[Apply[IntervalIntersection, 
      Map[Interval, 
       Subsets[Transpose[
            Transpose[objboxes, {1, 3, 2}]][[#]], {2}] & /@ 
        Range[dim], {3}], {2}]], {Repeated[
      Interval[_List], {dim}]}]]]

and apply this to an randomly generated dataset with some intersections.

 rval = 0.1;
 Nbox = 250;
 dim = 3;
 SeedRandom[35];
 p = Prepend[RandomReal[{-5 + rval, 5 - rval}, {Nbox, dim}], {0, 0, 0}];
 r = ConstantArray[rval, Length[p]];
 Sort[Sort /@ FindCollisionCandidates[{p - r, p + r}\[Transpose], 2, False]] // Timing
 Sort[Sort /@ FindCollisionCandidates[{p - r, p + r}\[Transpose], 2, True]] // Timing

yielding

{0.374402, {{4, 133}, {6, 199}, {7, 42}, {7, 52}, {7, 132}, {9, 98}, {10, 229}, {14, 117}, {15, 44}, {15, 100}, {15, 226}, {16, 49}, {21, 231}, {22, 237}, {26, 157}, {27, 207}, {34, 173}, {34, 185}, {35, 89}, {35, 119}, {41, 242}, {42, 52}, {42, 132}, {43, 154}, {44, 100}, {44, 226}, {46, 248}, {48, 245}, {50, 223}, {52, 132}, {53, 195}, {54, 194}, {55, 242}, {61, 94}, {67, 170}, {71, 155}, {71, 179}, {71, 216}, {74, 200}, {75, 141}, {81, 114}, {87, 159}, {89, 119}, {92, 131}, {104, 227}, {105, 225}, {112, 235}, {129, 176}, {142, 188}, {142, 222}, {144, 204}, {144, 210}, {146, 236}, {148, 207}, {164, 168}, {173, 185}, {179, 216}, {180, 210}, {183, 199}, {197, 221}, {214, 224}, {228, 244}}}

for the brute-force approach, and

{3.681624, {{4, 133}, {6, 199}, {7, 42}, {7, 52}, {7, 132}, {9, 98}, {10, 229}, {14, 117}, {15, 44}, {15, 100}, {15, 226}, {16, 49}, {21, 231}, {22, 237}, {26, 157}, {27, 207}, {34, 173}, {34, 185}, {35, 89}, {35, 119}, {41, 242}, {42, 52}, {42, 132}, {43, 154}, {44, 100}, {44, 226}, {46, 248}, {48, 245}, {50, 223}, {52, 132}, {53, 195}, {54, 194}, {55, 242}, {61, 94}, {67, 170}, {71, 155}, {71, 179}, {71, 216}, {74, 200}, {75, 141}, {81, 114}, {87, 159}, {89, 119}, {92, 131}, {104, 227}, {105, 225}, {112, 235}, {129, 176}, {142, 188}, {142, 222}, {144, 204}, {144, 210}, {146, 236}, {148, 207}, {164, 168}, {173, 185}, {179, 216}, {180, 210}, {183, 199}, {197, 221}, {214, 224}, {228, 244}}}

for the SAP approach.

This demonstrates a rather poor performance for the SAP implementation of almost 10x slower on my machine.

Now my question(s):

  1. Has anyone implemented a SAP (or similiar) collision detection algorithm using on-board Mathematica commands and in vectorized form?
  2. Is anyone aware of a potential "undocumented" feature for collision detection (e.g. for 2D or 3D graphics in Mathematica?
  3. Of course any suggestion to speed up the SAP code by vectorizing it, or extending it with temporal coherence exploitation are welcome too.
$\endgroup$
  • 1
    $\begingroup$ If your objects have reasonably small size bounds in all directions then you can almost certainly do better using Nearest. As for the implementation above, if you have many objects then Append and Delete will be slow. Might be better to emulate a linked list by nesting e.g. active = {active,newelement}. $\endgroup$ – Daniel Lichtblau Oct 7 '14 at 22:27
8
$\begingroup$

Edit: 4x faster after refactoring and using {# - #2, +##} &[+##, Abs[# - #2]] & instead of Sort/@.

Let we have the following 1D intervals

$$ [a_1,b_1], [a_2,b_2], [a_3, b_3]. $$

Then we can take all these edges of intervals and sort them. For example,

                                                             enter image description here

Now we see that the first interval intersects with all intervals, which have edges between $a_1$ and $b_1$.

We can easily implement it with Ordering

sap1[obj_] := Module[{ord = Quotient[Ordering@Flatten@obj + 1, 2]}, 
  Quotient[#, 2] &@ DeleteDuplicates[
    Transpose@{# - #2, +##} &[+##, Abs[# - #2]] &[
       Flatten[Unitize@# Range@Length@#], Flatten@#] &[
     ord[[# + 1 ;; #2 - 1]] & @@@ Partition[Ordering@ord, 2]]]]

sap1@objects // Sort
(* {{1, 3}, {1, 5}, {1, 6}, {1, 7}, {1, 8}, {1, 9}, {2, 6}, {3,
   5}, {3, 7}, {3, 8}, {3, 9}, {4, 6}, {5, 7}, {5, 9}, {7, 8}, {7, 
  9}, {7, 10}} *)

For multidimensional collision detecting you can use

Intersection @@Table[sap1@objboxes[[;; , ;; , d]], {d, dim}]

Timings:

FindCollisionCandidates[objboxes_, dim_: 3, sap_: True] := 
 If[sap, Intersection @@Table[sap1@objboxes[[;; , ;; , d]], {d, dim}], 
  Extract[Subsets[Range[Length[objboxes]], {2}], 
   Position[Transpose[Apply[IntervalIntersection, 
      Map[Interval, Subsets[Transpose[
            Transpose[objboxes, {1, 3, 2}]][[#]], {2}] & /@ 
        Range[dim], {3}], {2}]], {Repeated[Interval[_List], {dim}]}]]]

    rval = 0.1;
Nbox = 250;
dim = 3;
SeedRandom[35];
p = Prepend[RandomReal[{-5 + rval, 5 - rval}, {Nbox, dim}], {0, 0, 0}];
r = ConstantArray[rval, Length[p]];
res1 = Sort[
    Sort /@ FindCollisionCandidates[{p - r, p + r}\[Transpose], 2, 
      False]]; // AbsoluteTiming
res2 = FindCollisionCandidates[{p - r, p + r}\[Transpose], 2, 
    True]; // AbsoluteTiming
res1 == res2
(* {0.577945, Null} *)
(* {0.008889, Null} *)
(* True *)

It is definitely faster!

$\endgroup$
  • $\begingroup$ Great answer. Thank you very much! This was exactly the kind of input I was looking for... $\endgroup$ – Rainer Oct 8 '14 at 13:16

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.