On 03/30/2018 08:21 AM, William F Pokorny wrote:
>.... 

Sturm / polysove() investigation. Chapter ??. Corruption of the Sturm Chain.

Let me start this post with the news I found something amiss in the 
recent bound estimate code recently added to:

https://github.com/wfpokorny/povray/tree/fix/polynomialsolverAccuracy

despite having run 200 plus scenes and image compares without a difference.

Caught the issue while playing with the idea of trading off the one root 
vs all roots method automatically. Not trusting the implementation; 
added an auto-bound check for lost roots and an assert and BOOM. So, I 
have more to do there even if at worst it's leaving a sanity check in 
place that resorts to an upper bound of MAX_DISTANCE. Back to that later.

---

Today let us elsewhere ponder. The function modp() does the bulk of the 
sturm chain set up. At the end of modp is this bit of code:

k = v->ord - 1;

while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH)
{
     r->coef[k] = 0.0;

     k--;
}

r->ord = (k < 0) ? 0 : k;

It individually prunes equations in the sturm chain at a depth of >= 2 
where the leading coefficients are small.

In Mr. Enzmann's original 1.0 code, 'SMALL_ENOUGH' was 1e-20. In the 
days of single floats likely meant lots of looping and looking while 
doing little. As of POV-Ray v3.0 'SMALL_ENOUGH' was drastically 
increased to 1e-10.

Over the past week or more I've been looking again at achieving better 
accuracy. Tried quite a lot of things finally coding up a __float128 
version of polysolve all aiming at more accuracy. Finally realized the 
remaining issues are far from all about accuracy. Sometimes the 
ray-surface equation coming in is simply messed up with respect to any 
sign change based methods for isolating the roots - unless we cheat.

For best accuracy while avoiding spurious roots with ill-formed 
equations, the value used should be right at the best accuracy for the 
floating point type. This would be 1e-7 for floats, 1e-15 for doubles 
and 1e-33 for 128 bit floats. With well behaved equations with no zero, 
near zero or cross term zero terms during the creation of the sturm 
chain equations can support much smaller minimum values, but in a 
raytracer we cannot count on the good behavior of the ray-surface equations.

So, why is that 'SMALL_ENOUGH' value sitting at a much larger 1e-10 
value for doubles. My guess is someone understood or stumbled upon a 
wonderful bit of serendipity in how the pruning of the sturm chain 
equations effectively work. Namely, when we have ill-formed equations 
with respect to sign chain based root isolation, the pruning value can 
be adjusted upward and by 'magic' some or all of the 'ill' in the 
incoming equation is pruned and the root isolation works well enough to 
hand the regula-falsi method a problem it can solve.

The downside is you cannot universally set this value above the minimum 
for the float size when you have well behaved incoming equations or you 
prune off sign change information you need to accurately see the 
intervals in which roots exist. We are using a fixed value today.

To better demonstrate attaching two images. In both the left column is 
what the current master renders. The remaining columns to the right are 
rendered with a float128 version of polysolve to better drive home the 
point it our primary need isn't better overall accuracy, but more 
accurate root isolation.

In the latheQuadratic.png image the top row is the point set from last 
falls lathe issue used in a quadratic spline. The top row is a 
perspective camera but zoomed way out and using a camera angle of 2. The 
bottom row is the same input except with an orthographic camera set 
up(1). In the middle column is what the current patched branch renders 
using the SMALL_ENOUGH 1e-10 value. The value has been set so as to 
behave well for the orthographic case, but this causes issues for the 
perspective case. In the right column SMALL_ENOUGH has been set to 
float128's minimum value of 1e-33. This orthographic case falls apart 
showing the underlying ill-formed for sign change isolation problem - 
while the perspective camera case is OK.

In the sphere_sweeps.png image showing the same progression(2) in 
columns 1-3 but in column 4 we dial in a higher SMALL_ENOUGH value which 
prunes off just the right about of the sturm chain for the root 
isolation to work well.

I'm busy with real life for a while starting later today, but I think we 
likely have a path to cleaner renders in these fringe cases given our 
sturm chain based method.

We're working in a ray tracer and it is the ray's dx,dy,dz going to zero 
which primarily causes the ill-formed for sign base root isolation 
issue. We'll need to detect problem directions in the object's code and 
create some calculation(3) which we'd pass all the way into the modp() 
pruning. Possible I think.

Bill P.

(1) - It's not just the orthographic camera which tends to pull out the 
worst case solver issues. Shadow rays with cylindrical and parallel 
lights too. There are too the occasional rays which line up orthogonally 
in one or compound dimensions.

(2) - Except running with jitter off. Our AA implementations are grid 
based which unfortunately tends to line up with the compound (diagonal 
type) missing root issues. Thinking something like a Halton sequence 
based AA never keeping the center pixel might really help in these cases 
by getting the samples more often off the problem ray / surface 
grid-axis alignments.

(3) - Perhaps coupled with some ability for users to apply a multiplier 
from the SDL?

Post a reply to this message

Attachments:
Download 'lathequadratic.png' (17 KB) Download 'sphere_sweeps.png' (63 KB)

Preview of image 'lathequadratic.png'

Preview of image 'sphere_sweeps.png'