News: Check out the originating "3d Mandelbulb" thread here  ## The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!

 Pages:    Go Down       Author Topic: Improved delta-DE algorithm  (Read 2446 times) Description: 0 Members and 1 Guest are viewing this topic.
Buddhi
Fractal Iambus   Posts: 895   « on: August 15, 2013, 03:38:02 PM »

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.  old.jpg (81.39 KB, 400x300 - viewed 679 times.)  inproved.jpg (50.42 KB, 400x300 - viewed 711 times.) Logged

Roquen
Iterator Posts: 180 « Reply #1 on: August 15, 2013, 05:19:30 PM »

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 Logged

All code submitted by me is in the public domain. (http://unlicense.org/)
eiffie
Guest « Reply #2 on: August 15, 2013, 05:41:10 PM »

Thanks for this improvement! Logged
lycium
Fractal Supremo     Posts: 1158   « Reply #3 on: August 15, 2013, 08:15:04 PM »

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!  Logged

Roquen
Iterator Posts: 180 « Reply #4 on: August 15, 2013, 08:40:17 PM »

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. Logged

All code submitted by me is in the public domain. (http://unlicense.org/)
Syntopia
Fractal Molossus  Posts: 681    « Reply #5 on: August 15, 2013, 10:13:46 PM »

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! 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 Logged
Roquen
Iterator Posts: 180 « Reply #6 on: August 15, 2013, 10:22:06 PM »

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. Logged

All code submitted by me is in the public domain. (http://unlicense.org/)
Buddhi
Fractal Iambus   Posts: 895   « Reply #7 on: August 16, 2013, 06:59:50 PM »

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  central-difference.jpg (118.9 KB, 400x300 - viewed 616 times.) Logged

lycium
Fractal Supremo     Posts: 1158   « Reply #8 on: August 16, 2013, 07:33:33 PM »

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

Roquen
Iterator Posts: 180 « Reply #9 on: August 16, 2013, 08:25:40 PM »

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. Logged

All code submitted by me is in the public domain. (http://unlicense.org/)
 Pages:    Go Down
 Related Topics Subject Started by Replies Views Last post  Delta DE limits 3D Fractal Generation David Makin 5 2416 October 24, 2009, 05:35:30 PM by David Makin   Delta Sunset Mandelbulb3D Gallery lenord 0 835 December 17, 2010, 01:12:15 AM by lenord  Delta Images Showcase (Rate My Fractal) Kali 6 1047 June 27, 2012, 06:28:37 AM by Kali  Hybri-Station on Delta 4 Images Showcase (Rate My Fractal) JoeFRAQ 1 723 October 31, 2014, 02:09:53 AM by JohnVV  Orange Delta Images Showcase (Rate My Fractal) C.K. 4 680 January 11, 2015, 06:46:21 PM by C.K.