Smooth colouring of convergent fractals

[zenzero-2001]

Hi All

I am trying to implement David Makin's formula (as described here: http://www.fractalforums.com/programming/help!-newton-fractal-smooth-shading/ and in his Ultra Fractal formula file) for smooth colouring of convergent fractals, but I am experiencing some difficulties. Hopefully he is reading this!

I have not yet tried applying it to the Magnet fractals (which are both convergent and divergent), but am currently trying to smoothly colour the Manget and Nova fractals with the auto-power variant of the formula. I can not get it to work at all with the Nova fractal, but with the Newton fractal it works properly with positive low powers (eg 3), but at higher powers the result becomes progressively less smooth. With negative powers, the inside is smooth but the outside is not. I know the formula should work because it does so under all these conditions with Ultra Fractal. Anyway here is some of my code, maybe somebody can see a silly mistake?

Code:

vector<Element>* Newton::calculateNoOrbit(vector<Point>* points)
{
    vector<Point>::const_iterator i = points->begin();

    int iterations;
    double zx;
    double zy;
    double z2x;
    double z2y;
    double z3x;
    double z3y;
    double zx_old = DBL_MAX;
    double zy_old = DBL_MAX;
    double zx_older = DBL_MAX;
    double zy_older = DBL_MAX;
    double zx_delta;
    double zy_delta;
    double denominator;

    double radius;
    double modulus;
    double phi;

    const double exponent2 = exponent - 1.0;
    const double temp1 = 0.5 * exponent;
    const double temp2 = 0.5 * exponent2;

    getCoord(i->x, i->y, zx, zy);

    for(iterations = 0; iterations < maxIterations; ++iterations)
    {
        // Calculate |z - z_old| for bailout check
        zx_delta = zx - zx_old;
        zy_delta = zy - zy_old;

        zx_delta *= zx_delta;
        zy_delta *= zy_delta;

        // Check |z - z_old| against bailout
        if(zx_delta + zy_delta < bailout)
        { break; }

        // Update z_older and z_old
        zx_older = zx_old;
        zy_older = zy_old;
        zx_old = zx;
        zy_old = zy;

        // Calculate z^exponent - 1
        radius = zx * zx + zy * zy;
        modulus = exp(log(radius) * temp1);
        phi = exponent * atan2(zy, zx);

        z3x = modulus * cos(phi) - 1.0;
        z3y = modulus * sin(phi);

        // Calculate exponent * z^(exponent - 1)
        modulus = exp(log(radius) * temp2);
        phi = exponent2 * atan2(zy_old, zx_old);

        z2x = modulus * cos(phi) * exponent;
        z2y = modulus * sin(phi) * exponent;

        // Calculate the division and subtract from z
        denominator = z2x * z2x + z2y * z2y;

        if(denominator == 0.0)
        { break; }

        zx -= (z3x * z2x + z3y * z2y) / denominator;
        zy -= (z3y * z2x - z3x * z2y) / denominator;
    }

    results.clear();

    // Calculate smooth color
    zx_delta = zx - zx_old;
    zy_delta = zy - zy_old;
    double m1 = log(sqrt(zx_delta * zx_delta + zy_delta * zy_delta));

    zx_delta = zx_old - zx_older;
    zy_delta = zy_old - zy_older;
    double m2 = log(sqrt(zx_delta * zx_delta + zy_delta * zy_delta));

    double fractional = iterations + (log(-0.5 * log(bailout)) - log(-0.5 * m1)) / log( m1 / m2 );
    writeResult(i->x, i->y, iterations, fractional);

    return &results;
}

Thanks!

[fractalrebel]

My favorite, which is nowhere as complicated as the algorithm above, is called general smoothing:

initialize with

iterexp = 0
zold = z

loop until bailout:
  iterexp =iterexp + exp(-cabs(z)-0.5/(cabs(zold - z)))
  zold = z

Then use iterexp to color the fractal. It works with both convergent and divergent fractals, so it iwll work for all regions of a magnet fractal

[fractalrebel]

cabs is an UltraFractal function which you may not know. For a complex number c = a + ib

cabs(c) = sqrt(a^2 + b^2)

[fractalrebel]

For those of you who have Ultrafractal the General Smoothing coloring algorithm is in reb.ulb and the implementation in reb5.ucl.

[zenzero-2001]

Thanks fractalrebel, that is a very nice algorithm. I had look at it in Ultra Fractal.

I would still like to get the other smoothing algorithm working however, because it more closely matches the contours of the level sets method.

[David Makin]

Ron actually sort-of gave the answer - in your implementation you are using the full magnitude calculation i.e. actually getting the square roots.

In both Ultra Fractal and Fractint (and some other fractal software) then:

|x+i*y| = x^2+y^2
abs(x+i*y) = abs(x)+i*abs(y)

and as Ron said:

cabs(x+i*y) = sqrt(x^2+y^2)

Whereas in formal math |x+i*y|=abs(x+i*y)=cabs(x+i*y)=sqrt(x^2+y^2)

To check on the above and indeed how most complex stuff is generally implemented in Fractint and Ultra Fractal see the "Trig identities" section here:

http://www.nahee.com/spanky/www/fractint/append_a_misc.html

So in your code just get rid of the square roots - including for the bailout test i.e. such that the bailout tests (zx-zoldx)^2+(zy-zoldy)^2 instead of the square root of this.
The reason for avoiding the square roots is that they are a major slow-down when used

edit: Also taking away the square roots should give you the reason for the "-0.5*log" in my calculations as here we do actually need the log of the square root and of course -0.5*log(x^2) = -log(sqrt(x^2)). (The minus signs being because we're testing convergent in this case)

[zenzero-2001]

Thanks David! I wondered what the 0.5 was for I will have a play with the code tomorrow as I'm in bed now and it is time for sleep.

I noticed that your formula does not seem to work so well for the magnet fractals in Ultra Fractal, especially when adjusting the colour density. Is there a way to mitigate this, perhaps by using the fudge factor?

I am developing an open source cross platform fractal explorer written in Qt, the smooth colour algorithm will be used in it. I hope to release the first version later this year.

[fractalrebel]

Here is another one which keeps the level sets. This also works well with the magnet formula. Set the divergent bailout on the fractal formula to 1e20 and the convergent bailout to 1e-20. The coloring bailout can be set to either 1e20 or 1e-20 - both work.

Harkonen Vespstas smoothing

loop until bailout:
    zold2 = zold
    zold = z

After bailout
    power = log(|z - zold|) / log(|zold - zold2|)
    f = (log(log(@bailout)) - log(log(1/(|(zold - z)|))))/log(power)

    color = 0.05*(iterations+f)

The |complex num| usage is the UltraFractal and Fractint usage, not the standard mathematical usage.

[David Makin]

Thanks David! I wondered what the 0.5 was for I will have a play with the code tomorrow as I'm in bed now and it is time for sleep.

I noticed that your formula does not seem to work so well for the magnet fractals in Ultra Fractal, especially when adjusting the colour density. Is there a way to mitigate this, perhaps by using the fudge factor?

I am developing an open source cross platform fractal explorer written in Qt, the smooth colour algorithm will be used in it. I hope to release the first version later this year.

Depends which formula/s of mine you are starting from - my best attempts for the Magnet formulas are UF5-class formulas in mmf.ulb, see this thread here at ff:

http://www.fractalforums.com/mathematics/convergent-distance-estimation-t1566/

However all my convergent DE colourings exhibit the same "problem" i.e. singularities (at the roots in the Newtons and equivalent locations in the Magnets and other convergent formulas/areas).
I did manage a work around to fix the issue in Newtons for the solutions of plain z^n-v at least where n is an integer, but extending that further was beyond my formal math knowledge !

[David Makin]

I should add that for point attractors (and indeed periodic attractors) then it's always a lot more accurate to test against known value/s for the attractor for a given orbit (pixel) rather than testing against previous values - for instance for a point attractor use:

z-attractor and zold-attractor instead of z-zold and z-zolder

For a standard Magnet formula that would be (1,0).

[David Makin]

Aargh - sorry, just remembered you were talking about level sets - the important thing is that the bailout values must match those used in the main formula even when using the "auto-power" method, AFAIK my UF5 class colouring "Extended Cilia" in mmf.ulb should work for any Magnet formula provided the bailout values are set correctly - note that to use it you should choose the colouring formula from "Standard.ucl" that allows plug-ins and plug-in "MMF Extended Cilia" from mmf.ulb.
Edit: And initially set the Cilia colouring parameter to "Just Iteration" for both divergent and convergent areas

[zenzero-2001]

The Cilia colouring is interesting but a bit complicated!

I have your formula working now for both Magnet fractals, Nova and Newton. It actually works better with my code on the Magnet fractals than it does in Ultra Fractal for some reason. I still have banding on the outside of the Newton fractal for negative powers, but this is because of trying to find the log of zero.

I see that Ron's method is similar to exponential smoothing, so there would be some overhead finding the exponent every iteration. The results are very nice though.

I don't have too many things left to do before releasing the first version of my program, the biggest being the gradient editor and a few bug fixes. Here are some of the features:

MDI Interface
Animated pan and zoom (similar to Fractal eXtreme)
Arbitrary Precision (using MPIR and MPFR)
Smooth colouring
Fullscreen mode
Colour cycling
Julia animation
Currently supported fractals: Mandelbrot, Mandelbrot3, Mandelbrot4, Mandelbrot5, Mandelbrot6, Mandelbrot Inverse, Mandelbrot General Power, Burning Ship, Lambda, Phoenix, Newton General Power, Nova, Magnet1, Magnet2

It utilises solid guessing, orbit detection and SSE2 intrinsics (for Mandelbrot) for speed. Hopefully, it will be a useful and enjoyable program to use.

I know just how much the World needs another fractal explorer.

Thanks again for the help!

John

[David Makin]

Currently supported fractals: Mandelbrot, Mandelbrot3, Mandelbrot4, Mandelbrot5, Mandelbrot6, Mandelbrot Inverse, Mandelbrot General Power, Burning Ship, Lambda, Phoenix, Newton General Power, Nova, Magnet1, Magnet2

It utilises solid guessing, orbit detection and SSE2 intrinsics (for Mandelbrot) for speed. Hopefully, it will be a useful and enjoyable program to use.

I know just how much the World needs another fractal explorer.

Thanks again for the help!

John

Of course if you really want patronage of your software to take off then it needs to have it's own parser/compiler - and now I guess both CPU and GPU support plus built-in support for 3D+ rendering (ray traced, voxel and 2D slices) with export to STL, VRML etc. as well as an animation engine for both 2D and 3D
And of course it should be available for Windows, OSX and Linux at a minimum but ideally also for iOS, Android and Blackberry.

[zenzero-2001]

Of course if you really want patronage of your software to take off then it needs to have it's own parser/compiler - and now I guess both CPU and GPU support plus built-in support for 3D+ rendering (ray traced, voxel and 2D slices) with export to STL, VRML etc. as well as an animation engine for both 2D and 3D wink
And of course it should be available for Windows, OSX and Linux at a minimum but ideally also for iOS, Android and Blackberry.

Yep, I should have all those done by tomorrow evening - no problem!