Print Page - Mandelbuld render glitch

Title: Mandelbuld render glitch
Post by: lapinot on February 10, 2015, 04:14:11 PM

Hello there.

I am programming a raymarcher in C++. It is still 100% CPU but planning to move to OpenGL as soon as I can.

So far I tried rendering a mandelbulb (variation) using Quilez formula (http://www.iquilezles.org/www/articles/mandelbulb/mandelbulb.htm - the non-optimized one) but can't figure out where to look exactly to solve this problem.

Here is a small render of it - please note the red color is the background/sky color.

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423476821.png)

It looks like some kind of shiva/elephant to me - using a power = 4 (code below). The only problem being, most of the fractal surface isn't displayed at all but gets overshooten. I guess it should look something like :

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423476821-bounded.png)

with the area approximately inside the green bound fully drawn. Worth mentioning too : the white pixels are actually shaded using the raymarching steps ambiant occlusion method. They "appear" white because for some reason, all the "white" pixels are intersected at either step 1 or 2. Any step more and I get geometry overshoot and no intersection at all.

So I dumped the raymarch steps of a random "red" pixel and I get :

Code:

#STEP 0# x=0 y=1 z=-15
#DISTANCE# 17.8838

#STEP 1# x=-0.107225 y=-2.16106 z=2.60185
#DISTANCE# 3.00974

#STEP 2# x=-0.12527 y=-2.69305 z=5.56415
#DISTANCE# 9.87785

#STEP 3# x=-0.184495 y=-4.43902 z=15.2863
#DISTANCE# 24.4502

#STEP 4# x=-0.33109 y=-8.76075 z=39.3511
#DISTANCE# 74.5188

#STEP 5# x=-0.777879 y=-21.9324 z=112.695
#DISTANCE# 272.294

#EXIT#

The weird thing is at Step 0 the distance is already 17.8 when the point of view is only at z=-15 so it looks to me the fractal geometry got overshot. I get the same kind of results with a "real" power 8 mandelbulb formula.

Is my DE function not correct ?

Code:

long double DistanceEstimators::CPPDE (D3Vec w) {

	D3Vec z = w;
	float r = z.getLength();
	float dr = 1.0;

	int escape = 20;
	int Power = 4;
	double Scale = 1;

	// Iterate
	for (int i = 0; i != 10 && r < escape * Scale; i++) {

		// extract polar coordinates
		float wr = std::sqrt(w * w);
		float wo = std::acos(w.y / wr);
		float wi = std::atan2(w.x, w.z);

		// scale and rotate the point
		wr = std::pow(wr, Power);
		wo = wo * Power;
		wi = wi * Power;

		// convert back to cartesian coordinates
		w.x = std::sin(wo) * std::sin(wi);
		w.y = std::cos(wo);
		w.z = std::sin(wo) * std::cos(wi);

		z = z + w;
		r = Scale * z.getLength();
	}

	// Back
	return 0.5 * std::log(r) * r / dr;
}

Thanks for any hint!

-Lapinot

Title: Re: Mandelbuld render glitch
Post by: eiffie on February 10, 2015, 06:58:12 PM

2 things:
1) the distance estimate for bulbs is inaccurate at long distances so just use something like DE=min(DE,1.0) or start your ray closer to the bulb.
2) I don't see where you set dr to anything? just dr=1.0
dr is the running derivative so for power 4 it looks something like...
dr=dr*4.0*pow(r,3.0)+1.0;
..inside the loop

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 11, 2015, 12:19:06 PM

I added this dr computation line in my code :

Code:

	// Iterate
	for (int i = 0; i != 10 && r < escape * Scale; i++) {

		//if (i > 0) std::cout << i << std::endl;

		// extract polar coordinates
		float wr = std::sqrt(w * w);
		float wo = std::acos(w.y / wr);
		float wi = std::atan2(w.x, w.z);
		//float wi = std::atan2(w.z, w.x); // <-- test

		// scale and rotate the point
		wr = std::pow(wr, Power);
		wo = wo * Power;
		wi = wi * Power;

		// convert back to cartesian coordinates
		w.x = std::sin(wo) * std::sin(wi);
		w.y = std::cos(wo);
		w.z = std::sin(wo) * std::cos(wi);
		
		// running derivative
		dr = std::pow(r, Power - 1) * Power * dr + 1;

		z = z + w;
		r = Scale * z.getLength();
	}

Something is wrong. My DR ends up reaching infinity at step 10. Here is a sample :

Code:

step #0: 13591.1
step #1: 1.98131e+008
step #2: 2.46404e+012
step #3: 3.6136e+016
step #4: 4.3393e+020
step #5: 4.47056e+024
step #6: 5.11905e+028
step #7: 6.40744e+032
step #8: 9.51215e+036
step #9: inf

Of course my returned distance is then 0 for all DE evals. Anything obvious ? I didn't change my PoV (x=0 y=1 z=-15).

Title: Re: Mandelbuld render glitch
Post by: knighty on February 11, 2015, 04:28:13 PM

Try this: (not tested)

Code:

double DistanceEstimators::CPPDE (D3Vec w) {

	D3Vec z = w;
	double wr = w.getLength();
	double dr = 1.0;

	double escape = 20;
	int Power = 4;

	// Iterate
	for (int i = 0; i < 10 && wr < escape ; i++) {

		// extract polar coordinates
		
		double wo = std::acos(w.y / wr);
		double wi = std::atan2(w.x, w.z);

		dr= dr * Power * std::pow(wr, Power-1) + 1.;//"scalar" running derivative update
		
		// scale and rotate the point
		wr = std::pow(wr, Power);
		wo = wo * Power;
		wi = wi * Power;

		// convert back to cartesian coordinates
		w.x = wr * std::sin(wo) * std::sin(wi);
		w.y = wr * std::cos(wo);
		w.z = wr * std::sin(wo) * std::cos(wi);

		w = w + z; 
		wr = w.getlength();
	}

	// Back
	
	return 0.5 * std::log(wr) * wr / dr;
}

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 13, 2015, 10:58:26 AM

I get the same problem with your DE Knighty.

For the same sample point I get :

Code:

#STEP 0# x=0 y=1 z=-15
#DISTANCE# 20.3728

#STEP 1# x=-0.122148 y=-2.60101 z=5.05162
#DISTANCE# 4.96099

#STEP 2# x=-0.151893 y=-3.47789 z=9.93441
#DISTANCE# 12.3968

#STEP 3# x=-0.226219 y=-5.66909 z=22.1357
#DISTANCE# 35.7509

#STEP 4# x=-0.44057 y=-11.9883 z=57.3231
#DISTANCE# 119.184

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 13, 2015, 06:33:31 PM

Bump I need help :(

Here is a sample power 8 bulb:

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423848275-bulb.png)

It definitely looks like the bulb to me, except it's not filled ofc.

I used Quilez formula and no DR. When computing a DR I get:

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423848640-bulb-dr.png)

a full white image because the DR gets infinite so the returned distance is 0 for all DE evals.

Code below for the first image:

Code:

long double DistanceEstimators::CPPDE (D3Vec w) {

	D3Vec z = w;
	float r = z.getLength();
	float dr = 1.0;

	int escape = 20;
	int Power = 8;
	double Scale = 1; // or 0.5

	// Iterate
	for (int i = 0; i != 10 && r < escape * Scale; i++) {

		// extract polar coordinates
		float wr = std::sqrt(w * w);
		float wo = std::acos(w.y / wr);
		float wi = std::atan2(w.x, w.z);

		// scale and rotate the point
		wr = std::pow(wr, Power);
		wo = wo * Power;
		wi = wi * Power;

		// convert back to cartesian coordinates
		w.x = std::sin(wo) * std::sin(wi);
		w.y = std::cos(wo);
		w.z = std::sin(wo) * std::cos(wi);

		z = z + w;
		r = Scale * z.getLength();
	}

	// Back
	return 0.5 * std::log(r) * r / dr;
}

I think the DR is the culprit here and it shouldn't reach infinity so quick ? The DR I used in the second image is:

Code:

dr = std::pow(r, Power - 1) * Power * dr + 1;

Title: Re: Mandelbuld render glitch
Post by: eiffie on February 13, 2015, 07:27:22 PM

It is getting worse! The last bit of code posted is backwards. You are doing the triplex math on w then adding it to z. Go back to adding z to w like the previous code. The distance estimate is bad at long distances (look at your steps - you step right over it, that is why dr goes to infinity). Do something like this...
return min(0.5 * std::log(wr) * wr / dr,1.0);

...and post the complete code using dr. We all know it won't work without a running derivative.

You'll get it! oh and drop the scale until you get it working at scale 1.0 :)

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 13, 2015, 07:31:21 PM

Thx. Working on it ill keep you updated :) Funny thing tho, I still get a mandelbulb-like render with the wrong math :alien:

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 13, 2015, 08:31:56 PM

OK I reverted back + fixed to:

Code:

long double DistanceEstimators::CPPDE (D3Vec w) {

	D3Vec z = w;
	float wr = w.getLength();
	float dr = 1.0;

	int escape = 20;
	int Power = 4;

	// Iterate
	for (int i = 0; i != 10 && wr < escape; i++) {

		// extract polar coordinates
		float wo = std::acos(w.y / wr);
		float wi = std::atan2(w.x, w.z);
		dr = std::pow(wr, Power - 1) * Power * dr + 1;

		// scale and rotate the point
		wr = std::pow(wr, Power);
		wo = wo * Power;
		wi = wi * Power;

		// convert back to cartesian coordinates
		w.x = std::sin(wo) * std::sin(wi);
		w.y = std::cos(wo);
		w.z = std::sin(wo) * std::cos(wi);
		
		// add starting point to new point
		w = w + z;
		wr = w.getLength();
	}
	
	// Back
	return std::min(1.0, 0.5 * std::log(wr) * wr / dr);
}

I still get the infinite DR aka DE evals = 0 :(
http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/DR_sample.png (http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/DR_sample.png)

Title: Re: Mandelbuld render glitch
Post by: knighty on February 14, 2015, 03:20:37 PM

You forgot to multiply w by wr in "convert back to cartesian coordinates" section. As eiffie said the DE obtained by using 0.5*r*log(r)/dr is always too big when far away from the fractal.
You are also using floats inside the function. better use doubles for now because it doesn't overflow (giving infinty) as easily as floats. This may happen when the point where the DE is computed is very close to the surface of the fractal. Could we see the raymarching code? Have you tried your raymarching code with a simple object (say a sphere: DE=w.length()-radius;)?
Beside this I don't understand why your function returns long double.

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 14, 2015, 04:48:56 PM

Yes I can render scene from a long time ago.

For example sphere + ground

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423928350.png)

I return long doubles to get as much precision as I can. I just pretty much copied the DE code from around the web (Quilez/Hvidfeldts) hence the floats. I can as well use long doubles instead I guess.

As you can see from the picture above, I think my raymarching code is OK. I'm not keen on posting it all there tho ;D

Should I just do a scalar multiply of w by wr ? I don't pretend to understand the math here hehe.

Code:

		// convert back to cartesian coordinates
		w.x = std::sin(wo) * std::sin(wi);
		w.y = std::cos(wo);
		w.z = std::sin(wo) * std::cos(wi);
		w *= wr; // ?

Thanks.

Title: Re: Mandelbuld render glitch
Post by: eiffie on February 14, 2015, 06:11:12 PM

Yes multiply that way. It seems like that is all you are missing.

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 15, 2015, 12:46:05 AM

Sooooooooooooooooooo. Here are my first renders of the bulb!!

Thx a lot guys for the help so far ;D

Too high DE bailout
(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423954752.png)

Smaller DE bailout
(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423954833.png)

Double res.
(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423954995.png)

Close-up view
(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423955632.png)

Lightning test
(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1423956071.png)

Need to work on the lighting ;D

Title: Re: Mandelbuld render glitch
Post by: cKleinhuis on February 15, 2015, 01:47:16 AM

congrats, wonderful ;) now for the colors,lighting and stuff ;) say hello to pandorras box you just opened!

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 15, 2015, 03:55:07 PM

Thx ;D

I have no idea about the colouring yet. I did implement lighting tho but need to retwink it. The last image is with lighting but is way too dark.

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 15, 2015, 05:32:31 PM

I just tried to improve the rendering speed using Quilez optimization and I get:

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1424017245.png)

The (working) code I used for the images above was:

Code:

long double DistanceEstimators::CPPDE (D3Vec w) {

	D3Vec z = w;
	float wr = w.getLength();
	float dr = 1.0;

	int escape = 20;
	int Power = 4;
	//double Scale = 1;

	// Iterate
	for (int i = 0; i != 10 && wr < escape; i++) {

		// extract polar coordinates
		//float wr = std::sqrt(w * w);
		float wo = std::acos(w.y / wr);
		float wi = std::atan2(w.x, w.z);
		dr = std::pow(wr, Power - 1) * Power * dr + 1;

		// scale and rotate the point
		wr = std::pow(wr, Power);
		wo = wo * Power;
		wi = wi * Power;

		// convert back to cartesian coordinates
		w.x = std::sin(wo) * std::sin(wi);
		w.y = std::cos(wo);
		w.z = std::sin(wo) * std::cos(wi);
		w = w * wr;
		
		// add starting point to new point
		w = w + z;
		wr = w.getLength();
		
		// running derivative
		//std::cout << "#step " << i << ": " << dr << std::endl;
	}
	
	// Back
	return std::min(1.0, 0.5 * std::log(wr) * wr / dr);
}

And the one with Quilez optimization:

Code:

long double DistanceEstimators::CPPDE (D3Vec w) {

	D3Vec c = w;
	long double wr = w.getLength();
	long double dr = 1.0;

	int escape = 20;
	int Power = 4;
	//double Scale = 1;

	// Iterate
	for (int i = 0; i != 10 && wr < escape; i++) {
	
		// running derivative
		dr = std::pow(wr, Power - 1) * Power * dr + 1;
		
		// quilez optimization
		long double x = w.x; long double x2 = x*x; long double x4 = x2*x2;
		long double y = w.y; long double y2 = y*y; long double y4 = y2*y2;
		long double z = w.z; long double z2 = z*z; long double z4 = z2*z2;

		long double k3 = x2 + z2;
		long double k2 = 1.0L / std::sqrt(k3*k3*k3*k3*k3*k3*k3);
		long double k1 = x4 + y4 + z4 - 6.0*y2*z2 - 6.0*x2*y2 + 2.0*z2*x2;
		long double k4 = x2 - y2 + z2;

		w.x =  64.0*x*y*z*(x2-z2)*k4*(x4-6.0*x2*z2+z4)*k1*k2;
		w.y = -16.0*y2*k3*k4*k4 + k1*k1;
		w.z = -8.0*y*k4*(x4*x4 - 28.0*x4*x2*z2 + 70.0*x4*z4 - 28.0*x2*z2*z4 + z4*z4)*k1*k2;
	
		// scale w
		wr = std::pow(wr, Power);
		w = w * wr;
		
		// add starting point to new point
		w = w + c;
		wr = w.getLength();
	}
	
	// Back
	return std::min(1.0L, 0.5 * std::log(wr) * wr / dr);
}

I basically tried to insert the modification into the previous code. And checked that my vars are initialized in the correct order regarding previous code but there's a mistake.

Title: Re: Mandelbuld render glitch
Post by: eiffie on February 15, 2015, 06:45:13 PM

Glad you got it working! I would stick with the version that lets you adjust the power easily. When you switch from CPU to GPU it will be plenty fast.

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 15, 2015, 09:43:33 PM

Yes you're right.. Actually I even mixed Powers 4 and 8 in the formula. With only power 8 I still get the same kind of image tho. I'll save it for later maybe :)

Title: Re: Mandelbuld render glitch
Post by: lapinot on February 17, 2015, 06:24:17 PM

Trying to improve my raymarcher's escape condition and I got:

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1424193130.png)

(http://grid.turtlespeak.net/upload/LMo6pu2bZ3pYC5QU/OUmytkQpaGkM6C5e/Render-1424193703.png)

Funny pics ;D

Welcome to Fractal Forums

Fractal Software => Programming => Topic started by: lapinot on February 10, 2015, 04:14:11 PM