jcs@crash.cts.com (John Schultz) (12/31/90)
Given two 3D vectors, compute a third vector normal to the first
and in the plane of both. The two given vectors are unit length, and
the result must be unit length. I've come up with two solutions:
Given unit vectors A and B, result is C:
1. C = normalize(A cross (A cross B)), [15 muls, 8 add/subs, 1 sqrt,
and 3 divs].
2. C = (A*(A dot B) - B) / (sin(acos(A dot B))), [6 muls, 5 add/subs,
2 table lookups,
Here's what it looks like: and 3 divs].
C
| B
| /
| /
| /
|/_____A
Method two is much faster, and produces correct results with error
increasing as the angle between A and B decreases. I've tried both
methods in double precision floating point and 32/64 bit integer, and
get similar results. Method one has much less error. Is there a faster
way to do this, or with less error?
Johnlevine@well.sf.ca.us (Ron Levine) (01/01/91)
Newsgroups: comp.graphics Subject: Re: 2 3D Vectors: Get a 3rd perp. to 1st? Summary: Expires: References: <6613@crash.cts.com> Sender: Followup-To: Distribution: Organization: Whole Earth 'Lectronic Link, Sausalito, CA Keywords: In article <6613@crash.cts.com> jcs@crash.cts.com (John Schultz) writes: > Given two 3D vectors, compute a third vector normal to the first >and in the plane of both. The two given vectors are unit length, and >the result must be unit length. I've come up with two solutions: > > Given unit vectors A and B, result is C: > > 1. C = normalize(A cross (A cross B)), [15 muls, 8 add/subs, 1 sqrt, > and 3 divs]. > > 2. C = (A*(A dot B) - B) / (sin(acos(A dot B))), [6 muls, 5 add/subs, > 2 table lookups, > Method two is much faster, and produces correct results with error >increasing as the angle between A and B decreases. ... >... Method one has much less error. Is there a faster >way to do this, or with less error? > > The source of your errors is probably in the table look ups for the trig functions. Why not just use C = normalize(B - (A dot B)A) It requires 9 multiplies, 7 adds, 1 sqrt and 3 divides. This is well-known as the "Gram-Schmidt procedure".
jcs@crash.cts.com (John Schultz) (01/02/91)
In <22365@well.sf.ca.us> levine@well.sf.ca.us (Ron Levine) writes: >In article <6613@crash.cts.com> jcs@crash.cts.com (John Schultz) writes: >> Given two 3D vectors, compute a third vector normal to the first >>and in the plane of both. The two given vectors are unit length, and >>the result must be unit length. I've come up with two solutions: >> >> Given unit vectors A and B, result is C: >> >> 1. C = normalize(A cross (A cross B)), [15 muls, 8 add/subs, 1 sqrt, >> and 3 divs]. >> >> 2. C = (A*(A dot B) - B) / (sin(acos(A dot B))), [6 muls, 5 add/subs, >> 2 table lookups, >> Method two is much faster, and produces correct results with error >>increasing as the angle between A and B decreases. ... >>... Method one has much less error. Is there a faster >>way to do this, or with less error? >The source of your errors is probably in the table look ups >for the trig functions. Why not just use Nope. Tried it using the 68882 and still got large errors, although slightly lower than the tabled version (but not much). Double precision floats were used. When drawn on paper, you can see that as A dot B approaches 1.0, the length of the new vector (C) becomes shorter, and sin(acos(A dot B)) grows smaller. It appear that the errors are due to the accumulation of error of the short vector being scaled by a length derived using transcendentals, which only has so many bits of accuracy. > C = normalize(B - (A dot B)A) I tried that as well, before deriving the non-sqrt version. On a machine where muls.l is 44 cycles and sqrt is (60+107+100) 267 cycles, sqrt is to be avoided! I just now looked up Gram-Schmidt (process), of which the above is a small part. There should be a book titled "Well Known Procedures for Computer Graphics", subtitled "so you don't have to re-invent the wheel". (Kind of like Graphics Gems, only all procedures. All of them known to man at the time of printing. A phone book of tools, with an intelligent index making the knowledge the reader seeks easy to find). But, deriving the stuff is fun... >It requires 9 multiplies, 7 adds, 1 sqrt and 3 divides. 3 extra muls, 2 extra adds, and 1 sqrt. Thats about 415 extra cycles, which is too expensive (my work is real-time). >This is well-known as the "Gram-Schmidt procedure". Yup, I figured what I derived was probably well-known... But they didn't use the trig-projection trick (which is probably obvious). The best solution then, looks like a two part heuristic: When A dot B is large, use the normalization method, else use the much faster table lookup method. John
ccplumb@rose.uwaterloo.ca (Colin Plumb) (01/16/91)
jcs@crash.cts.com (John Schultz) wrote: >In <22365@well.sf.ca.us> levine@well.sf.ca.us (Ron Levine) writes: >>In article <6613@crash.cts.com> jcs@crash.cts.com (John Schultz) writes: >>> Given two 3D vectors, compute a third vector normal to the first >>>and in the plane of both. The two given vectors are unit length, and >>>the result must be unit length. I've come up with two solutions: >>> >>> Given unit vectors A and B, result is C: >>> >>> 2. C = (A*(A dot B) - B) / (sin(acos(A dot B))), [6 muls, 5 add/subs, >>> 2 table lookups, >>> Method two is much faster, and produces correct results with error >>>increasing as the angle between A and B decreases. ... >>>... Method one has much less error. Is there a faster >>>way to do this, or with less error? Why two table lookups? sin(acos(x)) = sqrt(1-x^2) and the division can be done in one lookup or, if you prefer, near x=0, 1/sqrt(1-x^2) = 2 4 6 35 8 63 10 231 12 429 14 1 + 1/2 x + 3/8 x + 5/16 x + --- x + --- x + ---- x + ---- x 128 256 1024 2048 6435 16 12155 18 46189 20 88179 22 676039 24 + ----- x + ----- x + ------ x + ------ x + ------- x 32768 65536 262144 524288 4194304 1300075 26 5014575 28 30 + ------- x + -------- x + O(x ) 8388608 33554432 This is similar to a problem I'm trying to solve. I want to do very fast (real-time animation) incremental rotation of an orthonormal basis in 3-space. E.g. rotate 1/20 degree up, then 1/20 degree right, then 1/20 degree right, etc. The user has a joystick, and is flying over a surface I'm texture-mapping. The step size is always the same (1/20 degree is about right) but just multiplying by the matrix doesn't work if you're composing it 10,000+ times. I have to correct so the vectors stay orthonormal. In my application, accuracy is not important, as long as I don't get sudden jumps in any of the vectors, and they never drift too far from perpendicular. Does anyone know a recurrence that takes an almost orthonormal basis and returns a more nearly orthonormal basis that I can execute once per step or something? (If it matters, I'm using scaled integers, not floating point.) Your help is greatly appreciated. -- -Colin
falk@peregrine.Sun.COM (Ed Falk) (01/19/91)
In article <7012@crash.cts.com> jcs@crash.cts.com (John Schultz) writes: >In <1991Jan15.233218.7646@watdragon.waterloo.edu> ccplumb@rose.uwaterloo.ca (Colin Plumb) writes: > > Yuck! First, I do two table lookups because I have sin and asin tabled >already for all possible cases. Can't do that for sqrt... Sure you can. Just have a table of all your sqrts from sqrt(1) to sqrt(4) at any resolution you want. Then: i = 0 ; while(x < 1.) x *= 4. , ++i ; while(x > 4.) x *= .25 , --i ; x = sqrt_table_lookup(x) ; while( i > 0 ) x *= .5 ; while( i < 0 ) x *= 2. ; (Optimization left as an excercise to the student. A smart compiler should be able to convert the multiplies into increments or decrements of the exponent.) -ed falk, sun microsystems sun!falk, falk@sun.com To be loyal to rags, to shout for rags, to worship rags, to die for rags -- that is a loyalty of unreason, it is pure animal (Mark Twain).