Improved delta-DE algorithm

[Buddhi]

Hi
Because rendering on GPU has limitation for calculation accuracy (mostly possible to use only floats), I had to redesign a little delta-DE algorithm. Standard 4-point delta-DE needs delta lower about 1000 times than distance threshold. With single rescission floating point numbers the lowest acceptable delta is anout 1e-6. Lower value produces big numeric errors. Higher gives artifacts in areas where DE is not continuous function. With this limitation there was not possible to get images with magnification higher than 1e3.
To improve delta-DE algorithm I decided to use additional 3 points. It's better to calculate deltas in positive and negative direction and use lower value of them. Then non-continuous areas of DE are eliminated.

Below there is a code which I use actually in Mandelbulber-OpenCL
 

Code:

formulaOut CalculateDistance(__constant sClInConstants *consts, float3 point, sClCalcParams *calcParam)
{
	float distance;
	float delta = 1e-6;
	float3 dr = 0.0;
	formulaOut out = Fractal(consts, point, calcParam);

	if (out.iters == calcParam->N)
	{
		distance = 0.0f;
	}
	else
	{
		float r = out.distance;
		float r11 = Fractal(consts, point + (float3)	{ delta, 0.0, 0.0}, calcParam).distance;
		float r12 = Fractal(consts, point + (float3)	{ -delta, 0.0, 0.0}, calcParam).distance;
		dr.x = min(fabs(r11 - r), fabs(r12 - r)) / delta;
		float r21 = Fractal(consts, point + (float3)	{ 0.0, delta, 0.0}, calcParam).distance;
		float r22 = Fractal(consts, point + (float3)	{ 0.0, -delta, 0.0}, calcParam).distance;
		dr.y = min(fabs(r21 - r), fabs(r22 - r)) / delta;
		float r31 = Fractal(consts, point + (float3)	{ 0.0, 0.0, delta}, calcParam).distance;
		float r32 = Fractal(consts, point + (float3)	{ 0.0, 0.0, -delta}, calcParam).distance;
		dr.z = min(fabs(r31 - r), fabs(r32 - r)) / delta;
		float d = length(dr);
	}
	out.distance = distance;
	return out;
}

where Fractal().distance returns last z value from fractal formula iterations.

[Roquen]

It would seem desirable to have the delta be a value with no significant bits (a power of two) so you can convert the divide into a multiple and remove one rounding error from the process.  Simple to implement since CL support hex constants.  The closest to the current would be: 0x1.0-20f

[eiffie]

Thanks for this improvement!

[lycium]

Using central-differences as you do turns the error in Taylor series from O(h) to O(h^2). Roquen's suggestion suggests he's pretty experienced in numerical methods, right on the money!

More info here: http://en.wikipedia.org/wiki/Finite_difference

[Roquen]

Another random though is to make the delta adaptive.  So:  norm = dot(p,p); delta = nextAfter(norm)-norm;  then multiple delta some small power of two to allow for compounding errors in the computation of the function.

[Syntopia]

Using central-differences as you do turns the error in Taylor series from O(h) to O(h^2). Roquen's suggestion suggests he's pretty experienced in numerical methods, right on the money!

More info here: http://en.wikipedia.org/wiki/Finite_difference

I don't think it is a central-difference. Buddhi calculates both a forward and a backwards difference, and chooses the minimum of these two.

Buddhi, could you make a central-differences image for comparison? I think that would be just:

Code:

dr.x = (r11 - r12) / (2.0*delta); // sign doesn't matter

[Roquen]

My concrete mathematics is weak, but isn't the error analysis similar?  BTW: My thinking about the adaptive delta is if that were to work, I'd be curious if an intermediate difference might be sufficient...and thus drastically increase performance.

[Buddhi]

Central-difference doesn't help at all. I have changed delta-DE algorithm not to reduce numeric errors but to eliminate artifacts caused by steps of z value. In places where z (in function of fractal coordinates) is not continuous the length(dz) is quite high. Then in those places calculated estimated distance is very low. Ray-marching algorithm sees in those places fake objects (lower distance than threshold) so we have artifacts. When I use minimum value of two directions then the algorithm uses more proper value of delta. Even if there is a step of z value, one of deltas has to be not in the place where is this step (will be just behind or just in front of it).
In attached file there is a result with central-difference

[lycium]

That's very interesting, thanks for the render Buddhi! And good catch Syntopia, I didn't look closely enough at the not-central-differences

[Roquen]

My thoughts are always worth about what you paid for the them...and I have very little (read none) experience at the real topic at hand.  My thinking wasn't really about precision for precision sake, but as a potential means to increase robustness of the computation and if by chance computational complexity could be reduced...that would be added bonus.

An example to be slightly clearer to a wider audience.  If you have delta = 1e-6f, then any value X with |X|>= 16, X + delta = X.  And conversely if X is small, then delta becomes quite large.  If |X| is on [8,16) then the lowest bit will be changed, [4,8) the second to lowest, etc.  Now it's very easy to formula simple equations which burn through valid bits, so just because the input has been effected by the sum (or difference) doesn't imply that the returned result isn't more noise than signal.