Newton Fractal Smoothing

[fodinabor]

Hey Guys,

I'm writing myself a Newton Fractal Generator. The images all formerly looked like this:

But I actually would like it to look a bit smoother - so I've done some research and I ran over http://www.hiddendimension.com/FractalMath/Convergent_Fractals_Main.html and some threads here and no it looks nearly fine, except that there are some issues at the edges of the basins..


This is my generation loop (it's ran in a OpenCL Kernel):

Code:

    while (i < 6000 && fabs(z.r) < 10000 && !found){
		f = computeFunction(z, params, paramc[0]);
		d = computeFunction(z, paramsD, paramc[1]);

		iterexp = iterexp + exp(-fabs(z.r) - 0.5 / (fabs(subComplex(zo, z).r)));

		zo = z;

		z = subComplex(z, divComplex(f, d));
		
		i++;
		
		for (int j = 0; j < paramc[0] - 1; j++){
			if (compComplex(z, zeros[j], RESOLUTION)){
				resType[x + xRes * y] = j;
				result[x + xRes * y] = iterexp;
				found = true;
				break;
			}
		}

		if (compComplex(z, zo, RESOLUTION/100)){
			resType[x + xRes * y] = 12;
			break;
		}
	}

The coloration:

Code:

	const int xRes = res[0];
	const int yRes = res[1];

	for (int y = 0; y < fraktal->getHeight(); y++){
		for (int x = 0; x < fraktal->getWidth(); x++){
			int type, it;
			double conDiv;
			if (genCL && genCL->err == CL_SUCCESS){
				conDiv = genCL->result[x + y * xRes];
				type = genCL->typeRes[x + y * xRes];
				it = genCL->iterations[x + y * xRes];
			} else {
				type = 3;
				conDiv = runNewton(std::complex<double>((double)((x - (double)(xRes / 2)) / zoom[0]), (double)((y - (double)(yRes / 2)) / zoom[1])), type);
			}

			if (type < 15){
				Color col;
				col.setColorHexRGB(colors[type]);
				col.setColorHSV(col.getHue(), col.getSaturation(), 1-conDiv);
				fraktal->setPixel(x, y, col);
			} else {
				fraktal->setPixel(x, y, conDiv, conDiv, conDiv, 1);
			}
		}
	}

The full source code: http://git.war-of-universe.com/fodinabor/newton-fraktale

I appreciate any help to actually smooth this ;-)  
Thanks,
- fodinabor

[claude]

I'd not seen the method you linked before, I'll have to try it sometime.  I use the ones described here:
http://www.chiark.greenend.org.uk/~sgtatham/newton/

a good ad-hoc approximation is obtained simply by looking at the last iteration, in which the point first comes within the specified distance of a root. We look at the distance D0 from the previous point to the root and the distance D1 from the new point to the root, and we know that the threshold distance T is somewhere in between the two. I've found that simply looking at (log T-log D0)/(log D1-log D0), in other words whether the log of the threshold radius was near to the start or the end of the inward distance travelled by the point (on a logarithmic scale), produces a perfectly acceptable result

You either add or subtract that expression from the integer iteration count, can't remember which.

Here's my GLSL implementation:
http://code.mathr.co.uk/fractaloids/blob/HEAD:/webgl/fractaloids.html#l39

Some example videos using variations of that code:
http://mathr.co.uk/blog/2012-12-24_fractal_juggling.html

You might also be interested in distance estimation, here's a post about that too:
http://mathr.co.uk/blog/2013-06-22_distance_estimation_for_newton_fractals.html

[fodinabor]

I translated: (log T-log D0)/(log D1-log D0) into (log(RESOLUTION) - log(fabs(subComplex(zo, zeros[j]).r))) / (log(fabs(subComplex(z, zeros[j]).r)) - log(fabs(subComplex(zo, zeros[j]).r)));

zo is the old z, and zeros[j] is the root the Newton Method is converging to, .r is the abs value (sqrt(t.re*t.re + t.im*t.im))

[claude]

looks good - now you just need to add the iteration count and adjust the colouring range - something like

Code:

col.setColorHSV(col.getHue(), col.getSaturation(), 1- (it + conDiv)/maxIters);

[fodinabor]

thanks!

The coloration now looks like this:

Code:

int maxIters = 0;
	for (int y = 0; y < fraktal->getHeight(); y++){
		for (int x = 0; x < fraktal->getWidth(); x++){
			if (genCL->typeRes[x + y * xRes] < 14){
				maxIters = MAX(maxIters, genCL->iterations[x + y * xRes]);
			}
		}
	}

	maxIters = maxIters / polynom->getNumCoefficients()+1;

	for (int y = 0; y < fraktal->getHeight(); y++){
		for (int x = 0; x < fraktal->getWidth(); x++){
			int type, it;
			double conDiv;
			if (genCL && genCL->err == CL_SUCCESS){
				conDiv = genCL->result[x + y * xRes];
				type = genCL->typeRes[x + y * xRes];
				it = genCL->iterations[x + y * xRes];
			} else {
				type = 3;
				conDiv = runNewton(std::complex<double>((double)((x - (double)(xRes / 2)) / zoom[0]), (double)((y - (double)(yRes / 2)) / zoom[1])), type);
			}

			if (type < 15){
				Color col;
				col.setColorHexRGB(colors[type]);
				conDiv = conDiv + (double)it;
				if (conDiv > maxIters)
					conDiv = maxIters;
				if (conDiv < 0.3)
					conDiv = 0.3;
				col.setColorHSV(col.getHue(), col.getSaturation(), (1.0 - mapCL(conDiv, 0, maxIters, 0.0, 1.0)));
				fraktal->setPixel(x, y, col);
			} else {
				fraktal->setPixel(x, y, 0, 0, 0, 1);
			}
		}
	}

and the result is:

[fodinabor]

hmm... I'd still get the exponential smoothing thing get to work.. it looks a little bit nicer and you don't have to say what to use as maxIters (as it does not smooth anything when using the real max val...

[Pauldelbrot]

Convergent basins are probably best smooth-colored by first finding a point belonging to the attractor and the attractor's period p, then examining every pth step of the orbit of a point that goes to that attractor and:

• If the attractor is superattracting (which occurs with the finite attractors in Newton's method), using the usual divergent smooth coloring applied to 1/z for z the points in the (every pth point) orbit.
• If the attractor isn't superattracting, or whether it is or not varies from point to point (Mandelbrot views do this), then once the orbit is close enough to the attractor it will spiral in in a predictable way: z_{n + p} - z_a = k(z_n - z_a), where z_a is the attracting point this 1/p sub-orbit converges to, and k is a complex constant with |k| < 1. This means that one point in the orbit will land in a target annulus around z_a whose inner diameter is |k| times its outer diameter. The logarithm of the proportion of the way this point landed between the inner diameter and the outer one is the fraction to add to the integer number of iterations (divided by p) taken to hit the annulus. (Thus, if the 6th iteration lands in the annulus, and it just barely does so, landing just inside the outer rim, that's maybe a 6.9; if it lands near the inner rim it's barely above 6.0; if it lands just past the inner rim, the fifth iteration should have landed just inside the outer rim, and we actually have a 5.9. So it's smooth and continuous.)

I use the second method in many of my own images to produce smooth basin coloring for, for instance, Julia set interiors.