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[], n[], 0}; (* n and s are orthogonal *)
p = RandomReal[{-1, 1}, 3];


then

Timing[RT /@ data][]
(*    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. – Verbeia Aug 11 '12 at 21:36
• Are n,s,p already orthogonal? – Dr. belisarius Aug 11 '12 at 22:00
• n and s are. – Pragabhava Aug 11 '12 at 22:13
• There's an ampersand missing in the definition of RT – Sjoerd C. de Vries Aug 11 '12 at 22:23
• @SjoerdC.deVries Thanks for all the debugging! – Pragabhava Aug 11 '12 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[], n[], 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][]
Timing[r2 = RT2 /@ data][]


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? – Pragabhava Aug 12 '12 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. – Oleksandr R. Aug 12 '12 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}? – Pragabhava Aug 12 '12 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. – Oleksandr R. Aug 12 '12 at 5:21