OK. Been running both approaches within the c++ 'vectors.h' framework of 
POV-Ray. Initial idea is both forms could be in one inbuilt function and 
therefore "munction" too.

// Original macro approach:

#macro VAngle(V1, V2)
     acos(max(-1,min(1,vdot(vnormalize(V1),vnormalize(V2)))))
#end


// Tor Olav's proposal for better accuracy:
#macro AngleBetweenVectors(v1, v2)

     #local v1n = vnormalize(v1);
     #local v2n = vnormalize(v2);

     (2*atan2(vlength(v1n - v2n), vlength(v1n + v2n)))

#end // macro AngleBetweenVectors

-----------------------
Over millions of randomly generated vectors with component values in the 
+-25.0 range.

Max dot value seen: +0.99999999991558441
Min dot value seen: -0.99999999134553685

Largest difference between approaches was 0.00000000000995055 (1e-11).
The vectors for this difference were:

<6.67080701406617393,
9.26765558974242509,
20.94262630955438453>

and

<7.10376197570178647,
9.86894333496174880,
22.30091129769613900>

Angles near zero and pi (near same-direction and nearly opposing) 
showing the larger differences. The 0 angle case running larger for 
differences than those near pi.

The performance difference in the context of all programming framework
VAngle to AngleBetweenVectors.

127.15 -> 164.98 ---> +29.75%  slower

Taking a specific measure of the overhead out of both times:

100.5  -> 138.33 ---> +37.64%  slower

I have still a few other tests in mind(1), but busy for a couple days - 
starting now.

(1) - Different optimizations, etc.

Bill P.

Post a reply to this message

On 5/11/21 4:50 AM, William F Pokorny wrote:
> OK. Been running both approaches within the c++ 'vectors.h' framework of 
> POV-Ray. Initial idea is both forms could be in one inbuilt function and 
> therefore "munction" too.
> 
> // Original macro approach:
> 
> #macro VAngle(V1, V2)
>      acos(max(-1,min(1,vdot(vnormalize(V1),vnormalize(V2)))))
> #end
> 
> 
> // Tor Olav's proposal for better accuracy:
> #macro AngleBetweenVectors(v1, v2)
> 
>      #local v1n = vnormalize(v1);
>      #local v2n = vnormalize(v2);
> 
>      (2*atan2(vlength(v1n - v2n), vlength(v1n + v2n)))
> 
> #end // macro AngleBetweenVectors
> 
...
> 
> I have still a few other tests in mind(1), but busy for a couple days - 
> starting now.
> 
> (1) - Different optimizations, etc.
...

All more complicated than expected and I'm not normally concise, but 
I'll try...

I now have an inbuilt f_vangle() which supports both the old and Tor 
Olav's method - and where both are faster when used within VAngle making 
it now a "munction." I'll post more about f_vangle later in povray 
unofficial patches.

---
Using more targeted random vectors near the noisy parts:

Largest diff 2 ver    0.00000003650005986  ~ 1e-7
Max dot prod sum     +1.00000000000000067  (3x DBL_EPSILON)
Min dot prod sum     -1.00000000000000067
Min abs dot prod sum  0.00000000000000000

---
With the results above, one might ask how did we not have domain errors 
on the -1 side with all those acos(dot(... based macros we have?

Ans: It turns out the parser version of acos tests the inputs, clamps if 
need be and issues a warning on each bad-domain call. This means any 
math.inc clamping is really there to suppress the warning messages.

Why not a parse error over a warning is a good question - given warnings 
can be ignored. It would be better for the user to learn and correct code.

---
The raw acos result and so the vm version of acos returns nan (not a 
number) for bad inputs. Peripherally related compiler optimizations 
related to the -ffastmath optimizations cause internal routines like
POV_ISINF, POV_ISNAN (std::isnan and the like too) to always return 
false. And, POV_ISFINITE and IsFinite() to always return true. These are 
not much used, but this caught me up as I was trying to use them to pick 
up domain errors initially.

---
Related to this VAngle look we have in math.inc a macro called 
VCos_Angle which I believe should be clamped and it now is in my version.

// Cosine of angle between V1 and V2
// #macro VCos_Angle(V1, V2) vdot(vnormalize(V1), vnormalize(V2)) #end
    // +14% over million calls to clamp, but we should clamp.
    #macro VCos_Angle(V1, V2)
        max(-1,min(1,vdot(vnormalize(V1),vnormalize(V2))))
    #end

---
By experiment only it looks like the old method gives us about 7 
significant digits, Tor Olav's 12 significant digits. This is at 
extremes, however.

---
I ran both approaches posted at top against over a million random vector 
calls to measure the performance impact of the more accurate method.

One million calls to VAngle

9.65user 5.52system 0:15.70elapsed  <-- Original acos(dot..)macro.
10.01user 5.24system 0:15.79elapsed     with both min and max clip.
9.69user 5.35system 0:15.58elapsed
9.69user 5.36system 0:15.58elapsed
   "9.65 9.69 9.69"

One million calls to AngleBetweenVectors

11.21user 5.31system 0:17.05elapsed
11.35user 5.26system 0:17.14elapsed
11.57user 5.48system 0:17.57elapsed
11.43user 5.45system 0:17.42elapsed
   "11.21 11.35 11.43"                 29.03 -> 33.99 ---> +17.09%

Expect similar performance differences in any v3.8 based branch.

Bill P.

Post a reply to this message

Am 13.05.2021 um 21:35 schrieb William F Pokorny:

> Ans: It turns out the parser version of acos tests the inputs, clamps if 
> need be and issues a warning on each bad-domain call. This means any 
> math.inc clamping is really there to suppress the warning messages.
> 
> Why not a parse error over a warning is a good question - given warnings 
> can be ignored. It would be better for the user to learn and correct code.

I don't entirely agree.

There is a long-standing tradition that POV-Ray does not abandon parsing 
of a scene file if there's a sensible way to deal with the situation.

There are plenty of conceivable cases where the result of some 
computation should mathematically be within valid range for `acos`, but 
due to a numerical precision issue is actually not. Leaving it up to the 
user to protect against such issues would effectively mean that they'd 
have to safeguard pretty much each and every invocation of `acos` in 
such a way, which costs tons of time when done in POV-Ray's scene 
description language, compared to hard-coding such a mechanism into the 
`acos` function itself.

It could be argued that it would make sense to issue an error instead of 
a warning if the argument is out of range by more than just a bit. But 
for just barely out of range arguments, issuing a warning is the maximum 
POV-Ray should do about it.


And then there's also the long-standing tradition of not breaking 
existing scenes unless it can't be avoided. If we change the warning 
into a hard error, I guess there _will_ be legacy scenes broken. If we 
do this even for "minor infractions", breakage of legacy scenes is very 
much guaranteed.

Post a reply to this message