In continuing to look at media I returned to a 'blob-media' issue Gail 
Shaw originally posted back in 2005. The bright spots at the very edge 
are indeed caused by the 'intersection depth < small_tolerance' or 
missed intersection issue as Slime guessed at the time.

However, the sturm result Gail showed never made any sense unless sturm 
is broken for blobs. Gail's original belief something changed 3.5 to 3.6 
was correct. Optimizations were done involving bounding and polynomial 
reformulation for 3.6 which occasionally inhibit 4th to 3rd order 
reduction(1). The bright spots not on the edge happen where 
Solve_Polynomial passes the original 4th order blob equation to 
polysolve and we get bad roots back.

(1) - Happens mostly with transparency.

In the lower left in the attached image I hacked the code to prevent 
order reduction with blobs. Surprised to find it looks like blobs only 
work as formulated due order reduction - even in v3.5.

All got me digging again into the solver code. Since v1.0 there has been 
a regula_falsa function called once the rule of signs stuff has 
determined a sub-interval contains a single root. The regula_falsa is 
buggy once polynomial order is >= 4 given how many of our polynomials 
are formulated - the evaluated values can be very small.

I've coded up a branch modifying the regula_falsa function so it sanity 
checks the root found is a real root - within some tolerance - otherwise 
it invalidates the root and tosses control back to the interval 
bisection code.

Branch fixes this blob issue, the lathe issue #335 from last fall and 
many remaining #147 (sphere_sweep) solver issues - except one b_spline / 
orthogonal camera case. The patch also fixes or makes better many of the 
accuracy issues I have in hand, however, I have too still some scenes 
with 'accuracy' issues to dig through. The work is tedious and I need a 
break - so putting out the fix I have.

A branch of master with the code update is available at:

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

for those compiling their own POV-Ray versions. I plan to create a pull 
request once I've looked at some of the remaining accuracy issues.

Aside: The regula_falsa function isn't strictly needed for root finding, 
but the solver runs faster with it than not.

Bill P.

Post a reply to this message

Attachments:
Download 'blobsturmbroken.jpg' (129 KB)

Preview of image 'blobsturmbroken.jpg'

On 03/30/2018 08:21 AM, William F Pokorny wrote:
> 
> Branch fixes this blob issue, the lathe issue #335 from last fall and 
> many remaining #147 (sphere_sweep) solver issues - except one b_spline / 
> orthogonal camera case. The patch also fixes or makes better many of the 
> accuracy issues I have in hand, however, I have too still some scenes 
> with 'accuracy' issues to dig through. The work is tedious and I need a 
> break - so putting out the fix I have.
> 
> A branch of master with the code update is available at:
> 
> https://github.com/wfpokorny/povray/tree/fix/polynomialsolverAccuracy
> 
> for those compiling their own POV-Ray versions. I plan to create a pull 
> request once I've looked at some of the remaining accuracy issues.
> 

I've run down another type of solver issues related to our sor object. 
In this case the problem is with the sor polynomial set up often not 
crossing through 0.0 which is a requirement for the sturm / polysolve 
methods.

The ray / sor set up sometimes creates sets of coefficients which are 
are all positive or negative which sturm / polysolve sees as meaning 
there are no roots to find. The fix is to perturb the direction vector 
by +EPSILON during the set up which is in my testing enough to get the 
crossing 0.0 criteria always.

I started with a sor scene shown in the top row of the attached image 
which Severi Salminen posted in 2004. v381 shown on the left and the fix 
on the right.

(The remaining three rows show differences at 4x the real magnitude.)

The second row is our shipped sor1 sample scene which doesn't use sturm 
at all. I expected no difference, but the 'sturm off', solve_cubic() 
function is also failing to find roots on reflected rays which I didn't 
previously notice. The +EPSILON change perturbs these noisy solution 
regions.

The third row is the same sor1 scene, but with sturm turned on for all 
sor objects meaning polysolve gets used. The difference on the right of 
row 3 shows how a noisy original solution is cleaned up with the fix.

The last row shows the updated code's solve_cubic solution vs the 'sturm 
on' / polysolve solution.  Note that 'sturm on' is needed for clean 
roots & a clean sor1.pov image.

Aside: Yes, it's looking like continued rays whether transmitted or 
reflected tend to be harder for the solvers to handle. This echoes back 
to early lemon issues with noise when the texture had a reflective 
component. Also true errors in continued rays are often harder to see 
given the contributions are mixed with any ambient or diffuse components.

I've more solver accuracy issues to run down, but I've updated the 
branch at:

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

with this sor polynomial set up fix for those interested.

Bill P.

Post a reply to this message

Attachments:
Download 'sorsoverstory.jpg' (317 KB)

Preview of image 'sorsoverstory.jpg'

On 03/30/2018 08:21 AM, William F Pokorny wrote:
> ...
> 
> A branch of master with the code update is available at:
> 
> https://github.com/wfpokorny/povray/tree/fix/polynomialsolverAccuracy
> 
> for those compiling their own POV-Ray versions. I plan to create a pull 
> request once I've looked at some of the remaining accuracy issues.
> ...

I've continued to work on the polysolve 'sturm' solver in recent weeks 
and another set of updates is available at:

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

A couple quick notes followed at the bottom by more detailed change and 
performance information. With respect especially to sphere_sweeps these 
changes fix a number of bad or slightly bad root issues, but do not fix 
missing root issues as shown in the attached image.

Another recent revelation has to do with the previous observation 
continued rays seem to be more difficult for the sturm solver. There are 
two major reasons for this.

One. The regula falsi method, while usually performing better than 
bisection and guaranteeing a solution on proper set up, can - with 
certain equations - be very, very slow to converge. Our transmitted, 
internally reflected and self shadowing test rays all tend to create 
polynomials which are a struggle for the regula_falsa() code.

Two. The polysolve code was always set up to jump to regula_falsa() once 
it was known only one root existed on a 'ray' interval. The ray types 
mentioned previously start with only one root as do rays with respect to 
some surfaces. In these situations we do the work to set up the sturm 
chain, use it to verify there is just one root - and immediately hand 
often difficult to solve equations to regula_falsa() with the original 
0.0 to MAX_DISTANCE interval. In these cases we are really using the 
regula falsi method.

The changes below re-work things enough so as to be able to eliminate 
the previously added sanity check at the end of regula_falsa(). Other 
updates in the works.

Bill P.


Performance info:
------------ /usr/bin/time povray -j -wt2 -fn -p -d lemonSturm.pov

0)  master    30.29user 0.02system 0:15.83elapsed
0)  3f2e1a9   30.31user 0.00system 0:15.81elapsed

1)  7be07ff   29.99user 0.01system 0:15.65elapsed  -1.02%
2)  7f7eddd   29.91user 0.01system 0:15.61elapsed
3)  a621c04   19.63user 0.01system 0:10.45elapsed  -34.46%
4)  b2042e6   19.67user 0.02system 0:10.45elapsed
5)  58a1933   20.36user 0.01system 0:10.80elapsed   +3.51%
5a) 1c13a17   20.24user 0.00system 0:10.75elapsed
6)  7258366   19.33user 0.02system 0:10.31elapsed   -5.06%
7)  f313a5b   19.19user 0.02system 0:10.25elapsed   -0.72%
8)  cbd3af5   19.34user 0.01system 0:10.30elapsed   +0.78%
9)  687f57b   19.16user 0.01system 0:10.21elapsed   -0.93%   -36.77%


1) Fix for earlier regula_falsa sturm chain based sanity check where the 
sturm chain was shortened during set up and the sanity check was using 
the full length.

2) Removed a normalization of the base polynomial added in the late 90s.

3) Delaying regula_falsa use until ray domain range < 1.0.

Important particularly where only 1 root. Previously the sturm chain 
used only to verify the single root after which regula_falsa was called 
with an initial range of 0 to MAX_DISTANCE.

4) Updating regula_falsa to remove two initial SMALL_ENOUGH fast paths.

These filters return quickly where the regula_falsa method performs 
poorly - very small values at one end of the range and much larger at 
the other. However, the cost is the roots often being very inaccurate 
where the polynomial values are small over more than one end of the range.

5) Updating regula_falsa to use ray domain for the RELERROR conditions.

Previously a mix but primarily the polynomial value domain. Led to 
different root values at about 5-6 digits in the ray domain depending 
upon whether sturm chain bisection came up with the root or regula_falsa 
did. With update results identical out to 8-10 digits.

5a) Adding missing else condition in previous commit.

6) Removing regula_falsa root sanity check.

With changes in recent commits the sanity check is no longer needed.

7) Increasing the MAX_ITERATIONS limits in sbisect and regula_falsa.

Even with recent updates found in practice sbisect sometimes needs 
nearly 60 iterations and regula_falsa just over 100. Previously used 
identical limits of 50. Now 65 and 130 respectively.

8) Update sbisect to return no roots when ray domain range is tiny.

Roots being very closely packed as happens most often when a ray is 
nearly tangent to a surface cause instability in the sturm change sign 
counting (numchanges) results. The instability, if not avoided, leads to 
incomplete and sometimes inaccurately ordered roots. One of the triggers 
for the CEY patch added in 1997 and where we get the diagonal lines of 
no roots in certain objects.

9) Update sbisect to return a single root with closely spaced roots.

This is the current practice in the non-sturm chain based solvers and 
expected by at least some primitive objects. Root returned is the middle 
of the tiny range at which we stop bisection.

Post a reply to this message

Attachments:
Download 'lipka_gh147_fs81.png' (50 KB)

Preview of image 'lipka_gh147_fs81.png'

William F Pokorny <ano### [at] anonymousorg> wrote:

[Quite a lot of stuff]

First of all - great work!
I can appreciate how much you need to pick and poke and test and plow through in
order to unravel the source code, design test objects, and tease out the little
problematic details and then conceive of and implement solutions.

This is definitely where Stephen shall christen this work as the plural MATHS.

I only have the most rudimentary and tenuous grasp of what you're working on,
having never been taught these methods [that I can recall], and only getting a
sample of how they're used from the many and varied excursions into math videos
on the web.
I went to an American public indoctrination center.   You should see the
abomination they call "common core".  It's worse than the New Math.

https://thefederalistpapers.org/education-2/kid-looks-like-a-genius-shows-exactly-why-common-core-math-doesnt-add-up

https://www.youtube.com/watch?v=UIKGV2cTgqA

Just for clarification, POV-Ray uses:
bisection
regula falsa
Newton's method, and
the Sturmian method
?
Or are some of those the same?

and - where is this all located in the source code?
Because I'm a hopeless intellectual masochist like that.

Post a reply to this message

On 04/26/2018 10:16 AM, Bald Eagle wrote:
> William F Pokorny <ano### [at] anonymousorg> wrote:
> 
...
> 
> Just for clarification, POV-Ray uses:
> bisection
> regula falsa
> Newton's method, and
> the Sturmian method
> ?
> Or are some of those the same?
> 
> and - where is this all located in the source code?
> Because I'm a hopeless intellectual masochist like that.
> 

Thanks. Suppose I've been exposed to some of the concepts in the past, 
but trust me, I'm climbing the learning curve/mountain(a) with this 
stuff. There is much I don't understand. My standard set of test scenes 
number a 195 at present which have corrected many of my missteps.

You can find most of the code in source/core/math/polynomialsolver.cpp.

There are three custom solvers in solve_quadratic(), solve_cubic() and 
solve_quartic() for orders 2 through 4. I've hardly looked at these.

In polysolve / 'sturm' - a sturm chain based bisection method is used in 
combination with regula_falsa(). Value based bisection isn't used at 
all. Newton–Raphson isn't used at all - though the first two equations 
in the sturm chain are, in my present version, the original and it's 
normalized derivative. Who knows, might play with it as an alternative 
or adjunct to regula_falsa(). The spin is it's really fast when it works 
- but sometimes it doesn't.

Near term looking more at the missing root issue - some(b) of which was 
addressed with (9), but the rest comes to the roots having collapsed 
toward each other so severely the sturm chain root counting / isolation 
falls apart. Not looking like it will be easy to mitigate...

Bill P.

(a) - Trying to get out of the media whirlpool! :-)

(b) - Where the sturm chain root counting and isolation stuff still 
works despite how close the roots are.

Post a reply to this message

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

While studying the missing root problem I finished off some initial 
lower and upper root bound estimation code aimed at polysolve() 
performance. The update is at:

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

The update is directly related to a sub-discussion of issue:

https://github.com/POV-Ray/povray/issues/236

Currently the upper bound method(d) is hard coded to run always. Longer 
term, and given I've found some POV-Ray objects such as blob already 
calculate a 'max_bounds' value, I'm thinking about a change to the 
Solve_Polynomial(), polysolve() calls to add a parameter where :

a) if > 0 value, we'd be passing a known upper bound for all roots.

b) if == 0,  we'd do what we do today.

c) if == -1, we'd run the original Cauchy's bound estimate method good 
for the lowest value root. In some applications, that first root might 
be all we need.

d) if == -2, we'd run an estimator based on Cauchy's method from a 2009 
paper(z) found via:

https://en.wikipedia.org/wiki/Properties_of_polynomial_roots#Other_bounds

and good for the upper bound of all roots.

There are better bound estimators, but all I've looked at are more 
complex - and several I don't understand. In any case I doubt we'd buy 
back the performance given our equation orders are tiny compared to that 
aimed at by the more generic solvers for which these better estimator 
were created. Many others also use complex numbers.

Bill P.

(z) - Had to tweak it slightly due two sor fails in testing. Due the 
surprise with sturm-cubic it's on my list to enable sturm-quadratic in 
Solve_Poly() to do some more testing.

Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -c lemonSturm.pov

0)  master    30.22user 0.04system 0:30.89elapsed

9)  687f57b   19.44user 0.01system 0:20.06elapsed  -35.67%
10) 0f6afcc   14.77user 0.02system 0:15.36elapsed  -24.02%  -51.13%

Post a reply to this message

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'

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

Sturm / polysove() investigation. Those Difficult Coefficients.

When we have 4th order equations we are today running this piece of code:

   /* Test for difficult coeffs. */

   if (difficult_coeffs(4, c))
   {
       sturm = true;
   }

in the Solve_Polynomial wrapper for sturm/polysolve() and three specific 
solvers including the one for 4th order equations in solve_quartic(). 
The difficult_coeffs() function scans the coefficients looking for large 
differences in values. If found, user end up running sturm for some rays 
no matter what was specified in the scene SDL. The code existed in the 
original POV-Ray 1.0, but in a slightly modified form and used with an 
older version of solve_quartic().

I've had the feeling for a while the code was probably doing nothing 
useful for us today. I was almost completely right. Close enough it's 
going into the bit bucket with this set of commits.

The code was set up to crash anytime difficult_coeffs() returned true. 
Scanned hundreds of scenes. Turned up helix.pov and a handful of scenes 
with blobs. Ran those hard coded to use only solve_quartic() and not a 
single difference in result. Changed the code above to return 
immediately with no roots on difficult_coeffs() returning true. The 
helix scene was unchanged as the 'rays' were not part of the resultant 
shape.

The blobs were a different story. Many rays are being re-directed to 
sturm / polysolve for no benefit given the current solve_quartic code. 
To oblivion with that code you all say!

I thought the same, but remembered we have known blob accuracy issues in 
https://github.com/POV-Ray/povray/issues/187, maybe half a dozen related 
problem scenes in my test space. Plus to support the subsurface feature 
Chrisoph added this bit of code in blob.cpp to avoid visible seams:

   DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);

I've been running with DEPTH_TOLERANCE at 1e-4 for a long time in my 
personal version of POV-Ray which fixes all the blob accuracy issues 
except for the subsurface one.

Given difficult_coeffs() is almost exclusively active with blobs, 
updated blob.cpp to run with accuracy addressing the issues in the 
previous paragraph. With this update in place, we do occasionally see 
cases where the bump into to sturm / polysolve method is necessary. 
However, the automatic difficult_coeffs() method doesn't currently catch 
all of the fringe case rays.

It could be adjusted so as to push more solve_quartic rays into sturm / 
polysolve, but it's aleady the case most difficult-coefficient rays are 
unnecessarily run in sturm / polysolve when solve_quartic would work 
just fine. The right thing to do is let users decide when to use 'sturm' 
and this set of commit dumps the difficult_coeffs() code.

Examples.

The attached image is a scene originally created to explore default 
texturing related to negative blobs, but it demonstrates nicely the 
issues described above.

In the top row we have on the left the image rendered with the blob 
accuracy unmodified and all the difficult_coeffs() code running as it 
does today. In the middle image the 'if difficult_coeffs()' test has 
been hard coded to return no roots. On the right the difference image 
shows all the rays running in the sturm / polysolve method which all run 
fine - and faster given the blob set up - with solve_quartic().

In the second row we are running with the new more accurate blob.cpp. 
Further, the blobs in the scene have been changed to have a threshold of 
0.001 instead of 0.1 - values this low are not common because you don't 
get much blobbing.

The first column of the second row was rendered with a version of 
POV-Ray compiled to always force the use of solve_quartic(). The middle 
column restores the 'if difficult_coeffs()' test and internal swap to 
sturm / polysolve. It does help as shown in column 3. However column 2 
shows the already too aggressive use of sturm / polysolve with blobs is 
not enough to automatically 'fix' all the solve_quartic issues. Running 
with sturm is clean in all cases.

Aside: The potentials as glows image I did sometime back had what I 
thought were media related speckles. I applied very aggressive and 
expensive AA to hide them as I usually do with media. Turns out using 
sturm would have been the better choice. I wanted no blobbing for the 
blob cylinders and used a threshold of 1e-5. It was a case of outrunning 
the effectiveness of the 'if difficult_coeffs() patch' with existing 
blob accuracy, but I didn't realize it at the time.

Updates at:

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

Performance and specific commit details below.

Bill P.

Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -c lemonSturm.pov

0)  master    30.22user 0.04system 0:30.89elapsed

10) 0f6afcc   14.77user 0.02system 0:15.36elapsed  -51.13%
11) 0f9509b   14.77user 0.01system 0:15.36elapsed
12) c590c0e   15.08user 0.01system 0:15.71elapsed   +2.10%
13) 4ab17df   14.88user 0.02system 0:15.51elapsed   -1.33%  -50.76%
14) 78e664c   --- NA ---
15) 5d8cedc   --- NA ---
16) 21387e2   --- NA ---

---- povray subsurface.pov -p +wt2 -j -d -c -fn
master  sturm on       540.89user 0.06system 4:31.15elapsed
21387e2 sturm on       347.29user 0.07system 2:54.42elapsed  -35.79%

master  sturm off      286.99user 0.06system 2:24.28elapsed
21387e2 sturm off      278.54user 0.05system 2:20.12elapsed   -2.94%

If blob container needs sturm.  278.54 -> 347.29  +24.68%


11) Changing polysolve visible_roots and numchanges zero test to <=0.

Depending upon value used for reducing and trimming the sturm chain 
equations in modp, the sphere_sweep b_spline functions can return 
negative root counts when there are no roots.

12) Adding sanity test of upper bound Cauchy method result.

In practice have occasionally found roots slightly above upper bound. 
Perhaps related to how sturm chain is tweaked during pruning. Sort it later.

13) Removing pointless operations in polysolve and Solve_Polynomial.

visible_roots() passed a couple integer pointers which it set prior to 
return but the values were never used in polysolve itself.

When fourth order equations seen a test call to difficult_coeffs() was 
made which if true would set sturm true - even when sturm was already true.

14) Updating blob accuracy so better synchronized with other shapes.

Fix for GitHub issue #187

Returned minimum intersection depths now set to MIN_ISECT_DEPTH.

Internal determine_influences() intersect_<sub-element> functions using 
 >0.0 as was already the subsurface feature. Anything more opens up gaps 
or jumps where the sub-element density influences are added too late or 
dropped too earlier. The change does create ray equations slighly more 
difficult to solve which was likely the reason for the original, 
largish, v1.0 1e-2 value. The sturm / polysolve solver is sometimes 
necessary over solve_quartic().

15) Removing difficult_coeffs() and related code in Solve_Polynomial().

Code not completely effective where sturm / polysolve needed and very 
costly in most scenes with 4th order equations due use of sturm over 
solve_quartic() where solver_quartic works fine and is much faster.

16) Reversing change to sor.cpp as unneeded with updated polysolve().

Not verified, but suspect bump made in commit 3f2e1a9 to get sign change 
was previously needed where regula_falsa() handed single root problem 
with 0 to MAX_DISTANCE range. In any case, the patch is no longer needed 
with updated polysolve() code.

Post a reply to this message

Attachments:
Download 'storydifficultcoeff.png' (110 KB)

Preview of image 'storydifficultcoeff.png'

On 05/18/2018 01:00 PM, William F Pokorny wrote:
> On 03/30/2018 08:21 AM, William F Pokorny wrote:
> ....
> 
> Sturm / polysove() investigation. Those Difficult Coefficients.
> 

I'll add quickly this set of updates fixes the remaining media speckles 
for the scene which started this thread - see attached image. The 
general media speckling issue is better with more accurate roots, but 
not completely fixed.

Bill P.

Post a reply to this message

Attachments:
Download 'gailshawexampleclean.png' (167 KB)

Preview of image 'gailshawexampleclean.png'

William F Pokorny <ano### [at] anonymousorg> wrote:
> On 03/30/2018 08:21 AM, William F Pokorny wrote:
>
[snip]
> I've had the feeling for a while the code was probably doing nothing
> useful for us today. I was almost completely right. Close enough it's
> going into the bit bucket with this set of commits.
> ...
>
> The blobs were a different story. Many rays are being re-directed to
> sturm / polysolve for no benefit given the current solve_quartic code.
> To oblivion with that code you all say!
>
[snip]

Your detailed research into these arcane issues is really appreciated. I can't
say I understand some of it (most of it!), but we are all lucky to have your
keen eye probing the code.

Post a reply to this message

Am 18.05.2018 um 19:00 schrieb William F Pokorny:

> Chrisoph added this bit of code in blob.cpp to avoid visible seams:
> 
>   DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);
> 
> I've been running with DEPTH_TOLERANCE at 1e-4 for a long time in my
> personal version of POV-Ray which fixes all the blob accuracy issues
> except for the subsurface one.

For subsurface rays it's safer to keep it at exactly 0.

The `depthTolerance` mechanism, to all of my knowledge, is in there to
make sure that transmitted, reflected or shadow rays don't pick up the
original surface again.

This problem is absent in rays earmarked as subsurface rays (presuming
we have high enough precision) because such rays are never shot directly
from the original intersection point, but a point translated slightly
into the object. So subsurface rays already carry an implicit
depthTolerance with them.

> It could be adjusted so as to push more solve_quartic rays into sturm /
> polysolve, but it's aleady the case most difficult-coefficient rays are
> unnecessarily run in sturm / polysolve when solve_quartic would work
> just fine. The right thing to do is let users decide when to use 'sturm'
> and this set of commit dumps the difficult_coeffs() code.

Would it make sense to, rather than ditching the test entirely, instead
have it trigger a warning after the render (akin to the isosurfaces
max_gradient info messages) informing the user that /if/ they are seeing
artifacts with the object they may want to turn on sturm (provided of
course sturm is off)?

Post a reply to this message

On 05/18/2018 10:16 PM, clipka wrote:
> Am 18.05.2018 um 19:00 schrieb William F Pokorny:
> 
>> Chrisoph added this bit of code in blob.cpp to avoid visible seams:
>>
>>    DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);
>>
>> I've been running with DEPTH_TOLERANCE at 1e-4 for a long time in my
>> personal version of POV-Ray which fixes all the blob accuracy issues
>> except for the subsurface one.

Formally a blob's 'DEPTH_TOLERANCE' value was used for both the returned 
from ray origin intersection tolerance and the internal 
determine_influences() 'mind_dist' tolerances - though they control 
different things. The internal sub-element influences mechanism now uses 
 >0.0. Should have always I believe except it does sometimes create 
slightly more difficult ray-surface equations - which is my bet for why 
the original v1.0 code didn't.

> 
> For subsurface rays it's safer to keep it at exactly 0.
> 
> The `depthTolerance` mechanism, to all of my knowledge, is in there to
> make sure that transmitted, reflected or shadow rays don't pick up the
> original surface again.
> 

It depends about which intersection depth tolerance you speak. There are 
variations of the mechanism.

a) What you say is close to my understanding for the continued and 
shadow ray mechanisms which exists outside the code for each shape. 
These mechanisms live in the 'ray tracing' code. I'd slightly reword 
what you wrote to say these mechanisms prevent the use of ray - surface 
intersections within some distance of a previous intersection. This 
distance is today usually SMALL_TOLERANCE or 1e-3. SHADOW_TOLERANCE has 
also been 1e-3 for more than 25 years and was only every something 
different because Enzmann's initial solvers had a hard coded shadow ray 
filter at a larger the SMALL_TOLERANCE value. There is too 
MIN_ISECT_DEPTH.

b) Root / intersection depth tolerance filters exist in the code for 
each shape/surface too. These control only whether all roots / 
intersections for the shape are greater than a minimum intersection 
depth relative to the incoming ray's origin or less than MAX_DISTANCE. 
They do not control how closely ray surface intersections are spaced 
from each other(1)(2). True the shape depth tolerance can influence 
continued or shadow ray handling as seen by being 'different' than those 
used in (a). However, we should aim for (b) <= (a) depth tolerance 
alignment.

On seeing your subsurface code running the (b) variety at >0.0. I had 
the thought, why not do this with all shapes / surfaces? Let (a) filter 
the roots for its need. This would get us closer to a place where we'd 
be able to make SMALL_TOLERANCE & MIN_ISECT_DEPTH smaller so people 
don't have to so often do the ugly scale up their scene trick. There is 
no numerical reason we cannot support generally much smaller scene 
dimensions than we do.

Well it turns out returning all the roots >0.0 in each shapes (b) 
mechanism costs cpu time - about 1% of a 3% increase with the internal 
blob influence min_dist change for a glass blob scene. Time is saved 
where intersections which (a) would have to filter later are just never 
returned by (b). This experiment is why I set blob.cpp's DEPTH_TOLERANCE 
= MIN_ISECT_DEPTH - which only happens also be the 1e-4 I've long been 
using. I currently think the (b) depth for each shape / surface should 
be migrated to MIN_ISECT_DEPTH over time.

(1) - I expect CSG generally needs all the roots regardless of spacing, 
but this is code I've never reviewed.

(2) - Mostly true for what I've seen. But who knows, might be other 
inbuilt shape, surface implementations do spacial root filtering for 
internal reasons. The sphere_sweep's code seems to be doing a sort of a 
root correction thing - with which you are far more familiar than me.


> This problem is absent in rays earmarked as subsurface rays (presuming
> we have high enough precision) because such rays are never shot directly
> from the original intersection point, but a point translated slightly
> into the object. So subsurface rays already carry an implicit
> depthTolerance with them.
> 

I tested subsurface scenes with my updates and believe the changes are 
compatible without the need for the returned depth to be >0.0. Thinking 
about it another way, if the 'returned' min depth needed to be >0.0 for 
subsurface to work, all shapes would need that conditional for >0.0 and 
this is not the case. Believe we're good, but we can fix things if not.

>> It could be adjusted so as to push more solve_quartic rays into sturm /
>> polysolve, but it's aleady the case most difficult-coefficient rays are
>> unnecessarily run in sturm / polysolve when solve_quartic would work
>> just fine. The right thing to do is let users decide when to use 'sturm'
>> and this set of commit dumps the difficult_coeffs() code.
> 
> Would it make sense to, rather than ditching the test entirely, instead
> have it trigger a warning after the render (akin to the isosurfaces
> max_gradient info messages) informing the user that /if/ they are seeing
> artifacts with the object they may want to turn on sturm (provided of
> course sturm is off)?
> 

We should think about it, but currently I answer no to your specific 
question. The difficult coefficient method is wildly inaccurate in 
practice unlike the warning with isosurfaces.

A warning would be good if we can figure out something better... Mostly 
it will be obvious - and more often obvious with this update. If running 
blobs and you see speckles / noise - try sturm. Advice good today by the 
way. Advice I should have taken with my own potential glow scene over 
extreme AA - but the code dropped hid much of the solve_quartic issue 
from me and I thought media the cause.

I ran test scenes where knowing you need sturm is pretty much impossible 
to see without trying sturm and comparing images. For example, attached 
is an image where I was aiming for a fruit looking thing with a faked 
subsurface effect using multiple concentric blobs positive and negative 
- taking advantage of blob's texture blending.

In the left column we see the <= master (3.8) current result where the 
internal blob surfaces are not accurately resolved. The right column 
shows the difference between column 1 and 2 for each row.

The top row, middle shows the new all solve_quartic() solution which in 
fact has bad roots on a surface internal to the outer ones. Negative 
blobs in addition to tiny thresholds tend to trip the solve_quartic 
failures. We cannot 'see' we've got bad roots... The bottom row is the 
accurate - though it happens less attractive - sturm result.

Reasonable accurate warnings would be good. Wonder if there is a way to 
pick the issue up in solve_quartic itself or in blob.cpp knowing sturm 
not set....

Bill P.

Post a reply to this message

Attachments:
Download 'diffcoeffhidden.png' (160 KB)

Preview of image 'diffcoeffhidden.png'

On 05/18/2018 10:16 PM, clipka wrote:
> Am 18.05.2018 um 19:00 schrieb William F Pokorny:
> 
>> Chrisoph added this bit of code in blob.cpp to avoid visible seams:
>>
>>    DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);
>>
>> I've been running with DEPTH_TOLERANCE at 1e-4 for a long time in my
>> personal version of POV-Ray which fixes all the blob accuracy issues
>> except for the subsurface one.
> 
> For subsurface rays it's safer to keep it at exactly 0.
> 
...
> 
> This problem is absent in rays earmarked as subsurface rays (presuming
> we have high enough precision) because such rays are never shot directly
> from the original intersection point, but a point translated slightly
> into the object. So subsurface rays already carry an implicit
> depthTolerance with them.
> 

Aside: 'exactly 0' unless we change blob.cpp code means >0.0.

I woke up this morning thinking more about your description of SSLT and 
though I picked up no major differences in result, results were not 
identical - but I didn't expect them to be. I'd not run the returned 
root depth threshold at both 0.0 and 1e-4 (MIN_ISECT_DEPTH) and compared 
results.

The top row in the attached image is without sturm the bottom with. The 
differences are very subtle, but can be seen after being multiplied by 50.

The left column is less my recent changes so min_dist and the return 
depth >0.0 with sturm automatically getting used often in the upper left 
image.

The middle column is the updated code so the middle top now using 
strictly solve_quartic() making the differences on the right are 
greater. The bottom, sturm always row best shows the effect of not 
returning roots <= 1e-4 from the ray origin.

My position is still the difference isn't large enough to be of concern, 
but if so, maybe we should more strongly consider going to >0.0 or some 
very small values for all shapes. Doing it only in blob.cpp doesn't make 
sense to me.

Aside: MIN_ISECT_DEPTH doesn't exist before 3.7. Is it's creation 
related to the SSLT implementation?

Bill P.

Post a reply to this message

Attachments:
Download 'sslt_rtrnrootquestion.png' (460 KB)

Preview of image 'sslt_rtrnrootquestion.png'

Am 20.05.2018 um 18:10 schrieb William F Pokorny:

> My position is still the difference isn't large enough to be of concern,
> but if so, maybe we should more strongly consider going to >0.0 or some
> very small values for all shapes. Doing it only in blob.cpp doesn't make
> sense to me.

I guess the reason I changed it only for blobs was because there I
identified it as a problem - and I might not have been familiar enough
with the shapes stuff to know that the other shapes have a similar
mechanism.

After all, why would one expect such a mechanism if there's the "global"
MIN_ISECT_DEPTH?

Technically, I'm pretty sure the mechanism should be disabled for SSLT
rays in all shapes.


> Aside: MIN_ISECT_DEPTH doesn't exist before 3.7. Is it's creation
> related to the SSLT implementation?

No, it was already there before I joined the team. The only relation to
SSLT is that it gets in the way there.

My best guess would be that someone tried to pull the DEPTH_TOLERANCE
mechanism out of all the shapes, and stopped with the work half-finished.


BTW, one thing that has been bugging me all along about the
DEPTH_TOLERANCE and MIN_ISECT_DEPTH mechanisms (and other near-zero
tests in POV-Ray) is that they're using absolute values, rather than
adapting to the overall scale of stuff. That should be possible, right?

Post a reply to this message

On 05/20/2018 05:09 PM, clipka wrote:
> Am 20.05.2018 um 18:10 schrieb William F Pokorny:
> 
>> My position is still the difference isn't large enough to be of concern,
>> but if so, maybe we should more strongly consider going to >0.0 or some
>> very small values for all shapes. Doing it only in blob.cpp doesn't make
>> sense to me.
> 
> I guess the reason I changed it only for blobs was because there I
> identified it as a problem - and I might not have been familiar enough
> with the shapes stuff to know that the other shapes have a similar
> mechanism.
> 
> After all, why would one expect such a mechanism if there's the "global"
> MIN_ISECT_DEPTH?

Understand. Only recently are the root filtering and bounding mechanisms 
coming into some focus for me & I'm certain I still don't understand it 
all. There are too, the two, higher level bounding mechanism which are 
tangled in the root/intersection handling. I'm starting to get the 
brain-itch the current bounding is not always optimal for 'root finding.'

> 
> Technically, I'm pretty sure the mechanism should be disabled for SSLT
> rays in all shapes.
> 

OK. Agree from what I see. Work toward >0.0 I guess - though I have this 
feeling there is likely an effective good performance 'limit' ahead of 
 >0.0.

> 
>> Aside: MIN_ISECT_DEPTH doesn't exist before 3.7. Is it's creation
>> related to the SSLT implementation?
> 
> No, it was already there before I joined the team. The only relation to
> SSLT is that it gets in the way there.
> 
> My best guess would be that someone tried to pull the DEPTH_TOLERANCE
> mechanism out of all the shapes, and stopped with the work half-finished.
> 
> 
> BTW, one thing that has been bugging me all along about the
> DEPTH_TOLERANCE and MIN_ISECT_DEPTH mechanisms (and other near-zero
> tests in POV-Ray) is that they're using absolute values, rather than
> adapting to the overall scale of stuff. That should be possible, right?
> 

Sure, but... Not presently in POV-Ray is my 'anytime soon' answer. Too 
many magic values or defined values like EPSILON used for differing -or 
effectively differing - purposes. Remember, for example, one of the 
recent changes I made to the regula-falsi method used within polysolve() 
was to use universally the ray value domain instead of a mix of 
polynomial values and ray values for relative & absolute error stops.

Root / intersection work is sometimes done in a normalized space(1) and 
sometimes not, etc. The pool of issues, questions and possibilities in 
which I'm mentally drowning is already deep - and I'm far from seeing to 
the bottom. We can often scale scenes up for better result because 
POV-Ray code has long been numerically biased to 0+ (1e-3 to 1e7) (2). 
Near term think we should work toward something zero centered (1e-7 to 
1e7).

Would like to first get to where we've got a more accurate polysolve() 
against which any changes to other solvers / specialized shape solvers 
can be tested. After which we can perhaps work to get all the solvers & 
solver variant code buried in shape code to double / 'DBL' accuracy.

Underneath everything solver / intersection wise is double accuracy 
which is only 15 decimal digits. The best double value step 
(DBL_EPSILON) off 1.0 is 2.22045e-16. Plus we use fastmath which 
degrades that accuracy.

Aside: I coded up a 128bit polysolve() but it was very slow (+10x) and 
the compiler feature doesn't look to be supported broadly enough as a 
standard to be part of POV-Ray. A near term interesting idea to me - one 
Dick Balaska, I think, touched upon recently in another thread -  is 
coding up a version of polysolve() we'd ship using 'long double'. Long 
double is guaranteed at only 64bits (double), but in practice today 
folks would mostly get 80bits over 64bits - and sometimes more. It would 
need to be an option to the default 64 bit double as it would certainly 
be slower - but being generally hardware backed, performance should be 
<< the +10x degrade I see for 128bits on my Intel processor. Would be - 
for most - an improvement to 18 decimal digits and an epsilon step of 
1.0842e-19 off 1.0. Worthwhile...?

Bill P.

(1) Normalizing & the inverse, transforms & the inverse all degrade 
accuracy too.

(2) From the days of single floats - where the scale up would have 
worked less well at hiding numerical issues (ie coincident surface 
noise) because the non-normalized / global accuracy was then something 
like 6 decimal digits or 1e-3 to + 1e3.

Post a reply to this message

Am 21.05.2018 um 16:39 schrieb William F Pokorny:

> Underneath everything solver / intersection wise is double accuracy
> which is only 15 decimal digits. The best double value step
> (DBL_EPSILON) off 1.0 is 2.22045e-16. Plus we use fastmath which
> degrades that accuracy.
> 
> Aside: I coded up a 128bit polysolve() but it was very slow (+10x) and
> the compiler feature doesn't look to be supported broadly enough as a
> standard to be part of POV-Ray. A near term interesting idea to me - one
> Dick Balaska, I think, touched upon recently in another thread -  is
> coding up a version of polysolve() we'd ship using 'long double'. Long
> double is guaranteed at only 64bits (double), but in practice today
> folks would mostly get 80bits over 64bits - and sometimes more. It would
> need to be an option to the default 64 bit double as it would certainly
> be slower - but being generally hardware backed, performance should be
> << the +10x degrade I see for 128bits on my Intel processor. Would be -
> for most - an improvement to 18 decimal digits and an epsilon step of
> 1.0842e-19 off 1.0. Worthwhile...?

I think that's a path worthwhile exploring further. My suggestion would
be to use a macro - say "PRECISION_FLOAT" - to define a type to be used
in places where precision matters even more than usual, so that we can
get binaries with different performance characteristics out of the same
source code: That macro could simply be defined as `double` for the
status quo, `long double` for a portable attempt at achieving more
precision without losing too much speed, a compiler-specific type such
as `__float128` on gcc, or even a type provided by a 3rd party library
such as GMP or MPFR.

> (2) From the days of single floats - where the scale up would have
> worked less well at hiding numerical issues (ie coincident surface
> noise) because the non-normalized / global accuracy was then something
> like 6 decimal digits or 1e-3 to + 1e3.

POV-Ray has never used single-precision floats. I doubt even DKBTrace
ever did: At least as of version 2.01, even on the Atari it already used
double-precision floats.

As a matter of fact, the IBM-PC configuration for DKBTrace 2.01 (*)
actually defined DBL as `long double` when a coprocessor was present.
DKBTrace 2.12 (**) additionally set EPSILON to 1.0e-15 instead of 1.0e-5
in that case (provided the compiler used wasn't Turbo-C, presumably
because it would define `long double` as 64-bit). This was carried on at
least until POV-Ray v2.2, but must have been gotten lost between v2.2
for DOS and v3.6 for Windows. (My guess would be that the Windows
version never had it.)

(* See here:
https://groups.google.com/forum/#!topic/comp.sources.amiga/icNQTp_txHE)

(** See here: http://cd.textfiles.com/somuch/smsharew/CPROG/)

Post a reply to this message

Am 21.05.2018 um 22:35 schrieb clipka:

> I think that's a path worthwhile exploring further. My suggestion would
> be to use a macro - say "PRECISION_FLOAT" - to define a type to be used

Make that `PRECISE_FLOAT`; fits better with a series of similar types
already in use, `PreciseColourChannel`, `PreciseRGBColour` etc., which
are used for colour-ish data where `float` is not enough, most notably
the SSLT code. (Obviously, there the data type of choice is `double`,
and does suffice.)

Post a reply to this message

On 05/21/2018 04:35 PM, clipka wrote:
> Am 21.05.2018 um 16:39 schrieb William F Pokorny:
> 
> POV-Ray has never used single-precision floats. 

Taking you to mean never in the 'solvers'.

> I doubt even DKBTrace
> ever did: At least as of version 2.01, even on the Atari it already used
> double-precision floats.
> 
> As a matter of fact, the IBM-PC configuration for DKBTrace 2.01 (*)
> actually defined DBL as `long double` when a coprocessor was present.
> DKBTrace 2.12 (**) additionally set EPSILON to 1.0e-15 instead of 1.0e-5
> in that case (provided the compiler used wasn't Turbo-C, presumably
> because it would define `long double` as 64-bit). This was carried on at
> least until POV-Ray v2.2, but must have been gotten lost between v2.2
> for DOS and v3.6 for Windows. (My guess would be that the Windows
> version never had it.)
> 
> (* See here:
> https://groups.google.com/forum/#!topic/comp.sources.amiga/icNQTp_txHE)
> 
> (** See here: http://cd.textfiles.com/somuch/smsharew/CPROG/)
> 

Thanks. I've never dug as far back. True too I was speculating about how 
many might have actually been running early POV-Ray. Me taking the early 
as v1.0 "#define DBL double to be a mechanism those without FPUs and 
early small bus sizes were running with "define DBL float" instead. Not 
that I was running POV-RAY or even aware of POV-Ray that early.

Don't remember exactly, but I didn't have a personal machine with an 
actual FPU until the late 90s. First 64 bit machine early 2000s I guess. 
My very first personally owned machine was an 8bit z80 based custom 
built thing...

Bill P.

Post a reply to this message

On 05/21/2018 05:13 PM, clipka wrote:
> Am 21.05.2018 um 22:35 schrieb clipka:
> 
>> I think that's a path worthwhile exploring further. My suggestion would
>> be to use a macro - say "PRECISION_FLOAT" - to define a type to be used
> 
> Make that `PRECISE_FLOAT`; fits better with a series of similar types
> already in use, `PreciseColourChannel`, `PreciseRGBColour` etc., which
> are used for colour-ish data where `float` is not enough, most notably
> the SSLT code. (Obviously, there the data type of choice is `double`,
> and does suffice.)
> 

OK. Will go for PRECISE_FLOAT and plan to define it in 
./base/configbase.h along with DBL and SNGL.

Further - for now - I plan to use it only in polysolve() and its 
sub-functions with conversions from and to DBL on entry and exit as I 
implemented the 128bit polysolve().

Bill P.

Post a reply to this message

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

Sturm / polysove() investigation. Coefficients wrongly seen as zeroes.

In this chapter we venture up into the Solve_Polynomial() function which 
calls polysove(), solve_quadratic(), solve_cubic() and solve_quartic(). 
Consider this code added back in 1994 and part of POV-Ray v3.0.

     /*
      * Determine the "real" order of the polynomial, i.e.
      * eliminate small leading coefficients.
      */

     i = 0;

     while ((i < n) && (fabs(c0[i]) < SMALL_ENOUGH))
     {
         i++;
     }

     n -= i;

     c = &c0[i];

My last commit deletes it! Having so large a SMALL_ENOUGH (1e-10) and 
any fixed value given equation order varies as does distance along the 
ray where the equations not normalized is a bad idea(1).

Perhaps reason enough that when any even order polynomial with no roots 
is reduced to an odd order polynomial it turns into a polynomial with 
roots. All odd order polynomials have at least one root.

In the attached image we have two examples - neither using sturm / 
polysolve - where the code above goes bad.

In the top row we've got the shipped primativ.pov scene and specifically 
a poly shape which is not explicitly bounded. It creates 4th order 
equations so should use solve_quartic(). However the code up top 
sometimes strips the leading coefficient and in this case we go from a 
quartic which has no roots to an odd third order cubic which has one. 
This is the vertical bar artifact in the upper left.

In the bottom row we have the shipped subsurface.pov scene. The 
differences on the bottom right (at 4x actual difference values) happen 
again where the code at top drops a single leading coefficient. In this 
case for pixels and rays where I've examine actual roots, subsurface 
related rays go from 4th order equations with >=2 roots with one of 
those roots inside the blob interval being kept, to a 3rd order equation 
with a single root outside the active blob interval. In other words we 
are missing subsurface roots/intersections where we see differences.

Notes:

- The subsurface effect here is visible with no difference multiplier 
and many times larger than the returned root threshold differences 
recently discussed.

- Because the ray-surface equations are based upon scene geometry the 
'zeroes-not-zeroes-dropped' issue here is often scene-object or axis 
aligned in some fashion. I'm reminded of a rounded box canyon subsurface 
scene Robert McGregor posted where I noticed one box face having a 
strange look. Someone joked it was perhaps a door to Narnia or similar 
after nobody had an idea what caused the effect. I'm thinking perhaps 
the code above - or the still present order reduction by ratio code - 
played a part...

NOTE! This update corrupts some previous scenes while fixing some 
previous bad ones. Some planned but not yet defined or implemented 
equation tuning on ray DX,DY,DZ values is to come.

Updates at:

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

Performance and specific commit details below.

Bill P.

(1) - Yes, this strongly suggests order reduction based upon the largest 
two coefficients is also a bad idea. However leaving that feature in 
place for now as some primitives require this type of reduction to 
function at all. Perhaps in these cases push the reductions out of the 
solvers and into the shape code.


Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -c lemonSturm.pov

0)  master    30.22user 0.04system 0:30.89elapsed

16) 21387e2   14.88user 0.02system 0:15.51elapsed    -50.76%
17) 80860cd   --- NA --- (1e-15 constant calc fix)
 
(PRECISE_FLOAT)
18) aa6a0a6   --- Bad Result(a) ---                           (float)
               14.95user 0.02system 0:15.56elapsed     +0.47%  (double)
               25.12user 0.02system 0:25.73elapsed    +68.82%  (long double)
               282.97user 0.13system 4:44.13elapsed +1801.68%  (__float128)
19) 4e16623   14.94user 0.01system 0:15.56elapsed      ----   -50.56%

(a) - Not exhaustively evaluated, but some shapes render OK at 'float' 
though those which did were slower (+150% ballpark). Unsure why so slow, 
but don't plan to dig further. Sturm / polysolve is for those looking to 
be more accurate, not less.


17) Correcting two double 1e-15 constant calculations.

In recent commits started setting two 1e-15 values with the calculation:

= 1/std::numeric_limits<DBL>::digits10;

when I meant to code:

= 1/pow(10,std::numeric_limits<DBL>::digits10);

Both forms bad practice in not being specific as to the types for the 
explicit values. The first returns integer 0, later cast to double and I 
got away with it in all but my full set of test cases. Corrected:

= (DBL)1.0/pow((DBL)10.0,std::numeric_limits<DBL>::digits10);

18) Implementing suggested PRECISE_FLOAT macro mechanism for polysolve().

Primarily enabling 'long double' support for the polysolve() / 'sturm' 
solver given 'long double' is part of the C+11 standard and gets most 
users 80 bits or more as opposed to 64. Further, the additional bits are 
almost always hardware backed making for reasonable performance.

With GNU g++ 64 bit environments '__float128' is a valid setting with no 
additional library required - though it's very slow. Other extended 
accuracy floats such as 'double double' enabled, via the PRECISE_ 
mechanisms, but they involve additional libraries and settings. See code.

Note! With this update the polysolve sturm chain is no longer pruned 
with a default 1e-10 value, but rather a value appropriate for the 
floating point type. Means better accuracy for some previously bad root 
results and less accuracy for others. Coming updates will better tune 
the chain equation pruning.

19) Removing leading, near zero coefficient reduction in Solve_Polynomial().

Leading coefficients of < SMALL_ENOUGH (1e-10) were being dropped 
reducing the order of the incoming polynomial when the near zero values 
had meaning with respect to roots. Led to both false roots and missed roots.

Post a reply to this message

Attachments:
Download 'notzeroesdropped.png' (326 KB)

Preview of image 'notzeroesdropped.png'

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

Sturm / polysove() investigation. Other solvers. Root polishing.

I've continue to look at the root solvers over the past month. At a good 
place to write another post about findings and changes.

A couple preliminaries before the main topic.

First, while most results are better, the polynomialsolverAccuracy 
branch still lacks planned updates. For example, scenes with lathes and 
sphere_sweeps using orthographic cameras currently render with more 
artifacts, not less.

Second, this set of commits updates the solve_quadratic function in a 
way both more accurate and faster.

---
The main topic is root finding tolerance and scaling.

All the solvers find roots / surface intersections to some tolerance. 
They do this in a coordinate space which has been normalized with 
respect all scene transforms. Additionally, objects like blobs 
internally normalized to internal 'unit' coordinate spaces.

Suppose in the solver coordinate space a ray/surface equation has a root 
essentially at the ray origin - at distance 0.0 - as happens with self 
shadowing rays. All solvers due their tolerances have the potential to 
allow such near-zero roots to drift into the >0.0 value space as 
ray-surface intersections. This numerical reality is why objects filter 
roots to a minimum returned intersection depth larger than zero(1).

Further, this minimum returned intersection depth must account for the 
translation of every root and its tolerance back into the original 
coordinate space. A wrongly positive 1e-8 root in solver space is a much 
more positive 1e-4 root after a scale up by 10000. The solvers each have 
different tolerances. There are many shape custom 'solvers' in use 
beyond those I've been working upon and differing internal 
normalizations. Why we've got many differing returned depth values. Why 
blob's have long used a minimum returned intersection depth of 1e-2 
while a sphere, for example, uses 1e-6.

It's the case that even blob's large 1e-2 filtering value is inadequate 
given solve_quartic()'s tolerance is relatively large. Continuing to use 
Gail Shaw's blob scene from 2005 to demonstrate, the attached image was 
rendered with the current 3.8 master branch without sturm - so using 
solve_quartic(). The original scene render is on the left. On the right 
is the same scene scaled up by 10000. The vertical bars showing up on 
the right happen due slightly ~=0.0 roots drifting positive in 
solve_quartic(). The wrongly positive roots are then scaled up to a 
value greater than the current 1e-2 small root/intersection depth 
filter. Not being filtered, the roots corrupt results.

While zero crossing roots due tolerance are the most serious tolerance 
multiplied by scale issue, roots well positive can drift quite a lot too 
in the global space where the scale up is large.

The improvement adopted for the solve_quartic() and polysolve()/sturm 
solvers - where 'scaled tolerance' issues have been seen - was to add a 
Newton-Raphson root polishing step to each. With blobs this looks to 
allow a returned depth on the order of 4e-8 over 1e-2 for a 1e-7 to 1e7 
global working range at DBL = 'double'.

Aside for thought: Might we be able to determine when we are evaluating 
a self-shadowing or same-shape-terminating ray-surface equation? If so, 
it should be we can use the knowledge we have a root 'at 0.0' to always 
deflate / reduce the order of these polynomials prior to finding roots.

Updates at:

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

Performance and specific commit details below.

Bill P.

(1) - The idea of going as small as >0.0 in each shape's implementation 
fly is not presently possible. Trying instead a new much smaller 
MIN_ISECT_DEPTH_RETURNED value - see commit comments.


Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -c lemonSturm.pov

0)  master    30.22user 0.04system 0:30.89elapsed

19) 4e16623   14.94user 0.01system 0:15.56elapsed   -50.56%
20) 0d66160   15.22user 0.03system 0:15.80elapsed    +1.87%
               (hilbert_curve linear sphere_sweep
                scene with new solve_quadratic())     (-7.5%)
21) 761dd0b   15.05user 0.04system 0:15.71elapsed    -1.12%
22) 552b625   (quartic polish non-sturm lemon scene) (+0.82%)
23) 63e0456   15.16user 0.02system 0:15.76elapsed    +1.20%
24) ec28851   NA Code documentation only.
25) 8962513   15.43user 0.02system 0:16.05elapsed    +1.78%  -48.94%
               (quartic polish non-sturm lemon scene) (+1.02% +2.70%)


19) Removing leading, near zero coefficient reduction in Solve_Polynomial().

Leading coefficients of < SMALL_ENOUGH (1e-10) were being dropped 
reducing the order of the incoming polynomial when the near zero values 
had meaning with respect to roots. Led to both false roots and missed roots.

20) New zero coefficient stripping and solve_quadratic implementation.

General effort to implement better 'effective zero' coefficient 
handling. Created new POV_DBL_EPSILON macro value which is 2x the C++ 
standards <float type>_EPSILON value and updated PRECISE_EPSILON to be 
2x the single bit epsilon as well.

Zero filtering in polysolve now looks at all polynomial coefficients and 
either sets 'effective zeros' to exactly 0.0 or strips them if they are 
leading coefficients.

The much more accurate near zero coefficients drove the need for a 
better solve_quadratic implementation included with this commit. Note it 
supports the PRECISE_FLOAT options like polysolve.

Zero filtering in solve_quadratic, solve_cubic and solve_quartic now 
standardized both in implementation and use of POV_DBL_EPSILON.

21) Created new constexpr DBL variable MIN_ISECT_DEPTH_RETURNED.

Near term need to use a value not MIN_ISECT_DEPTH in blob.cpp to test 
Jérôme's github pull request #358. Longer term aim is to drive all 
returned intersection depths from shape code to 
MIN_ISECT_DEPTH_RETURNED. The value is automatically derived from DBL 
setting and at double resolves to about 4.44089e-08. On the order of the 
square root of POV_DBL_EPSILON.

Moving blob.cpp's inside test value INSIDE_TOLERANCE to POV_DBL_EPSILON 
over previous recent re-calculation. Over time plan is to move EPSILONs 
best nearer a double's step to POV_DBL_EPSILON.

Cleaning up doxygen documentation added during recent solver related 
updates.

22) Adding root polishing step to solve_quartic function.

Newton-Raphson step added to polish initial roots found by the 
solve_quartic function. The core solve_quartic tolerance allows roots to 
drift from <=0.0 to >0.0 values with the latter causing artifacts and 
additional root filtering. This the reason for the long too large 1e-2 
intersection depth value in blob.cpp now reduced to 
MIN_ISECT_DEPTH_RETURNED.

As part of this change created a new FUDGE_FACTOR4(1e-8) constant DBL 
value to replace previous use of SMALL_ENOUGH(1e-10) within 
solve_quartic. Looks like the value had been smaller to make roots more 
accurate, but at the cost of missing roots in some difficult equation 
cases. With the root polishing can move back to a larger value so as to 
always get roots. Yes, this better addresses most of what the already 
remove difficult_coeffs() function and bump into polysolve was trying to do.

23) Adding root polishing step to polysolve function.

24) Cleaning up a few comments and doxygen documentation.

25) Moving to more conservative root polishing implementations.

Post a reply to this message

Attachments:
Download 'tolerancestory.png' (97 KB)

Preview of image 'tolerancestory.png'

On 3/30/18 8:21 AM, William F Pokorny wrote:
> In continuing to look at media I returned to a 'blob-media' issue Gail 
> Shaw originally posted back in 2005. 

A long time since I posted here on the progress I guess... Still hacking 
- just proving hard to get everything dependent on polysolve to 
something mostly 'fixed.'

Yesterday, Johannes (j13r) posted in newusers a scene where he was 
trying to create a scene using torus shapes and media to create 
something like solar flares on the surface of a disk. A slightly 
modified version without lights makes for a pretty good poster for 
current media issues compared to my current working branch with better 
solvers/shape code, my current version of the change proposed in: 
https://github.com/POV-Ray/povray/pull/358 and the internal small 
tolerance at 1e-6 instead of 1e-3.

The current v3.8 on the left, my patched branch in the middle and the 
differences (artifacts) shown on the right. The scene was modified for 
no lights and run with the flags: +w600 +h300 Corona.pov -j +p +a in 
both cases.

Bill P.

Post a reply to this message

Attachments:
Download 'coronaa_to_coronab.jpg' (21 KB)

Preview of image 'coronaa_to_coronab.jpg'

Am 25.11.2018 um 17:18 schrieb William F Pokorny:

> A long time since I posted here on the progress I guess... Still hacking 
> - just proving hard to get everything dependent on polysolve to 
> something mostly 'fixed.'

I admire and appreciate your perseverance in this matter.

Post a reply to this message

On 3/30/18 8:21 AM, William F Pokorny wrote:
> In continuing to look at media I returned to a 'blob-media' issue Gail 
> Shaw originally posted back in 2005. The bright spots at the very edge 
> are indeed caused by the 'intersection depth < small_tolerance' or 
> missed intersection issue as Slime guessed at the time.
> 

Sturm / polysove(), solvers investigation. End of this thread.

Almost exactly a year ago started to dig into numerical issues. It's a 
good time to take a break. Before I do, five more commits to close out 
this thread with updates representing a reasonably complete set of 
improvements. I have in mind still many possible solver improvements and 
code improvements beyond these. I could spend the rest of my lifetime at 
it. I'm instead going to get back to some scene work near term.

First commit a set of changes I've been sitting on since last fall which 
to my testing makes most everything better solver-side wise. The one 
user oriented change was adding a sturm option to sphere_sweep where 
that sturm means a two pass sort of sturm. The sturm control tests a fix 
for a class of issues most often seen with orthogonal cameras to an in 
plane sweep as might be done for human signatures or presentation like 
images.

Second commit is my version of Jerome's pull request #358 which needs 
some of the changes in the first commit. Note this update causes scenes 
with photons to be somewhat slower due more photons being deposited.

Third commit addresses what might be called 2.5 issues with the shadow 
cache for both artifacts and performance. The update causes a shift in 
certain shadow artifacts. I suspect the state of the shadow cache might 
have come about as a way to hide (it didn't completely and causes other 
issues) certain shadowing issues - but guessing. The also cleans up 
SMALL_TOLERANCE and SHADOW_TOLERANCE uses in the shadow cache / trace 
shadow code. The values have been equivalent for 25+ years so the 
tangling had not mattered to result.

Fourth commit mostly a completion of the second and third eliminating 
SMALL_TOLERANCE in addition to MIN_ISECT_DIST. This partly follows up on 
a question Jerome asked either in #358 or some related issue. Namely, 
why not zero or near it for the Intersect_BBox_Dir calls in object.cpp. 
Larger values do indeed cause artifacts with secondary rays as 
especially noticeable in scenes with media. It's one of many causes for 
media speckles. Pull #358 as originally submitted moved from a 
MIN_ISECT_DIST value of 1e-4 to SMALL_TOLERANCE of 1e-3 making the 
inherent problem worse. SMALL_TOLERANCE had also been adopted a few 
other places in the code. In those cases moved to gkDBL_epsilon or 
gkMinIsectDepthReturned as appropriate.

Fifth commit changes Cone_Tolerance(1e-9) in cone.cpp to 
gkMinIsectDepthReturned(4.4e-8). Found to be necessary during testing 
due pull #358 changes.

Solvers in much better shape to my testing (2000+ cases now) but, so 
long as we run more or less directly on the floating point hardware for 
best speed, there will always be floating point accuracy issues cropping up.

All shapes where I made modifications test better and scale better, but 
I'll mention the sphere_sweeps still has substantial issues with media 
(even at 1x scale) and scaling. Due - reasons - one being the necessary 
updates are not easy. Rays perpendicular to a sweeps directions will 
sometimes show MORE artifacts due tightening up on solver accuracy(1) - 
this is where the new sturm (double sturm) control must be used.

The two attached images shows a kind of media scaling set I've generally 
adopted as it seems to be a pretty good way to test the numerical 
soundness of a shape's code - especially the secondary rays. One image 
is blobs for sturms the other for the updated solve_quartic code. Middle 
row being the updated and essentially speckle-less new result. Seventh 
column over is the 1x scene. Scales from the left at 1e6x to the right 
at 1e-7x.

Updates at (2):

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

Performance and specific commit details below.

Bill P.

(1) - As mentioned previously some of the solvers had been de-tuned 
(solve_quadratic and sturm) to help lathes and sphere_sweeps - I guess. 
The solve_quartic solvers was tuned - perhaps by accident given the 
previous use of SMALL_ENOUGH - to be more accurate at the expense of 
finding roots. In that latter case now tuned to find the most roots 
possible with root polishing. Argh! I'd need to write a book to describe 
anything close to all the details addressed and still open. For the 
record significant solver related conversation can be found in the 
github pull request comments at: https://github.com/POV-Ray/povray/pull/358.

(2) - The last set of changes in master forced some branch merges 
instead of the usual re-basing. The solver branch, for one, had to be 
merged to maintain compile-ability on checkout of previous branch commits.

Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -cc lemonSturm.pov

0)  master    30.22user 0.04system 0:30.89elapsed

25) 8962513   15.43user 0.02system 0:16.05elapsed     +1.78% -48.94%
               (quartic polish non-sturm lemon scene) (+1.02% +2.70%)
26) 75ddd88   17.53user 0.02system 0:18.11elapsed    +13.61%
               (quartic non-sturm lemon scene)        (+3.10%)
27) 4ea0c37   NA (slowdown due photons. benchmark +3.5%)
28) ff6cd8d   14.41user 0.02system 0:15.01elapsed    (-17.80%)
29) ef9538b   NA
30) 38b434d   NA
...
25) Moving to more conservative root polishing implementations.

Further making constant names consistent with recommened coding style. In
polynomialsolver.cpp results in the following name changes.

FUDGE_FACTOR2 now kSolveQuarticV1_Factor2
FUDGE_FACTOR3 now kSolveQuarticV1_Factor3
FUDGE_FACTOR4 now kSolveQuarticV2_Factor4
TWO_M_PI_3 now kSolveCubic_2MultPiDiv3
FOUR_M_PI_3 now kSolveCubic_4MultPiDiv3
MAX_ITERATIONS now kMaxIterations
SBISECT_MULT_ROOT_THRESHOLD now kSbisectMultRootThreshold
REGULA_FALSA_THRESHOLD now kRegulaFalsaThreshold
RELERROR now kRelativeError
SMALL_ENOUGH now kSolveQuadratic_SmallEnough

26) Initial, reasonably complete, update to common solvers.

Working on possible further improvements, but those likely quite far out 
in time. In total dozens of issues addressed. New solver call structure. 
Solvers themselves more accurate and aligned better with common 
practice. Additional root polishing. Corresponding updates to shape code 
while working to extend scaling range for shapes. Changes encompass 
updates needed to support a commit to follow which covers pull request 
#358 and its associated issues.

Generally sturm option now much faster. Also true due improvements to 
the fixed solvers that the sturm option is less often necessary.

The sphere_sweep now supports a sturm option. Here sturm runs the 
sturmian solver in a two pass approach which, for certain scenes, will 
work where the previous single pass sturmian solver did not. Unlike 
other updated shapes, the sphere_sweeps scaling accuracy range was not 
much improved due a decision to leave other parse time optimizations for 
run time performance in place.

27) My implementation of Jerome's pull request #358.

Fix for issues #121, #125 and several related newsgroup reports as well. 
Mostly it restores 3.6 behavior with respect to intersection depths 
filtered.

Note! Scenes using photons will often run slower and look somewhat 
different due additional photons being deposited. This includes our 
benchmark scene.

28) Shadow cache fixes. SHADOW_TOLERANCE vs SMALL_TOLERANCE cleanup.

SHADOW_TOLERANCE not used for cached results leading to artifacts though 
sometimes hiding others.

Shadow cache not invalidated on misses causing sometimes significant 
performance hit.

SMALL_TOLERANCE being used instead of SHADOW_TOLERANCE in some trace 
shadow related comparisons. Values had long (25+ years) been cleaned up 
ahead of SMALL_TOLERANCE removal.

Note! This the last commit where SMALL_TOLERANCE will exists.

29) Mostly details completing the previous two commits.

Eliminating SMALL_TOLERANCE in addition to MIN_ISECT_DIST. This partly 
follows up on a question Jerome asked either in #358 or some related 
issue. Namely, why not zero or near it for the Intersect_BBox_Dir calls 
in object.cpp. Larger values do indeed cause artifacts with secondary 
rays as especially noticeable in scenes with media. It's one of many 
causes for media speckles. Pull #358 as originally submitted moved from 
a MIN_ISECT_DIST value of 1e-4 to SMALL_TOLERANCE of 1e-3 making the 
inherent problem worse. SMALL_TOLERANCE had also been adopted a few 
other places in the code. In those cases moved to gkDBL_epsilon or 
gkMinIsectDepthReturned as appropriate.

30) In cone.cpp changing Cone_Tolerance to gkMinIsectDepthReturned.

During testing found the Cone_Tolerance in cone.cpp, which had been 
changed from 1e-6 to 1e-9 v3.6 to v3.7, was too small for some photon 
scenes. Secondary rays starting on the surface (at zero but numerically 
not) were not getting filtered with the pull request 358 like changes 
(MIN_ISECT_DIST to 0.0). Moved to new gkMinIsectDepthReturned (4.4e-8) 
used now in many other shapes and OK for test scenes I have.

Post a reply to this message

Attachments:
Download 'blobsturmstory.png' (78 KB) Download 'blobnosturmstory.png' (74 KB)

Preview of image 'blobsturmstory.png'

Preview of image 'blobnosturmstory.png'