It is the newton method, z = z - f(z) / f'(z)

where f(z) = cos(z)

f'(z) = -sin(z)

Every point on the plane was rotated 45 degrees.

z(0) = rotate(z(0), 0 + 0i, 45) , where 0 + 0i is the center of the rotation

Next every point was transformed with the circle inversion.

circle_inversion(z, center, radius) {

distance = norm(z - center)^2

radius2 = radius * radius

z = [re(center) + (re(z) - re(center)) * radius2 / distance] + [im(center) + (im(z) - im(center)) * radius2 / distance]*i

return z

}

for the current image every z(0) = circle_inversion(z(0), 1 + 1i, 3)

Next every point is iterated through the normal newton method process, until it converges to a root.

The escape criterion is norm(z - zold) <= convergent_bailout (zold is the previous value)

The coloring algorithm is a mix of the standard smooth escape time algorithm for converging fractals, along with a formula that tries to assign a different number depending on the final z value (the root).

For the smooth escape time algorithm I used:

power = log(norm(z - zold)) / log(norm(zold - zold2))

f = iterations + log(log(convergent_bailout) * 0.5 / log(norm(z - zold))) / log(power)

(zold2 is the second previous value)

For the root extraction I used the phase and the squared magnitude of the final value along with some scaling factors (I had to round the final re and im values because of small variations):

temp = floor(1000 * re(z) + 0.5) / 1000

temp2 = floor(1000 * im(z) + 0.5) / 1000

f2 = ((atan2(temp2, temp) / (2 * pi) + 0.75) * 59 * pi) + (temp * temp + temp2 * temp2) * 2.5

59 * pi is just a random scale factor, you can use whatever you want

This formula obviously does not uniquely define every root of a function (you can have almost the same value for different roots)

but it works for my purpose.

the final number is, res = f + f2

Having this floating point number (res) you need to find the coresponding color.

Lets assume that you have a palette with N colors (palette[N])

You need to use as your data

(int)res, palette[(int)res]

(int)(res + 1), palette[(int)(res + 1)]

and you must perform linear interpolation to find the intermediate value res, which as I mentioned was a floating point value.

The interpolation must be performed in every color channel (R,G,B) seperately.

The final image was processed with the bump mapping algorithm, and the fake distance estimation algorithm.

In the attachments you will find the saved settings of the image (Dont forget to disable the boundary tracing algorithm, because it glitches the image when you use bump mapping).