In creating some testing for include file macros I came across code for 
the two VAngle* macros (and two similar rotation related macros):

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

The min(1,..) is obviously trying to protect acos from domain errors.

Why only the clamp to the positive side? Aren't we as exposed to <-1.0?

Often enough in code - including POV_Ray's own - I've seen code clamping 
to a [-1..1] range after a dot product into acos() for acos domain 
concerns.

I cannot come up with reasoning to clamp only the positive side. Anyone 
else? If valid to do, we might be able to slightly streamline some 
internal code.

Bill P.

Post a reply to this message

Le 2021-05-06 à 09:06, William F Pokorny a écrit :
> In creating some testing for include file macros I came across code for 
> the two VAngle* macros (and two similar rotation related macros):
> 
> #macro VAngle(V1, V2)
>         acos(min(1, vdot(vnormalize(V1), vnormalize(V2))))
> #end
> 
> The min(1,..) is obviously trying to protect acos from domain errors.
> 
> Why only the clamp to the positive side? Aren't we as exposed to <-1.0?
> 
> Often enough in code - including POV_Ray's own - I've seen code clamping 
> to a [-1..1] range after a dot product into acos() for acos domain 
> concerns.
> 
> I cannot come up with reasoning to clamp only the positive side. Anyone 
> else? If valid to do, we might be able to slightly streamline some 
> internal code.
> 
> Bill P.
> 

The reason is quite simple : Normally, the vdot function of two 
normalized vectors can not return a value that is less than -1 nor more 
than 1.

So, the min() could probably be safely removed.

#macro VAngle(V1, V2)
	acos(vdot(vnormalize(V1), vnormalize(V2)))
#end

Post a reply to this message

William F Pokorny <ano### [at] anonymousorg> wrote:
> In creating some testing for include file macros I came across code for
> the two VAngle* macros (and two similar rotation related macros):
>
> #macro VAngle(V1, V2)
>          acos(min(1, vdot(vnormalize(V1), vnormalize(V2))))
> #end
>
> The min(1,..) is obviously trying to protect acos from domain errors.
>
> Why only the clamp to the positive side? Aren't we as exposed to <-1.0?
>
> Often enough in code - including POV_Ray's own - I've seen code clamping
> to a [-1..1] range after a dot product into acos() for acos domain
> concerns.
>
> I cannot come up with reasoning to clamp only the positive side. Anyone
> else? If valid to do, we might be able to slightly streamline some
> internal code.

Perhaps the one that made that macro experienced a problem with at value from
the dot product being slightly higher than one so he inserted that min()
function to deal with it, without thinking through the problem thoroughly.

You are right there should also be a max() function there to discard values
below -1.

Anyway the way that macro calculates the angle is not the best numerically.
I would rather suggest something like this:


#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 not tested this macro yet, but the angle() function in my scikit-vectors
Python library seems to work ok. It works in is similar way as the macro above.

https://github.com/t-o-k/scikit-vectors/blob/master/skvectors/cartesian_vectors.py



"Computing Cross-Products and Rotations in 2- and 3-Dimensional Euclidian
Spaces"
https://people.eecs.berkeley.edu/~wkahan/MathH110/Cross.pdf

My opinion is that several of the other math related macros in the include files
need some rework.

--
Tor Olav
http://subcube.com
https://github.com/t-o-k

Post a reply to this message

Alain Martel <kua### [at] videotronca> wrote:

> > In creating some testing for include file macros I came across code for
> > the two VAngle* macros (and two similar rotation related macros):
> >
> > #macro VAngle(V1, V2)

> > #end
> >
> > The min(1,..) is obviously trying to protect acos from domain errors.
> >
> > Why only the clamp to the positive side? Aren't we as exposed to <-1.0?
> >
> > Often enough in code - including POV_Ray's own - I've seen code clamping
> > to a [-1..1] range after a dot product into acos() for acos domain
> > concerns.
> >
> > I cannot come up with reasoning to clamp only the positive side. Anyone
> > else? If valid to do, we might be able to slightly streamline some
> > internal code.
> >
> > Bill P.
> >
>
> The reason is quite simple : Normally, the vdot function of two
> normalized vectors can not return a value that is less than -1 nor more
> than 1.

I think that this is only in theory. In practice there will sometimes be
roundoff errors that cause problems.


> So, the min() could probably be safely removed.
>
> #macro VAngle(V1, V2)
>  acos(vdot(vnormalize(V1), vnormalize(V2)))
> #end

--
Tor Olav
http://subcube.com
https://github.com/t-o-k

Post a reply to this message

On 5/6/21 10:58 AM, Tor Olav Kristensen wrote:
> William F Pokorny <ano### [at] anonymousorg> wrote:
>> In creating some testing for include file macros I came across code for
>> the two VAngle* macros (and two similar rotation related macros):
>>
>> #macro VAngle(V1, V2)
>>           acos(min(1, vdot(vnormalize(V1), vnormalize(V2))))
>> #end
>>
>> The min(1,..) is obviously trying to protect acos from domain errors.
>>
>> Why only the clamp to the positive side? Aren't we as exposed to <-1.0?
>>
>> Often enough in code - including POV_Ray's own - I've seen code clamping
>> to a [-1..1] range after a dot product into acos() for acos domain
>> concerns.
>>
>> I cannot come up with reasoning to clamp only the positive side. Anyone
>> else? If valid to do, we might be able to slightly streamline some
>> internal code.
> 
> Perhaps the one that made that macro experienced a problem with at value from
> the dot product being slightly higher than one so he inserted that min()
> function to deal with it, without thinking through the problem thoroughly.
> 
> You are right there should also be a max() function there to discard values
> below -1.
> 

Thanks. What I suspected too, but I don't have near your expertise.

> Anyway the way that macro calculates the angle is not the best numerically.
> I would rather suggest something like this:
> 
> 
> #macro AngleBetweenVectors(v1, v2)
> 
>      #local v1n = vnormalize(v1);
>      #local v2n = vnormalize(v2);
> 
>      (2*atan2(vlength(v1n - v2n), vlength(v1n + v2n)))
> 
> #end // macro AngleBetweenVectors
> 

Interesting, thanks. I'll play with it. There will be four sqrt calls 
with this vs two so I expect it's a slower approach.

Reminds me I captured fast approximation atan2 code I was thinking of 
making an inbuilt function, but I've not gotten to it...

And that thought reminds me of my plans to turn many macros into 
"munctions"(1) calling new faster, inbuilt functions. Done only a tiny 
bit of that so far.

(1) - Yes, thinking about creating a munctions.inc file breaking out 
what are really macro wrappers to functions into that file to make clear 
what they are. Maybe all as M_ prefixed so my current F_dents() would 
become M_dents(). Angle perhaps then M_Angle.

> 
> I have not tested this macro yet, but the angle() function in my scikit-vectors
> Python library seems to work ok. It works in is similar way as the macro above.
> 
> https://github.com/t-o-k/scikit-vectors/blob/master/skvectors/cartesian_vectors.py
> 
> 

> "Computing Cross-Products and Rotations in 2- and 3-Dimensional Euclidian
> Spaces"
> https://people.eecs.berkeley.edu/~wkahan/MathH110/Cross.pdf
> 
> My opinion is that several of the other math related macros in the include files
> need some rework.
> 

I'd consider alternatives for my branch if you want to offer details. 
Suppose I could compare to forms in skvectors?

Something I've done too with the inbuilt function changes is allow 
alternative calculations for the same thing. Compiled such conditionals 
not that expensive. Most frequently single float options where those are 
much faster, but different methods too. For one thing, the alternatives 
I can bang against each other as a way to test!

Bill P.

Post a reply to this message

William F Pokorny <ano### [at] anonymousorg> wrote:
> On 5/6/21 10:58 AM, Tor Olav Kristensen wrote:
> >...
> > My opinion is that several of the other math related macros in the include files
> > need some rework.
> >
>
> I'd consider alternatives for my branch if you want to offer details.

That sounds good.

Today I looked at the statistical macros in math.inc.

Here's how they could be rewritten:
(Please note that have not tested these thoroughly yet.)


// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#version 3.8;


// Does also work for vectors
#macro Mean(Values)

   #local N = dimension_size(Values, 1);
   #local Sum = Values[0];
   #for (I, 1, N - 1)
      #local Sum = Sum + Values[I];
   #end // for

   (Sum/N)

#end // macro Mean


// Does also work for vectors
#macro Variance(Values)

    #local N = dimension_size(Values, 1);
    #local M = Mean(Values);
    #local D = Values[0] - M;
    #local SquareSum = D*D;
    #for (I, 1, N - 1)
        #local D = Values[I] - M;
        #local SquareSum = SquareSum + D*D;
    #end

    (SquareSum/N)

#end // macro Variance


#macro Minimum(Values)

   #local N = dimension_size(Values, 1);
   #local Min = Values[0];
   #for (I, 1, N - 1)
      #local Min = min(Min, Values[I]);
   #end // for

   Min

#end // macro Minimum


#macro Maximum(Values)

   #local N = dimension_size(Values, 1);
   #local Max = Values[0];
   #for (I, 1, N - 1)
      #local Max = max(Max, Values[I]);
   #end // for

   Max

#end // macro Maximum


// For population
#macro StandardDeviation(Values)

   sqrt(Variance(Values))

#end // macro StandardDeviation


#macro Statistics(Values)

    #local N = dimension_size(Values, 1);
    #local Value = Values[0];
    #local Min = Value;
    #local Max = Value;
    #local Sum = Value;
    #local SquareSum = pow(Value, 2);
    #for (I, 1,  N - 1)
        #local Value = Values[I];
        #local Min = min(Min, Value);
        #local Max = max(Max, Value);
        #local Sum = Sum + Value;
        #local SquareSum = SquareSum + pow(Value, 2);
    #end // for
    #local M = Sum/N;
    #local Var = SquareSum/N - pow(M, 2);

    dictionary {
        .Count: N,
        .Min: Min,
        .Max: Max,
        .Mean: M,
        .Var: Var,
        .StdDev: sqrt(Var)
    }

#end // macro Statistics


#macro PrintStatistics(Stats)

    #local S = array[6] { "Count", "Min", "Max", "Mean", "Var", "StdDev" };
    #debug "\n"
    #for (I, 0, 5)
        #debug concat(S[I], " = ", str(Stats[S[I]], 0, -1), "\n")
    #end // for
    #debug "\n"

#end // macro PrintStatistics


#macro Histogram(Values, NoOfBins)

   #local Stats = Statistics(Values);
   #local N = Stats.Count;
   #local Min = Stats.Min;
   #local Max = Stats.Max;
   #local BinWidth = (Max - Min)/NoOfBins;
   #local Bins = array[NoOfBins][3];
   #local BinLimitLow = Min;
   #for (I, 0, NoOfBins - 1)
      #local BinLimitHigh = BinLimitLow + BinWidth;
      #local Bins[I][0] = BinLimitLow;
      #local Bins[I][1] = BinLimitHigh;
      #local Bins[I][2] = 0; // Set count for bin to zero
      #local BinLimitLow = BinLimitHigh;
   #end // for
   #for (I, 0, N - 1)
      // Find the right bin for the value
      #local BinNo = int((Values[I] - Min)/BinWidth);
      // Make sure rounding errors do not lead to a non existent bin
      #local BinNo = min(max(0, BinNo), NoOfBins - 1);
      // Increase the count for the bin
      #local Bins[BinNo][2] = Bins[BinNo][2] + 1;
   #end // for

   Bins

#end // macro Histogram

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7


> Suppose I could compare to forms in skvectors?
>...

I'm not sure what you mean here.


I'll comment further on the math related macros in the include files later.

--
Tor Olav
http://subcube.com
https://github.com/t-o-k

Post a reply to this message

William F Pokorny <ano### [at] anonymousorg> wrote:
> On 5/6/21 10:58 AM, Tor Olav Kristensen wrote:
> >...
>...
> > Anyway the way that macro calculates the angle is not the best numerically.
> > I would rather suggest something like this:
> >
> >
> > #macro AngleBetweenVectors(v1, v2)
> >
> >      #local v1n = vnormalize(v1);
> >      #local v2n = vnormalize(v2);
> >
> >      (2*atan2(vlength(v1n - v2n), vlength(v1n + v2n)))
> >
> > #end // macro AngleBetweenVectors
> >
>
> Interesting, thanks. I'll play with it. There will be four sqrt calls
> with this vs two so I expect it's a slower approach.

I prefer better accuracy before speed. If we need this to go faster we could
compile it into POV-Ray and make it a new command; e.g. vangle().


> Reminds me I captured fast approximation atan2 code I was thinking of
> making an inbuilt function, but I've not gotten to it...
>
> And that thought reminds me of my plans to turn many macros into
> "munctions"(1) calling new faster, inbuilt functions. Done only a tiny
> bit of that so far.
>
> (1) - Yes, thinking about creating a munctions.inc file breaking out
> what are really macro wrappers to functions into that file to make clear
> what they are. Maybe all as M_ prefixed so my current F_dents() would
> become M_dents(). Angle perhaps then M_Angle.

This sounds interesting. Can you please explain more ?
Are there any benefits from using these instead of new built in commands ?


> Something I've done too with the inbuilt function changes is allow
> alternative calculations for the same thing. Compiled such conditionals
> not that expensive. Most frequently single float options where those are
> much faster, but different methods too. For one thing, the alternatives
> I can bang against each other as a way to test!

Are you using this functionality to perform speed tests ?


--
Tor Olav
http://subcube.com
https://github.com/t-o-k

Post a reply to this message

On 5/8/21 7:06 PM, Tor Olav Kristensen wrote:
> William F Pokorny <ano### [at] anonymousorg> wrote:
>> On 5/6/21 10:58 AM, Tor Olav Kristensen wrote:
>>> ...
>>> My opinion is that several of the other math related macros in the include files
>>> need some rework.
>>>
>>
>> I'd consider alternatives for my branch if you want to offer details.
> 
> That sounds good.
> 
> Today I looked at the statistical macros in math.inc.
> 
> Here's how they could be rewritten:
> (Please note that have not tested these thoroughly yet.)
> 

Thank you. I'll have a look.

Aside: I've been working on adding self testing to all the includes 
(Ingo's suggestion and method) and in addition to macros I broke myself 
with povr branch changes, Histogram, was a macro long shipped which 
already wasn't working... ;-)

> 
> 
>> Suppose I could compare to forms in skvectors?
>> ...
> 
> I'm not sure what you mean here.

I was guessing what you have in your package would be what you consider 
solid/better implementations numerically. However, it helps me if you do 
the sorting of what might be better implementations of current code.

> 
> 

For reference I'm shipping only the following includes as 'core' ones in 
the povr branch.

arrays.inc
functions.inc
math.inc
rand.inc
shapes.inc
strings.inc
transforms.inc

Limit your look to these files for my participation. Further, I have on 
my list to move some of the macros in these to munctions.inc or other 
include files as 'not core' related / not commonly used.

arraycoupleddf3s.inc - This temporarily needed as it goes with some DF3 
related work I did years back which is in the povr branch. Some of it 
should end up in arrays.inc; some of it should go away. I've not sorted it.

Aside: My belief is future POV-Ray ray needs to move to a support 
structure which breaks things into core-code, scenes, includes, 
documentation, ini/configs, helper scripts / programs, (fonts?) and the 
web site (also the wiki?) as independent code control structures on 
github say ideally all with different primarys though perhaps still 
common POV-Ray owners like Chris. I think it would work better having 
multiple gate keepers to multiple POV-Ray domains. This model also 
matches to a degree what linux package developers are already going 
themselves breaking POV-Ray into multiple packages.

Bill P.

Post a reply to this message

On 5/9/21 2:22 AM, Tor Olav Kristensen wrote:
> William F Pokorny <ano### [at] anonymousorg> wrote:
>> On 5/6/21 10:58 AM, Tor Olav Kristensen wrote:
>>> ...
>> ...
>>> Anyway the way that macro calculates the angle is not the best numerically.
>>> I would rather suggest something like this:
>>>
>>>
>>> #macro AngleBetweenVectors(v1, v2)
>>>
>>>       #local v1n = vnormalize(v1);
>>>       #local v2n = vnormalize(v2);
>>>
>>>       (2*atan2(vlength(v1n - v2n), vlength(v1n + v2n)))
>>>
>>> #end // macro AngleBetweenVectors
>>>
>>
>> Interesting, thanks. I'll play with it. There will be four sqrt calls
>> with this vs two so I expect it's a slower approach.
> 
> I prefer better accuracy before speed. If we need this to go faster we could
> compile it into POV-Ray and make it a new command; e.g. vangle().
> 
> 
>> Reminds me I captured fast approximation atan2 code I was thinking of
>> making an inbuilt function, but I've not gotten to it...
>>
>> And that thought reminds me of my plans to turn many macros into
>> "munctions"(1) calling new faster, inbuilt functions. Done only a tiny
>> bit of that so far.
>>
>> (1) - Yes, thinking about creating a munctions.inc file breaking out
>> what are really macro wrappers to functions into that file to make clear
>> what they are. Maybe all as M_ prefixed so my current F_dents() would
>> become M_dents(). Angle perhaps then M_Angle.
> 
> This sounds interesting. Can you please explain more ?

Unsure how much you've been following my povr branch posts...

For starters, over the last year or two I've replaced most of the 
existing 'odd shape' inbuilt functions for which functions.inc defines 
the interface and documents the usage with more generic ones. I've also 
extended and cleaned up the function<->pattern interface as related 
work. New function related wave modifiers for example.

On the idea on munctions.
-------------------------
One old macro always a "munction" is the Supertorus in shapes.inc but in 
the povr branch it now calls a new inbuilt function f_supertorus() for 
backward compatibility. New users should call f_supertorus() directly.

---
We've long had simple "munctions" that look like inbuilt capability such as:

#declare clip = function (x,y,z) {min(z,max(x,y))}

which in povr I changed to

#declare F_clip = function (x,y,z) {min(z,max(x,y))}

and I'm leaning now toward calling it M_clip because I want to make 
clear it is today a munction and not an inbuilt 'clip' keyword.

Longer term F_clip/M_Clip is educational in being simple, and it will be 
left as is on the creation of a proper 'clip' implementation to create 
inbuilt include test cases comparing M_clip results to 'clip' ones.

> Are there any benefits from using these instead of new built in commands ?

Sometimes. What munctions can today do not at all doable with inbuilts 
of the straight keyword kind, or the inbuilt functions, is pre-condition 
a larger set of stuff. Stuff being data / partial solutions - or the 
particular configuration and inter operation of a collection of 
functions(1).

When you say inbuilt command, there are flavors. For performance inbuilt 
is 'almost' - but not always better.

I've gotten somewhat comfortable adding inbuilt functions, patterns and 
shapes and the associated keywords.

When you talk a something like 'acos' there is an implementation in the 
parser and there is a parallel implementation in the vm to create/run 
the acos compiled code. Not all things in the parser are in the vm and 
visa versa. Things only in the parser cannot as a rule be used with the 
vm. Where the vm things mostly (all?) work during parsing.

Though 'acos' like implementations might often be best, I don't 
understand the code well enough to pull these sorts of dual keywords at 
the moment. In any case, I'd probably go after only a few things in this 
way such as: the 'munctions' even, odd, max3, min3, sind,... (Adding the 
constant 'e' is another though that I think I can work through today - 
just not done).

I can implement inbuilt functions helping simplify/speed or replace many 
effective 'munctions.' I've been working largely in this direction of late.

---
(1) Aside: I want to get to the point where we can tell whether a 
function is being run from the parser (single thread) or during 
rendering so we can precondition functions too. I have some ideas/notes 
around thread ids and such, but I've not worked out how to really do it.

Even today in f_ridged_mf() there is an allocation of memory which I 
believe should be done only by pre-conditioning during parsing(1a). 
Function equivalents of ripples, waves, wood, etc need a parse time 
preconditioning / set up capability too, to best work.

(1a) Depending on how threaded rendering work hits f_ridged_mf() today, 
I suspect we are sometimes getting minor memory leaks. In this case to 
with respect to functionality / result I think it harmless.

> 
> 
>> Something I've done too with the inbuilt function changes is allow
>> alternative calculations for the same thing. Compiled such conditionals
>> not that expensive. Most frequently single float options where those are
>> much faster, but different methods too. For one thing, the alternatives
>> I can bang against each other as a way to test!
> 
> Are you using this functionality to perform speed tests ?
> 

Sometimes, yes.

Aside: I've had several false starts looking at integrating the work you 
and I beleive more Alain Ducharme did years ago with quaternions.inc - 
as inbuilt capability of some kind, but never gotten very far with the 
effort. Always so much on that todo list - and I'm not that smart, so it 
goes slowly. :-)

Bill P.

Post a reply to this message

On 5/9/21 9:31 AM, William F Pokorny wrote:
> Are there any benefits from using these instead of new built in commands ?

Dang..., one thing I was going to mention related to the 'munctions' 
concept and forgot.

I'm about 95% sure truly macro wrapped functions don't get 'compiled' by 
the vm until the macro is called. So, my thinking is there is a parse 
time reason to create 'macro wrapped' functions over pure function 
definitions of more complex math in our includes...

Not proven as in I haven't yet added throws to the vm code and tried 
functions defined both ways.

Bill P.

Post a reply to this message