Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2002
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2002

[Date Index] [Thread Index] [Author Index]

Search the Archive

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/



  • Prev by Date: Re: approximation for partitial binomial sum
  • Next by Date: Why this simply elliptic integral fails
  • Previous by thread: Trouble with Export
  • Next by thread: Why this simply elliptic integral fails