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'