RE: Re: Generating Two Unit Orthogonal Vectors to a 3D Vector
- To: mathgroup at smc.vnet.net
- Subject: [mg36413] RE: [mg36377] Re: [mg36352] Generating Two Unit Orthogonal Vectors to a 3D Vector
- From: "DrBob" <drbob at bigfoot.com>
- Date: Fri, 6 Sep 2002 03:16:42 -0400 (EDT)
- Reply-to: <drbob at bigfoot.com>
- Sender: owner-wri-mathgroup at wolfram.com
I timed Daniel's three solutions and Gary's one (plus a couple of my own
a little later):
perps1[v_] := If[v[[1]] ==
v[[2]] == 0, {{1, 0, 0}, {0, 1, 0}}, {{v[[2]], -v[[
1]], 0}, Cross[v, {v[[2]], -v[[1]], 0}]}]
perps2[v_] := With[{vecs = NullSpace[{v}]}, {vecs[[
1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]}]
perps2C = Compile[{{v, _Real, 1}}, Module[{vecs = NullSpace[{v}]}, {
vecs[[1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]}]]
helzer[v : {a1_, a2_,
a3_}] := With[{w = First[Sort[{{a2, -a1, 0}, {a3,
0, -a1}, {0, a3, -a2}}, OrderedQ[{Plus @@ Abs[#2], Plus @@ \
Abs[#1]}] &]]}, {w, Cross[v, w]}]
vecs = Table[Random[], {10000}, {3}];
Timing[perps1 /@ vecs; ]
Timing[perps2 /@ vecs; ]
Timing[perps2c /@ vecs; ]
Timing[helzer /@ vecs; ]
{1.7350000000000012*Second, Null}
{0.5619999999999994*Second, Null}
{0.219 Second, Null}
{2.7349999999999994*Second, Null}
I made a small change to Daniel's perps1, and got a solution as fast as
perps2c, WITHOUT compiling. Compiling tripled the speed again, so
"treatC" is the fastest solution I've seen so far.
treat[{a_, b_, c_}] :=
If[a == b == 0, {{1, 0, 0}, {0, 1, 0}}, {{b, -a, 0}, {a*c, b*c, -a^2 -
b^2}}]
treatC = Compile[{{v, _Real, 1}},
If[v[[1]] == v[[2]] == 0,
{{1, 0, 0}, {0, 1, 0}},
{{v[[2]], -v[[1]], 0},
{v[[1]]*v[[3]], v[[2]]*v[[3]],
-v[[1]]^2 - v[[2]]^2}}]]
vecs = Table[Random[], {10000}, {3}];
Timing[perps2c /@ vecs; ]
Timing[helzer /@ vecs; ]
Timing[treat /@ vecs; ]
Timing[treatC /@ vecs;]
{0.2190000000000083*Second, Null}
{2.7339999999999947*Second, Null}
{0.25*Second, Null}
{0.07800000000000296*Second, Null}
None of these solutions reliably return normalized vectors.
Bobby Treat
-----Original Message-----
From: Daniel Lichtblau [mailto:danl at wolfram.com]
To: mathgroup at smc.vnet.net
Subject: [mg36413] [mg36377] Re: [mg36352] Generating Two Unit Orthogonal Vectors
to a 3D Vector
David Park wrote:
>
> There are many cases in graphics, and otherwise, where it is useful to
> obtain two orthogonal unit vectors to a given vector. I know a number
of
> ways to do it, but they all seem to be slightly inelegant. I thought I
would
> pose the problem to MathGroup. Who has the most elegant Mathematica
> routine...
>
> OrthogonalUnitVectors::usage = "OrthogonalUnitVectors[v:{_,_,_}] will
return
> two unit vectors orthogonal to each other and to v."
>
> You can assume that v is nonzero.
>
> David Park
> djmp at earthlink.net
> http://home.earthlink.net/~djmp/
Some possibilities:
perps1[v_] := If [v[[1]]==v[[2]]==0,
{{1,0,0},{0,1,0}},
{{v[[2]],-v[[1]],0}, Cross[v,{v[[2]],-v[[1]],0}]}
]
perps2[v_] := With[{vecs=NullSpace[{v}]},
{vecs[[1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]}
]
This appears to be 2-3 times faster than perps1 for vectors of machine
reals. I get another factor of 2 using Compile, which is appropriate for
e.g. graphics use.
perps2C = Compile[{{v,_Real,1}},
Module[{vecs=NullSpace[{v}]},
{vecs[[1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]}
]]
In[61]:= vecs = Table[Random[], {10000}, {3}];
In[62]:= Timing[p2 = Map[perps1C,vecs];]
Out[62]= {0.49 Second, Null}
This is on a 1.5 GHz processor.
Daniel Lichtblau
Wolfram Research