# Very fast way to do a coordinate frame transform

I need a function that rotates and translates a huge amount of line segments.

For example, I have a set of line segments in the form {{x0,y0,z0},{x1,y1,z1}}

data = RandomReal[{-1, 1}, {20000, 2, 3}]


and three vectors

n = {n1,n2,n3}
s = {s1,s2,s3} (* Orthogonal to n *)
p = {p1,p2,p3}


I want to perform a rigid (right handed) transformation that sends p to the origin, n to the $z$ axis and $s$ to the $y$ axis, in order to apply it to data.

I know about RotationMatrix, RotationTransform, and I have coded a (very ugly) solution but, due to inexperience and the amount of data, I'm struggling with performance.

I would be very grateful if someone can provide a fast solution and, if possible, to explain why is fast.

EDIT

Here is the ugly code, for my own embarrassment

rotZ = RotationTransform[{n, {0, 0, 1}}];
rotY = RotationTransform[{s, {0, 1, 0}}];
RT[v_]:= rotY[rotZ[-p + #]] & /@ v;


and then

dataTrans = RT /@ data


Taking

n = RandomReal[{-1, 1}, 3];
s = {-n[[2]], n[[1]], 0}; (* n and s are orthogonal *)
p = RandomReal[{-1, 1}, 3];


then

Timing[RT /@ data][[1]]
(*    3.60422    *)


In an Intel Core2 Duo CPU T8100 @2.10GHz, 4gb ram Ubuntu 12.04, Mathematica 8 distribution.

• It might help to post the code, "ugly" or not. Aug 11, 2012 at 21:36
• Are n,s,p already orthogonal? Aug 11, 2012 at 22:00
• n and s are. Aug 11, 2012 at 22:13
• There's an ampersand missing in the definition of RT Aug 11, 2012 at 22:23
• @SjoerdC.deVries Thanks for all the debugging! Aug 11, 2012 at 22:43

Defining n, s and p with their numeric values before defining the transformation function helps considerably, but a huge difference is made by changing the definition of RT to use the numerical matrices instead of the RotationTransform-s:

data = RandomReal[{-1, 1}, {20000, 2, 3}];
n = RandomReal[{-1, 1}, 3];
s = {-n[[2]], n[[1]], 0};(*n and s are orthogonal*)
p = RandomReal[{-1, 1}, 3];
rotZ = RotationTransform[{n, {0, 0, 1}}];
rotY = RotationTransform[{s, {0, 1, 0}}];

m = TransformationMatrix[rotY][[1 ;; 3, 1 ;; 3]].TransformationMatrix[
rotZ][[1 ;; 3, 1 ;; 3]];

RT1[v_] := m.(-p + #) & /@ v;
RT2[v_] := rotY[rotZ[-p + #]] & /@ v;

Timing[r1 = RT1 /@ data][[1]]
Timing[r2 = RT2 /@ data][[1]]


0.328

2.824

• I didn't know that one could extract the numerical matrix from RotationTransform. Is there any difference (qualitatively) than simply use RotationMatrix? Aug 12, 2012 at 4:08
• Let RT3[segs_] := Transpose[m.(Transpose[segs, {3, 2, 1}] - p), {3, 2, 1}] and data = RandomReal[{-1, 1}, {200000, 2, 3}]. The timing for RT1 /@ data is 3.156 seconds; for RT3[data] (using efficient multithreaded Transpose and Dot) it is a matter of milliseconds. Aug 12, 2012 at 4:31
• @OleksandrR. I suspected the (-p + #) & /@ approach was inefficient, but I didn't know how else to do the sum. Can you explain why the program can do the sum when dimensions are {3,2,N} but it can't when they are {N,2,3}? Aug 12, 2012 at 5:14
• @Manuel that's just the way Mathematica is--which is my way of saying I don't know. Personally, I'd find exhaustive array broadcasting rules (like one finds in NumPy) very useful, but for whatever reason there hasn't been much effort in this direction. There are only so many features WRI can work on at once, I suppose, and Transpose usually gets the job done after all. Aug 12, 2012 at 5:21