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

Post a reply to this 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.

Post a reply to this message

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.

Post a reply to this message

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.

Post a reply to this 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

Post a reply to this message

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.

Post a reply to this 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

Post a reply to this 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

Post a reply to this message

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}

Post a reply to this message

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.
>

Post a reply to this message