Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
3  
It might help to post the code, "ugly" or not. –  Verbeia Aug 11 '12 at 21:36
    
Are n,s,p already orthogonal? –  belisarius Aug 11 '12 at 22:00
    
n and s are. –  Pragabhava Aug 11 '12 at 22:13
1  
There's an ampersand missing in the definition of RT –  Sjoerd C. de Vries Aug 11 '12 at 22:23
1  
@SjoerdC.deVries Thanks for all the debugging! –  Pragabhava Aug 11 '12 at 22:43
show 1 more comment

1 Answer 1

up vote 7 down vote accepted

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

Mathematica graphics

share|improve this answer
    
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
3  
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
1  
@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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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