POV-Ray : Newsgroups : povray.general : Roots of a function Server Time: 21 Sep 2020 15:22:19 GMT
 Roots of a function (Message 1 to 10 of 15)
 From: IGM Subject: Roots of a function Date: 28 Sep 2019 08:45:00 Message:
```Hi,
I have a monotone function y=(x) defined as spline (it describes the position of
an object in the scene), and I need to find its root, i.e. where it crosses x
axis. Is there any elegant way to do this inside povray? Maybe I could transform
the function in a prism surface, and using "trace" to find the intersection of a
ray at height y=0 with the surface. Other ideas?

Thanks!
igmar
```
 From: Bald Eagle Subject: Re: Roots of a function Date: 28 Sep 2019 12:10:01 Message:
```"IGM" <iar### [at] gmailcom> wrote:
> Hi,
> I have a monotone function y=(x) defined as spline (it describes the position of
> an object in the scene), and I need to find its root, i.e. where it crosses x
> axis. Is there any elegant way to do this inside povray? Maybe I could transform
> the function in a prism surface, and using "trace" to find the intersection of a
> ray at height y=0 with the surface. Other ideas?
>
> Thanks!
> igmar

You could do it the typical algebraic way and solve for x, where y=0.
You could do it numerically, where you use a #for loop and choose an increment
for x, evaluating the function and checking each iteration to see if y<= some
small amount.
If you actually want to harness [one or more of] POV-Ray's internal root
solvers, that's currently William Pokorny's domain.

If you want to do it "graphically", maybe you could do something like define a
cylinder for the x-axis, and a sphere-sweep for the function, and do an
intersection{}.  Then maybe use trace () and/or eval_pigment to find those
regions, replace with spheres encompassing those regions, redo the intersections
and find the center of the bounding box.   Iterate to a small value.

There are also online equation root solvers that you can just plug in your
equation, and it will spit the answer out from the black box.  :D

That's all I've got [so far] pre-1st-coffee.
```
 From: William F Pokorny Subject: Re: Roots of a function Date: 28 Sep 2019 13:26:40 Message: <5d8f5f90\$1@news.povray.org>
```On 9/28/19 8:08 AM, Bald Eagle wrote:
>
> "IGM" <iar### [at] gmailcom> wrote:
>> Hi,
>> I have a monotone function y=(x) defined as spline (it describes the position of
>> an object in the scene), and I need to find its root, i.e. where it crosses x
>> axis. Is there any elegant way to do this inside povray? Maybe I could transform
>> the function in a prism surface, and using "trace" to find the intersection of a
>> ray at height y=0 with the surface. Other ideas?
>>
>> Thanks!
>> igmar
>
> You could do it the typical algebraic way and solve for x, where y=0.
> You could do it numerically, where you use a #for loop and choose an increment
> for x, evaluating the function and checking each iteration to see if y<= some
> small amount.
> If you actually want to harness [one or more of] POV-Ray's internal root
> solvers, that's currently William Pokorny's domain.
>
> If you want to do it "graphically", maybe you could do something like define a
> cylinder for the x-axis, and a sphere-sweep for the function, and do an
> intersection{}.  Then maybe use trace () and/or eval_pigment to find those
> regions, replace with spheres encompassing those regions, redo the intersections
> and find the center of the bounding box.   Iterate to a small value.
>
>
> There are also online equation root solvers that you can just plug in your
> equation, and it will spit the answer out from the black box.  :D
>
> That's all I've got [so far] pre-1st-coffee.
>
>
Good suggestions though the first hard I think if not a linear spline.
Nothing simpler than what's been suggested comes to mind for the
question posed. I should first go grab my second cup...

POV-Ray's internal common solvers to which Bill refers are not directly
accessible in the scene language except I think in one of the hgpovray
releases. If remembering correctly, I'm not sure the feature is today in
the most current hgpovray38. If it is there, you'd have to be careful
using these even if your function was set up as a uni-variate power
polynomial and it's not - you have some spline. Further, the internal
solvers are set up for ray tracing, not general root solving (and there
are bugs too). I'd not recommend using them as general solvers even with
the right form of polynomials.

To your ideas I'd add using the function in an isosurface with the
default threshold of 0. Something as simple as FnIso being z-Fn(x)
should create shape edges where the spline based function crosses zero -
which you'd then need to find with trace or some other method. This can
be complicated too if your spline isn't 'well behaved' with respect to
gradients and continuity. Some sort of pigment based solution might be
more robust.

A site and web search might be worthwhile. I remember Tor Olav
Kristensen creating a simplex solver (for min, max) some while ago in
the SDL. Expect somebody, at some time past, has gone after a problem
similar to yours. Quickly finding previous work of a type can be a
difficult trick though.

Bill P.
```
 From: William F Pokorny Subject: Re: Roots of a function Date: 28 Sep 2019 15:11:32 Message: <5d8f7824\$1@news.povray.org>
```On 9/28/19 9:26 AM, William F Pokorny wrote:
> On 9/28/19 8:08 AM, Bald Eagle wrote:
>>
>> "IGM" <iar### [at] gmailcom> wrote:
...
>
> A site and web search might be worthwhile. I remember Tor Olav
> Kristensen creating a simplex solver (for min, max) some while ago in
> the SDL. Expect somebody, at some time past, has gone after a problem
> similar to yours. Quickly finding previous work of a type can be a
> difficult trick though.
>

Had the thought for what you want, given you have one root/zero-crossing
and can probable specify known values of x for lo and hi where y is
negative and positive (or visa versa) a bisection root finder probably OK.

Here is the core of my tcl bisection root finding command:

for {set i 0} {\$i < \$interationLimit} {incr i} {
set mid  [expr {(\$lo + \$hi)/2.0}]
set midV [UPolyVal \$mid coefList]
if { \$midV == 0.0 } {
return [list \$mid]
} else {
if { [expr {(\$loV * \$midV)>=0.0}] } {
set lo \$mid
} else {
set hi \$mid
}
if { [expr {(\$hi - \$lo) < \$limitVal}] } {
return [list [expr {(\$lo + \$hi)/2.0}]]
}
}
}

return [list];

Replace [UPolyVal \$mid coefList] with your Fspline(x) evaluation while
otherwise translating the tcl into SDL and, I think, you'll be good to
go. You might want to add a check the intersection limit was not reached
and kick out and error in that case instead of just returning no roots
as I do above.

Bill P.
```
 From: IGM Subject: Re: Roots of a function Date: 29 Sep 2019 07:10:01 Message:
```This is my poor-man solution:

// root(s) of function F = Y(X)-thr
// Example:
// #declare X = array [10] {0,10,20,30,40,50,60,70,80,90};
// #declare C = array [10]
{0.000,21.000,44.000,69.000,97.000,124.000,158.000,181.000,206.000,232.000};
// root(X,C,35)
// #debug concat("root: ", str(macro_value,5,3),"\n")
#macro root(X,Y,thr)

#include "math.inc"
GetStats(Y);

#local f_prism = prism {
linear_spline
-1, 1, dimension_size(Y,1)+3,

#for (k,0,dimension_size(Y,1)-1)
<X[k],Y[k]>,
#end
<X[dimension_size(Y,1)-1],StatisticsArray[2]-1>, // add lower point 1
<X[0],StatisticsArray[2]-1>, // add lower point 2
<X[0],Y[0]> // close prism
// rotate x*-90 // For plot in x,z plane. In this case P0 = <0,thr,0>
}

#object {f_prism pigment {rgb <0,2,0>} } // no_image?

#local P0 = <0,0,thr>; // Starting point
#local V0 = <1,0,0>; // Incident vector
// Trace to object
#local N = <99,99,99>; // init
#local P1 = trace(f_prism, P0, V0, N); // Loop here to find multiple roots
// Plot traced ray:
// #cylinder {P0, P1, 1 pigment {rgb <0,1,0>} no_shadow no_reflection}

// For multiple roots create array!
#declare macro_value = P1.x;
#end
```
 From: William F Pokorny Subject: Re: Roots of a function Date: 29 Sep 2019 12:48:53 Message: <5d90a835\$1@news.povray.org>
```On 9/29/19 3:07 AM, IGM wrote:
> This is my poor-man solution:
>
> // root(s) of function F = Y(X)-thr
> // Example:
> // #declare X = array [10] {0,10,20,30,40,50,60,70,80,90};
> // #declare C = array [10]
> {0.000,21.000,44.000,69.000,97.000,124.000,158.000,181.000,206.000,232.000};
> // root(X,C,35)
> // #debug concat("root: ", str(macro_value,5,3),"\n")
> #macro root(X,Y,thr)
>
>      #include "math.inc"
>      GetStats(Y);
>
>      #local f_prism = prism {
>          linear_spline
>          -1, 1, dimension_size(Y,1)+3,
>
>          #for (k,0,dimension_size(Y,1)-1)
>              <X[k],Y[k]>,
>          #end
>          <X[dimension_size(Y,1)-1],StatisticsArray[2]-1>, // add lower point 1
>          <X[0],StatisticsArray[2]-1>, // add lower point 2
>          <X[0],Y[0]> // close prism
>          // rotate x*-90 // For plot in x,z plane. In this case P0 = <0,thr,0>
>          }
>
>      #object {f_prism pigment {rgb <0,2,0>} } // no_image?
>
>      #local P0 = <0,0,thr>; // Starting point
>      #local V0 = <1,0,0>; // Incident vector
>      // Trace to object
>      #local N = <99,99,99>; // init
>      #local P1 = trace(f_prism, P0, V0, N); // Loop here to find multiple roots
>      // Plot traced ray:
>      // #cylinder {P0, P1, 1 pigment {rgb <0,1,0>} no_shadow no_reflection}
>
>      // For multiple roots create array!
>      #declare macro_value = P1.x;
> #end
>
>

Looks workable, though a few things caught my eye.

1)
If thr is 0 the Starting point will be on the prism edge and I'm not
sure what the trace result might be in that case though expect one would
want 0.0 as the returned root.

2)
IIRC you have to look at the normal returned by the trace command to
really know whether it found a valid intersection. It's set on valid
intersections only I believe.

3) And a lesson for me!
Saw #object (which you shouldn't need for the trace) and #cylinder;
Thought, "Does that really work?" It does! You can also stick '#' on the
end on some lines or code #declare CylinderZ = #cylinder {} and nary a
peep or problem from the parser - new or old versions. The usual and
recommended syntax is object {} and cylinder {}.

A declare without the leading # works too but we get:

File 'tmp.pov' line 40: Parse Warning: Should have '#' before 'declare'.
File 'tmp.pov' line 40: Possible Parse Error: 'declare' should be changed to
'#declare'. Future versions may not support 'declare' and may require
'#declare'.

Given this warning, I'd think we should be getting one on #object {} and
the like too as the other side of the syntax change push - now that I
understand there is another side.

Suppose this behavior related to Christoph's recommendation to me back
in February to use default {} rather than #default {} though it's the
latter which is still in the documentation. The idea, as I understand
it, is to get to where the #<token> is used only with SDL language
control elements.

Twenty years playing with POV-Ray and still learning.

I'm also not the brightest light bulb. I code in bash and tcl which use
'#' style comments. I've long been aware

# camera { Camera00 }

doesn't comment the camera and isolated #s are ignored because sometimes
I comment with '#' by muscle memory. I'd personally like a parse error
for floating '#'s. It would save me time and I can see no reason to have
them though '# declare Widget12 = ...' works fine today (it's treated as
#declare).

Bill P.
```
 From: IGM Subject: Re: Roots of a function Date: 29 Sep 2019 15:15:01 Message:
```> 1)
> If thr is 0 the Starting point will be on the prism edge and I'm not
> sure what the trace result might be in that case though expect one would
> want 0.0 as the returned root.

This was actually happening in my script, and the result was that the intercept
was at the other side of the prism, i.e. the last point! It deserve a dedicated
conditional check.

Thanks
igmar
```
 From: Tor Olav Kristensen Subject: Re: Roots of a function Date: 29 Sep 2019 22:25:01 Message:
```"IGM" <iar### [at] gmailcom> wrote:
> This is my poor-man solution:
>
> // root(s) of function F = Y(X)-thr
>...

If your function is strictly monotone you can just define another linear spline
with the values switched at each index, like this:

#version 3.7;

#declare NoOfPoints = 10;

#declare X =
array [NoOfPoints] {
0,
10,
20,
30,
40,
50,
60,
70,
80,
90
}
;
#declare Y =
array [NoOfPoints] {
0.0,
21.0,
44.0,
69.0,
97.0,
124.0,
158.0,
181.0,
206.0,
232.0
}
;
#declare Spline =
spline {
linear_spline
#for (I, 0, NoOfPoints-1)
X[I], Y[I]*y
#end // for
}

#declare ReverseSpline =
spline {
linear_spline
#for (I, 0, NoOfPoints-1)
Y[I], X[I]*x
#end // for
}

#declare Thr = 35;

#declare Root = ReverseSpline(Thr).x;

#debug "\n\n"
#debug concat("Root: ", str(Root, 0, -1))
#debug "\n\n"

#error "To stop text output and prevent render window"

--
Tor Olav
http://subcube.com
https://github.com/t-o-k
```
 From: Alain Martel Subject: Re: Roots of a function Date: 30 Sep 2019 22:58:13 Message: <5d928885@news.povray.org>
```Le 2019-09-29 Ã  03:07, IGM a Ã©critÂ :
> This is my poor-man solution:
>
> // root(s) of function F = Y(X)-thr
> // Example:
> // #declare X = array [10] {0,10,20,30,40,50,60,70,80,90};
> // #declare C = array [10]
> {0.000,21.000,44.000,69.000,97.000,124.000,158.000,181.000,206.000,232.000};
> // root(X,C,35)
> // #debug concat("root: ", str(macro_value,5,3),"\n")
> #macro root(X,Y,thr)
>
>      #include "math.inc"
>      GetStats(Y);
>
>      #local f_prism = prism {
>          linear_spline
>          -1, 1, dimension_size(Y,1)+3,
>
>          #for (k,0,dimension_size(Y,1)-1)
>              <X[k],Y[k]>,
>          #end
>          <X[dimension_size(Y,1)-1],StatisticsArray[2]-1>, // add lower point 1
>          <X[0],StatisticsArray[2]-1>, // add lower point 2
>          <X[0],Y[0]> // close prism
>          // rotate x*-90 // For plot in x,z plane. In this case P0 = <0,thr,0>
>          }
>
>      #object {f_prism pigment {rgb <0,2,0>} } // no_image?
>
>      #local P0 = <0,0,thr>; // Starting point
>      #local V0 = <1,0,0>; // Incident vector
>      // Trace to object
>      #local N = <99,99,99>; // init
>      #local P1 = trace(f_prism, P0, V0, N); // Loop here to find multiple roots
>      // Plot traced ray:
>      // #cylinder {P0, P1, 1 pigment {rgb <0,1,0>} no_shadow no_reflection}
>
>      // For multiple roots create array!
>      #declare macro_value = P1.x;
> #end
>
>

Why that Â«#Â» in front of Â«objectÂ» and that other one in front of
Â«cylinderÂ» ?

The Â«#Â» is needed in front of directives, never in front of a primitive.

It should be :
object {f_prism pigment {rgb <0,2,0>} }

and
cylinder {P0, P1, 1 pigment {rgb <0,1,0>} no_shadow no_reflection}
```
 From: Alain Martel Subject: Re: Roots of a function Date: 30 Sep 2019 23:04:52 Message: <5d928a14\$1@news.povray.org>
```Le 2019-09-29 Ã  08:48, William F Pokorny a Ã©critÂ :

>
> Looks workable, though a few things caught my eye.
>
> 1)
> If thr is 0 the Starting point will be on the prism edge and I'm not
> sure what the trace result might be in that case though expect one would
> want 0.0 as the returned root.
>
> 2)
> IIRC you have to look at the normal returned by the trace command to
> really know whether it found a valid intersection. It's set on valid
> intersections only I believe.

Trace return a nul vector if there is no intersection, as well as a
location of <0,0,0>.

>
> 3) And a lesson for me!
> Saw #object (which you shouldn't need for the trace) and #cylinder;
> Thought, "Does that really work?" It does! You can also stick '#' on the
> end on some lines or code #declare CylinderZ = #cylinder {} and nary a
> peep or problem from the parser - new or old versions. The usual and
> recommended syntax is object {} and cylinder {}.
>
> A declare without the leading # works too but we get:
>
> File 'tmp.pov' line 40: Parse Warning: Should have '#' before 'declare'.
> File 'tmp.pov' line 40: Possible Parse Error: 'declare' should be
> changed to
>  Â '#declare'. Future versions may not support 'declare' and may require
>  Â '#declare'.

Using declare instead of #declare can result in an error at parse time.

>
> Given this warning, I'd think we should be getting one on #object {} and
> the like too as the other side of the syntax change push - now that I
> understand there is another side.

#object should generate a warning.
Same for #Primitive_Name

>
> Suppose this behavior related to Christoph's recommendation to me back
> in February to use default {} rather than #default {} though it's the
> latter which is still in the documentation. The idea, as I understand
> it, is to get to where the #<token> is used only with SDL language
> control elements.
>
> Twenty years playing with POV-Ray and still learning.
>
> I'm also not the brightest light bulb. I code in bash and tcl which use
> '#' style comments. I've long been aware
>
> # camera { Camera00 }
>
> doesn't comment the camera and isolated #s are ignored because sometimes
> I comment with '#' by muscle memory. I'd personally like a parse error
> for floating '#'s. It would save me time and I can see no reason to have
> them though '# declare Widget12 = ...' works fine today (it's treated as
> #declare).
>
> Bill P.
>
```