Re: pfDotVec3 in combination with pfArcCos

New Message Reply Date view Thread view Subject view Author view

Don Hatch (hatch++at++hell.engr.sgi.com)
Sat, 22 May 1999 01:04:10 -0700

On Mar 22, 6:57pm, MLM Veraart wrote:
> Subject: pfDotVec3 in combination with pfArcCos
> Hello,
>
> It has taken me some time to find out why in some situations my
> performer application shoes some very wierd behaviour.
> The application use shadowmap-shadows, and it calculates a tight
> frustrum for the shadow lightsource. For this frustum calculation
> I have to determine a lot of angles between vectors.
> To to this I use the following code snippet
>
> pfVec3 v1,v2;
> float angle;
>
> pfNormaliseVec3(v1);
> pfNormaliseVec3(v2);
> angle = pfArcCos(pfDotVec3(v1,v2));
>
> it turns out that I sometime get a result angle of nan(0x2000000).
> And any further calculation with this gave wrong results.
> The cause of it was that pfDotVec gave as a result a float
> that was just larger than 1.0, in only a few cases.
>
> Solution:
> write your own arccos function that checks for boundaries
> and then calls pfArcCos.
>
> Mario
>-- End of excerpt from MLM Veraart

That works, but keep in mind that taking the arccosine of
a dot product is inherently unstable/lossy near angles of 0 and pi,
and the catastrophic result you encountered can be considered the extreme
case of this general problem which remains.
Roughly speaking, that method destroys about half the significant bits,
so if you're using doubles you end up with only float accuracy,
and if you're using floats you end up with only
about 3 or 4 digits of accuracy (including the zeros after the decimal point,
in this case).
This may be good enough for your application,
but if not, you may want to try something like this instead:

    pfNormalizeVec3(v1);
    pfNormalizeVec3(v2);
    float dotprod = pfDotVec3(v1,v2);
    if (dotprod*dotprod <= .5f) // then acos is stable
        angle = pfArcCos(dotprod);
    else // asin is stable
    {
        pfVec3 crossprod;
        pfCrossVec3(crossprod,v1,v2);
        if (dotprod >= 0)
            angle = pfArcSin(pfLengthVec3(crossprod));
        else
            angle = PF_PI - pfArcSin(pfLengthVec3(crossprod));
    }

or:

    pfNormalizeVec3(v1);
    pfNormalizeVec3(v2);
    if (pfDotVec3(v1,v2) >= 0)
        angle = 2 * pfArcSin(pfDistance3(v1,v2) * .5f);
    else
    {
        pfVec3 minus_v1;
        pfNegateVec3(minus_v1,v1);
        angle = PF_PI - 2 * pfArcSin(pfDistance3(minus_v1,v2) * .5f);
    }

(thanks to Lee Willis for the latter idea).

Another (less inspired) solution would be to calculate the normalization,
dot product and arccosine (with bounds check) in double precision
to get a single-precision result,
but then you'd have to switch to a completely different method if you ever
wanted an accurate double precision result.

Don

-- 
Don Hatch  hatch++at++sgi.com  (650) 933-5150  Silicon Graphics, Inc.

New Message Reply Date view Thread view Subject view Author view

This archive was generated by hypermail 2.0b2 on Sat May 22 1999 - 01:04:17 PDT

This message has been cleansed for anti-spam protection. Replace '++at++' in any mail addresses with the '@' symbol.