Error estimation of distance estimators (Mandelbulb improvement found!)

[mermelada]

Hey guys! I'm new

First of all, congrats for the amazing community you built! This place is amazing, I've been lurking for some time.

I'm doing some research on the Mandelbulb and wanted to know if someone did error estimations of the different DE algorithms to the actual surface of the Mandelbulb (or any arbitrary fractal), specifically in the context of sphere tracing.

I haven't found anything like this yet! If there's nothing built, would people be interested in such a tool?

Thanks!

[mclarekin]

Hi and welcome .

Error estimation is an interesting subject that could do with more investigation, so please investigate it more.  The software Mandelbulber, has some plans for more research in reducing error in  distance estimation , but this is not a high priority for us at the moment.

[mermelada]

I see, I'll probably investigate it a little bit more then!

I implemented a simple tool to give me some insights on the Mandelbulb estimator (with some per iteration Julia offset so we get floating bulbs). It's using iq's optimized CPU estimation.

It's an early result and probably has some bugs, but it definitely makes sense: the DE is overestimating inside the inner sphere, generating the overstepping you can see on the left. Red is more error.

I think it is very interesting because it may lead to a distance bound on these cases without being too conservative (like through an overall bias)



What I'm doing is essentially this:

• Sphere trace through the scene
• On each step, I generate 512 samples on the sphere with radius (distance-eps) that the DE computed.
• If I hit anything inside this sphere, it means that my DE is overestimating something.
• I add the proportional error based on the distance of this hit point vs the expected distance.
• Then normalize by some factor, but it doesn't change a lot.

If you guys have any ideas please share them! I'll probably work on this more, I'm sure there's a lot of stuff I'm missing!

[mermelada]

I think I found the solution to the bulb's problem! The trick was in tweaking the derivative.

I made a small shadertoy to show the result: https://www.shadertoy.com/view/MdSBDR

I also had the chance to test my estimator with this technique, and it confirms the better results:



This is the Mandelbulb with an improved estimator, and what I called derivative bias of 10.0. It took more time to render, but can be tweaked.

I'll research more on to why this happens later with more info later!

[3dickulus]

Hello mermelada and welcome to the forum I applied your tweak to the 'classic mandelbulb' in Fragmentarium and it makes a big difference  
I think I must have applied it properly  or it would look like

Code:

#info Mandelbulb Distance Estimator
#define providesInit

#include "MathUtils.frag"
#include "DE-Raytracer.frag"
#group Mandelbulb

// Number of fractal iterations.
uniform int Iterations;  slider[0,9,100]

// Number of color iterations.
uniform int ColorIterations;  slider[0,9,100]

// Mandelbulb exponent (8 is standard)
uniform float Power; slider[0,8,16]

// Bailout radius
uniform float Bailout; slider[0,5,30]
// mermelada's tweak Derivative bias
uniform float DerivativeBias; slider[0,1,2]

// Alternate is slightly different, but looks more like a Mandelbrot for Power=2
uniform bool AlternateVersion; checkbox[false]

uniform vec3 RotVector; slider[(0,0,0),(1,1,1),(1,1,1)]

uniform float RotAngle; slider[0.00,0,180]

uniform bool Julia; checkbox[false]
uniform vec3 JuliaC; slider[(-2,-2,-2),(0,0,0),(2,2,2)]

uniform float time;
mat3 rot;
void init() {
	 rot = rotationMatrix3(normalize(RotVector), RotAngle);
}

// This is my power function, based on the standard spherical coordinates as defined here:
// http://en.wikipedia.org/wiki/Spherical_coordinate_system
//
// It seems to be similar to the one Quilez uses:
// http://www.iquilezles.org/www/articles/mandelbulb/mandelbulb.htm
//
// Notice the north and south poles are different here.
void powN1(inout vec3 z, float r, inout float dr) {
	// extract polar coordinates
	float theta = acos(z.z/r);
	float phi = atan(z.y,z.x);

	// mermelada's tweak
	// http://www.fractalforums.com/new-theories-and-research/error-estimation-of-distance-estimators/msg102670/?topicseen#msg102670
	  dr =  max(dr*DerivativeBias,pow( r, Power-1.0)*Power*dr + 1.0);

	// scale and rotate the point
	float zr = pow( r,Power);
	theta = theta*Power;
	phi = phi*Power;

	// convert back to cartesian coordinates
	z = zr*vec3(sin(theta)*cos(phi), sin(phi)*sin(theta), cos(theta));
}

// This is a power function taken from the implementation by Enforcer:
// http://www.fractalforums.com/mandelbulb-implementation/realtime-renderingoptimisations/
//
// I cannot follow its derivation from spherical coordinates,
// but it does give a nice mandelbrot like object for Power=2
void powN2(inout vec3 z, float zr0, inout float dr) {
	float zo0 = asin( z.z/zr0 );
	float zi0 = atan( z.y,z.x );
	float zr = pow( zr0, Power-1.0 );
	float zo = zo0 * Power;
	float zi = zi0 * Power;

	// mermelada's tweak
	// http://www.fractalforums.com/new-theories-and-research/error-estimation-of-distance-estimators/msg102670/?topicseen#msg102670
	  dr = max(dr*DerivativeBias,zr*dr*Power + 1.0);

	zr *= zr0;
	z  = zr*vec3( cos(zo)*cos(zi), cos(zo)*sin(zi), sin(zo) );
}

// Compute the distance from `pos` to the Mandelbox.
float DE(vec3 pos) {
	vec3 z=pos;
	float r;
	float dr=1.0;
	int i=0;
	r=length(z);
	while(r<Bailout && (i<Iterations)) {
		if (AlternateVersion) {
			powN2(z,r,dr);
		} else {
			powN1(z,r,dr);
		}
		z+=(Julia ? JuliaC : pos);
		r=length(z);
		z*=rot;
		if (i<ColorIterations) orbitTrap = min(orbitTrap, abs(vec4(z.x,z.y,z.z,r*r)));
		i++;
	}

	return 0.5*log(r)*r/dr;
}

#preset Default
FOV = 0.62536
Eye = 1.65826,-1.22975,0.277736
Target = -5.2432,4.25801,-0.607125
Up = 0.401286,0.369883,-0.83588
EquiRectangular = false
FocalPlane = 1
Aperture = 0
Gamma = 2.08335
ToneMapping = 3
Exposure = 0.6522
Brightness = 1
Contrast = 1
Saturation = 1
GaussianWeight = 1
AntiAliasScale = 2
Detail = -2.84956
DetailAO = -1.35716
FudgeFactor = 1
MaxRaySteps = 164
Dither = 0.51754
NormalBackStep = 1
AO = 0,0,0,0.85185
Specular = 1.6456
SpecularExp = 16.364
SpecularMax = 10
SpotLight = 1,1,1,1
SpotLightDir = 0.63626,0.5
CamLight = 1,1,1,1.53846
CamLightMin = 0.12121
Glow = 1,1,1,0.43836
GlowMax = 52
Fog = 0
HardShadow = 0.3538500
ShadowSoft = 12.5806
Reflection = 0.0
DebugSun = false
BaseColor = 1,1,1
OrbitStrength = 0.14286
X = 1,1,1,1
Y = 0.345098,0.666667,0,0.02912
Z = 1,0.666667,0,1
R = 0.0784314,1,0.941176,-0.0194
BackgroundColor = 0.607843,0.866667,0.560784
GradientBackground = 0.3261
CycleColors = false
Cycles = 4.04901
EnableFloor = false
FloorNormal = 0,0,0
FloorHeight = 0
FloorColor = 1,1,1
Iterations = 12
ColorIterations = 8
Power = 8
Bailout = 6.279
AlternateVersion = true
RotVector = 1,1,1
RotAngle = 0
Julia = false
JuliaC = 0,0,0
#endpreset

#preset Octobulb
FOV = 0.62536
Eye = -0.184126,0.843469,1.32991
Target = 1.48674,-5.55709,-4.56665
Up = 0,1,0
AntiAlias = 1
Detail = -2.47786
DetailAO = -0.21074
FudgeFactor = 1
MaxRaySteps = 164
BoundingSphere = 2
Dither = 0.5
AO = 0,0,0,0.7
Specular = 1
SpecularExp = 27.082
SpotLight = 1,1,1,0.94565
SpotLightDir = 0.5619,0.18096
CamLight = 1,1,1,0.23656
CamLightMin = 0.15151
Glow = 0.415686,1,0.101961,0.18421
Fog = 0.60402
HardShadow = 0.7230800
Reflection = 0.0
BaseColor = 1,1,1
OrbitStrength = 0.62376
X = 0.411765,0.6,0.560784,-0.37008
Y = 0.666667,0.666667,0.498039,0.86886
Z = 0.666667,0.333333,1,-0.25984
R = 0.4,0.7,1,0.36508
BackgroundColor = 0.666667,0.666667,0.498039
GradientBackground = 0.5
CycleColors = true
Cycles = 7.03524
FloorNormal = 0,0,0
FloorHeight = 0
FloorColor = 1,1,1
Iterations = 14
ColorIterations = 6
Power = 8.18304
Bailout = 6.279
AlternateVersion = true
RotVector = 1,0,0
RotAngle = 77.8374
#endpreset

edit: btw thank you

[mermelada]

Awesome! I'm glad it worked, your implementation seems correct. I think this is a simple trick that can be easily integrated in all fractal software out there.

The trick lies in keeping the maximum of the derivatives when iterating in the DE function. From your example,

dr = max(dr*DerivativeBias, pow(r, Power-1.0) * Power * dr + 1.0);
where DerivativeBias can be any number >= 1, although the higher the number, the shorter the steps (and more iterations/time).

It seems to work better than just multiplying the whole DE by a number < 1.0, because it adapts to these fractal bubbles of overstepping.
I have some intuition on why this happens (I think dr converges to 0 on those places), but I'll try to work on a proof later.

[mclarekin]

I tried this yesterday in old Mandelbulber OpenCL v1.21.  From zoomed out it seems to work  similarly to lowering the detail level but I did not do  much. testing.   I will try it later in MandelbulberV2 where I can see the statistics of Percentage of Wrong Distance estimations.

I also have been wondering  if there might be any benefit in varying DerivativeBias during the iteration.   e.g  DerivativeBase  =  DerivativeBase  + or *   (a  parabolic or some compounding function) hmmm??


Thanks and keep up the good work