> I observe that the surface appears at values smaller than the threshold value if
> it is not reached. To be more precise: I have a function yielding only 0 or 1

Non-continuous functions like that are very tricky for the solver to 
work with. Look up "mandelbulb distance function" and use something like 
that rather than just a zero/one function. It will make the solver much 
more efficient (ie better results in a shorter time).

> depending on the membership to a given set. 0 indicates membership, 1 not (1
> minus the characteristic function of the set). Threshold values of 0.1, 0.2 and
> 0.9 all give the desired object. To my surprise 1.0 too, 1.1 gives the expected
> unit sphere. Threshold 0.0 yields nothing.
>
> May be a numeric issue at threshold 0.0 and 1.0.
>
> What do I miss here?

The solver tries to find a small interval (the size of which depends on 
"accuracy") where one side is below the interval and the other side is 
above the interval. It does this by looking at the function value at the 
first point on the ray, estimating how much further to look ahead based 
on that value, then repeating, until it gets a satisfactory interval.

If your function is always one before the object, the solver must step 
along very slowly taking many samples before it gets to a zero value. 
Whether your "threshold" is 0.0001 or 0.9999 or anywhere in betweween 
won't make any difference. Your function instantly flips from zero to 
one, so any interval will always contain either all the values from zero 
to one or none of them.

However if you use a distance field function, the function value might 
be 100 at the first point, so the solver can take a huge step towards 
the object. And then when it gets very close you'll get fractional 
distances, so the solver can rapidly converge on an accurate interval 
that contains the surface.

See below code from my GPU mandelbulb renderer that calculates the 
distance function at a point, if it's any help:

float MB(vec3 pos)
{
	pos.xyz = pos.xzy;
	vec3 z = pos;
	float dr = 1.0;
	float r = 0.0;
	for (int i = 0; i < 64 ; i++) {
		r = length(z);
		if (r>2.0)
			break;
			
			if(r==0.0)return 0.0;
					
		// convert to polar coordinates
		float theta = acos(z.z/r);
		float phi = atan(z.y,z.x);
		float zr = pow( r,N-1.0);
		
		dr =  zr*N*dr + 1.0;
		
		// scale and rotate the point
		theta = theta*N;
		phi = phi*N;
		
		// convert back to cartesian coordinates
		z = zr*r*vec3(sin(theta)*cos(phi), sin(phi)*sin(theta), cos(theta));
		z+=pos;
	}
	return 0.5*log(r)*r/dr;
}

Post a reply to this message