Strange black dots showed up in some of my Mars renderings made with POVRay when the light comes in a very flat angle. These dots could not be suppressed by changing the accuracy or gradient settings (see sample scene below). So, I decided to have a look at POVRay's isosurface code and see if I can patch it to make the problem go away.

It took me quite a while to get through the code (unfortunately it completely lacks descriptive comments). I first verified that the numerical differentiation when calculating the surface normal was not the source of the problem. Then, I saw that the isosurface root finder could be improved by adding a linear interpolation after the last bisection step. Actually, this is the most common way to go numerically, and I wondered why this was not yet done already. So, this was the hour of birth for the less than 20 lines of code below.

This sample scene is a low-res version made while adjusting the light for one of the Mars renderings.

isosurface {
    function { ... }
    accuracy 1e-6
    max_gradient 3.1
    contained_by { sphere { 0, planet_radius+max_topo_height } }
}

Note that changing the accuracy value will not significantly reduce the amount of black dots in the original-code-renderings (neither when using smaller nor using larger accuracy; when using smaller accuracy, however, the amount of black spots will increase).

Original code

Patched code

Note that we cannot compare the trace time because of all these black dots. Time comparisons must be done with scenes which render well in both versions: This is the same scene as above but the light is now coming from the top.

Original code Time: 47.77 sec

Patched code Time: 47.89 sec

Additional overhead can be neglected since the calculation is not done within the innermost loop. For more timings, see below, too.

In isosurf.cpp, function IsoSurface_Function_Find_Root_R(), replace:

if(EP2->f < 0)
{
    ISOSRF->tl = EP2->t;
    return true;
}
else
    return false;

by the following code:

if(EP2->f <= 0.0)
{
    // Experimental accuracy improvement using linear
    // interpolation. -- Wolfgang Wieser 02/2004
    if(EP1->f >= 0.0)
    {
        double df = EP1->f-EP2->f;
        // Need to calc (EP1->t*EP2->f - EP2->t*EP1->f) / df
        // in a numerically stable way.
        if(df>1e-14)
            ISOSRF->tl = EP2->t + t21*EP2->f/df;
        else
            ISOSRF->tl = 0.5*(EP2->t+EP1->t);
    }
    else
        ISOSRF->tl = EP2->t;
    return true;
}
return false;

Due to the linear interpolation after the last bisection step, accuracy is improved considerably. As an example, see the following scene which is a test scene for "Water on Mars" with the water being a blue sphere and Mars topography being an isosurface just like the one above.

Original code	Patched code
Accuracy: 1e-1 Time: 4.95 sec	Accuracy: 1e-1 Time: 5.69 sec
Accuracy: 1e-2 Time: 7.53 sec	Accuracy: 1e-2 Time: 8.09 sec
Accuracy: 1e-3 Time: 11.42 sec	Accuracy: 1e-3 Time: 11.48 sec
Accuracy: 1e-4 Time: 15.10 sec	Accuracy: 1e-4 Time: 15.07 sec

In the above renderings, the gradient was set to 1 (which is still above the max gradient of the topography function according to a test trace with max_gradient set to 3).

One sees that, in this example, the patched version will give comparable results at accuracy values which are about 2 orders of magnitude larger than those needed for the original version. This saves about half the time!

This patch will not make the solver find additional roots which were not found in the original version. It will "just" give a much better guess for the actual position (i.e. ray depth) where the surface is located. This means that the patch cannot make those black regions disappear which occur because of missed roots due to e.g. too small max_gradient.

While the original bisection solver will track the root down to within a depth range of specified accuracy, the linear interpolation done by the patch will improve this guess. For very smooth functions and especially for topography functions which are linearly interpolated themselves (such as the the one used in the Mars renderings), the improvement can be one or two orders of magnitude. Expect smaller gain for not so smooth functions. Of course, there is no gain for functions which have a characteristic surface grain size which is about as large as (or even smaller than) the accuracy. This case, however, has little numerical relevance since there is no point in looking for the isosurface of a function which has surface details which are smaller than the calculated accuracy -- but even if one does that, the result of the patched algorithm is not worse than the original version.