END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

If our found point is vp+alpha*vcentre has distance estimate decentre from inside where vp is the viewpoint, alpha is the distance stepped and vcentre is the direction vector then I get the 4 adjacent distance estimates at vp+alpha*vleft, vp+alpha*vright, vp+alpha*vtop, vp+alpha*vbottom which gives me 5 points at vp+(alpha+decentre)*vcentre, vp+(alpha+deleft)*vleft etc. and then use the 4 triangles from the 5 points to compute the normals by summing the 4 unnormalised normals and normalising the result. Note that I do do some limited extra ray stepping (backwards) if the adjacent points are found to be truly "inside" in order to get a valid DE value obviously adjusting the calculations accordingly. I do all 4 adjacent points to allow restriction of the extra ray-stepping which means that occaissionally one or more adjacent point values are invalid. The 4 adjacent rays I use by default are at 1/2 pixel offsets.
If not using UF, or if using a global buffer in UF then it would perhaps be slightly more optimum to use whole pixel offsets for the adjacent rays.
This method is vastly superior to using the actual surface found method

My first step is to calculate the array that forms the fractal surface. The normals are then calculated from the x, y and z vectors to the immediately neighboring points.

using extra rays is basically incorrect (and outrageously inefficient besides); the gradient of the potential function is the (un-normalised) normal, or at least the finite difference approximation to it.

i use a finite difference approach as mentioned, and it's the worst-behaved bit of numerical code i've ever written

to approximate the normal at a point p, you get the potential at p, then three potentials at a "small" (herein lies the problem) offset along each dimension. then your grad vector is <px - p0, py - p0, pz - p0> where pxyz are the potentials at the offset locations. this is a first order finite difference approximation, pick up any numerical methods book to find better appoximations (bearing in mind that they assume the existence of higher derivatives, i.e. greater smoothness).

if your potential function is p(x,y,z),

p0 = p(x,y,z)
px = p(x + d, y, z)
py = p(x, y + d, z)
pz = p(x, y, z + d)

then your normal vector is normalise(px - p0, py - p0, pz - p0)

the trouble with the delta d is that you want it to be as small as possible, but if you make it too small then you destroy all the floating point precision in the estimate. see for example http://en.wikipedia.org/wiki/Numerical_differentiation#Practical_considerations

this is exacerbated by the fact that our function p isn't some nice x^2 + 2x or whatever, it's this nasty fractal function  

OK, but what is the potential function?  How do we get from the XYZ point in space to the potential?

These posts have stirred a lot of interest so hopefully there will be someone reading with the smarts to work out how to automatically find the optimal d value?

OK Lycium, I think I know what you mean. See the sketch below for the 2D analogy. Is my assumption correct ?

So in 3D you need the analytical distance estimate in three extra points. I will try this..

What I have been doing is shooting 4 extra rays. If the point on the surface is VP+t.V (VP=viewpoint, V=vector), I take VP+(t-epsilon).Vi (Vi =vector with slight offset, epsilon=something small). So I go back just a small distance on the ray so that finding the intersection points for the vectors Vi goes very fast. (but I have to calculate the normalised vectors Vi of course)