Re: Generating Two Unit Orthogonal Vectors to a 3D Vector - Summary
- To: mathgroup at smc.vnet.net
- Subject: [mg36399] Re: [mg36352] Generating Two Unit Orthogonal Vectors to a 3D Vector - Summary
- From: "David Park" <djmp at earthlink.net>
- Date: Wed, 4 Sep 2002 21:22:12 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Thanks for the many replies.
My main application for this is in drawing a piece of a plane at a point
normal to a vector. For me, efficiency is not an important criterion. But I
want orthogonal vectors so the plane will be composed of squares and not
diamonds. I want unit vectors so that I can easily specify the size of the
piece I am drawing.
I posted a summary last evening and it doesn't seem to have appeared today.
In any case, I have realized that none of the routines really require the
LinearAlgebra`Orthogonalization package, and run faster without it. I have
modified the original routines to use a pure function instead of Normalize
or GramSchmidt.
Here is a summary of the responses so far. The timings were done with the
following exact vector. They are faster with a machine precision vector.
vtest = {1, 2, 1};
__________________________________________________________
Hugh Goyder...
The original routine was
OrthogonalUnitVectors[v : {_, _, _}] := GramSchmidt[NullSpace[{v}]]
but I changed it to
OrthogonalUnitVectors[v : {_, _, _}] := #/Sqrt[#.#] & /@ NullSpace[{v}]
OrthogonalUnitVectors[vtest]
{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {-(1/Sqrt[3]), 1/Sqrt[3],
-(1/Sqrt[3])}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{2.53 Second, Null}
This appears to be the shortest and fastest routine. I would give it the
prize for elegance. It is certainly easy to remember.
________________________________________________________
Ingolf Dahl...
OrthogonalUnitVectors[v:{_, _, _}] := (#1/Sqrt[#1 . #1] & ) /@
NullSpace[{v, {0, 0, 0}, {0, 0, 0}}]
OrthogonalUnitVectors[vtest]
{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {-(2/Sqrt[5]), 1/Sqrt[5], 0}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{2.59 Second, Null}
This is really the same as Hugh Goyder's routine. It isn't necessary to add
the zero rows to the matrix. NullSpace automatically does it.
________________________________________________________
Selwyn Hollis...
OrthogonalUnitVectors[v : {_, _, _}] :=
Module[{u, w},
u = Which[(w = {0, v[[3]], -v[[2]]}).w != 0, w,
(w = {v[[3]], 0, -v[[1]]}).w != 0, w,
(w = {v[[2]], -v[[1]], 0}).w != 0, w,
True, Return[$Failed]];
#/Sqrt[#.#] & /@ {u, Cross[u, v]} ]
OrthogonalUnitVectors[vtest]
{{0, 1/Sqrt[5], -(2/Sqrt[5])}, {Sqrt[5/6], -Sqrt[2/15],
-(1/Sqrt[30])}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{11.92 Second, Null}
I wish I understood the reason for using rows of the v Cross matrix as test
vectors. It does have the advantage that a dot product is used to give a
scalar test for zero. Why not...
OrthogonalUnitVectors[v : {_, _, _}] :=
Module[{u, w},
u = Which[(w = {1, 0, 0})\[Cross]v != {0, 0, 0}, w,
(w = {0, 1, 0})\[Cross]v != {0, 0, 0}, w,
True, Return[$Failed]];
#/Sqrt[#.#] & /@ {u, Cross[u, v]} ]
OrthogonalUnitVectors[vtest]
{{1, 0, 0}, {0, -(1/Sqrt[5]), 2/Sqrt[5]}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{10.33 Second, Null}
____________________________________________________________
John Browne...
OrthogonalUnitVectors[v : {_, _, _}] :=
Module[{r, v1, v2}, r = {Random[], Random[], Random[]}; v1 = Cross[v, r];
v2 = Cross[v1, v]; {v1/Sqrt[Dot[v1, v1]], v2/Sqrt[Dot[v2, v2]]}]
OrthogonalUnitVectors[vtest]
{{0.26601, -0.534209, 0.802408}, {-0.873254, 0.218984, 0.435286}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{9.17 Second, Null}
This has the problem that the r vector might be parallel to v. But in 100000
test cases it never happened.
____________________________________________________________
Daniel Lichblau...
(Daniel left out the normalization, which I added.)
OrthogonalUnitVectors[v:{_, _, _}] :=
With[{vecs = NullSpace[{v}]}, (#1/Sqrt[#1 . #1] & ) /@
{vecs[[1]], vecs[[2]] - vecs[[2]] . vecs[[1]]*vecs[[1]]}]
OrthogonalUnitVectors[vtest]
{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {0, 1/Sqrt[5], -(2/Sqrt[5])}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{3.35 Second, Null}
______________________________________________________________
Garry Helzer
OrthogonalVectors[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]}
]
OrthogonalUnitVectors[vtest]
{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {-(2/Sqrt[5]), 1/Sqrt[5], 0}}
Do[OrthogonalUnitVectors[vtest], {10000}] // Timing
{2.64 Second, Null}
__________________________________________________________________
Tom Burton
"A nice feature of the following simple scheme is that you know that the two
points of discontinuity are where r and -r cross the unit sphere."
triad[u_, r_:{1, 0, 0}] := Module[{v = Cross[u, r]}, #/Sqrt[# . #] & /@
{u, v, Cross[u, v]}]
triad[vtest]
{{1/Sqrt[6], Sqrt[2/3], 1/Sqrt[6]},
{0, 1/Sqrt[5], -(2/Sqrt[5])}, {-Sqrt[5/6], Sqrt[2/15],
1/Sqrt[30]}}
Do[triad[vtest], {10000}] // Timing
{17.02 Second, Null}
David Park
djmp at earthlink.net
http://home.earthlink.net/~djmp/
From: David Park [mailto:djmp at earthlink.net]
To: mathgroup at smc.vnet.net
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/