MATLAB - Image Processing - Douday's rabbit fractal

[mheslep]

For anyone interested, here's the same code in Python using the third party Numpy and Matplotlib modules. Linux or windows or Mac of course. I read that Numpy runs some orders of magnitude faster than Matlab and of course Python/Numpy requires no license fee.

import numpy as np
import pylab
import math

iterations = 80
N=500

x,y = np.meshgrid(np.linspace(-1.5, 1.5, N), np.linspace(-1.5, 1.5, N))

c = -0.123 + 0.745j
z = x+y*1j
map = np.zeros( (N,N))
sqrt5 = math.sqrt(5)

for k in range(iterations):
    print "iteration:",k
    z = z**2 + c
    a = np.where(np.abs(z)<sqrt5)
    map[a] = k

pylab.imshow(map)
pylab.show()

In addition to interactive sessions, the code above can be run with the debugging module pdb for line by line execution and inspection:

import pdb
...
#run when debugging desired
pdb.set_trace()

 

[uart]

For anyone interested, here's the same code in Python

Thanks for the info mheslep. Very interesting, I'll have to take a look at those python modules. In the past I know it's been common to translate MATLAB to fortran for faster execution. I have actually translated the above code in Fortran and it is indeed several times faster than the MATLAB code (though i did it more as a sanity check in this case) BTW. Also take a look at gnu-octave for an excellent freeware MATLAB clone.Back to the question of why greenprint is getting unexpected results. Would you be able to make a quick test for us with that python code to test the conjecture that the direction of the inequality in a = np.where(np.abs(z)<sqrt5) should cause either a multi-level or a
two-level map to be produced.

In particular could you please test if replacing the above line with, a = np.where(np.abs(z)>sqrt5), causes your program to produce only a two level map (and a corresponding monochrome image).

 

[mheslep]

...
Back to the question of why greenprint is getting unexpected results. Would you be able to make a quick test for us with that python code to test the conjecture that the direction of the inequality in a = np.where(np.abs(z)<sqrt5) should cause either a multi-level or a
two-level map to be produced.

In particular could you please test if replacing the above line with, a = np.where(np.abs(z)>sqrt5), causes your program to produce only a two level map (and a corresponding monochrome image).

Sure. The orginal, abs(z)<sqrt(5):

and the almost(?) monochrome image obtained from
abs(z)>sqrt(5):

 

[uart]

abs(z)>sqrt(5):

Arrh. That's exactly the image that Greenprint is getting and that we're having difficulty explaining. You see that's actually far from a two level map, there are lots of different shades and colors in that image, which implies that the map "matrix" must contain may different numbers, but I can only explain two of them!

Anyway I'm currently installing those python packages and am just about to try and reproduce you code. Hopefully I'll be able to get to the bottom of it once I see the results for myself. :)

 

[uart]

Mystery solved !

Ok once I got the python code running it didn't take long to figure out what was happening. In the end it turned out to be fairly simple, the gnu-Octave I'm using thinks infinity > sqrt(5), as do I, but Matlab and Python think otherwise.

Well actually it's the way the respective programs handle numerical overflow. There are regions on the map that begin with |z|>2, so not surprisingly when you keep repetitively squaring this you soon get a numerical overflow, (((2^2)^2)^2)^…, it only takes about 10 iterations to exceed even a double precision floating point capability. When this happens Octave returns "+inf" for the absolute value, while MATLAB and python return "nan" (not a number).

So that's it. All the lovely shading that Greenprint (and his MATLAB textbook) have been getting in their images have been purely the result of bad programming and a glitch in the way the MATLAB handles numerical overflow.

 

[mheslep]

Ah, thanks.

Numpy has a hack around for the overflow that (efficiently) forces NaN's to zeros:
http://www.scipy.org/Numpy_Example_List_With_Doc?highlight=%28NaN%29#head-c98ac710ae88aadee85e953af821e560ab316ef3

So that

   ...
    z = z**2 + c
    z = np.nan_to_num(z)
    a = np.where(np.abs(z)>sqrt5)
    map[a] = k
    ...

produces the attached image as you would expect

 

[mheslep]

Python-Numpy comparison with Matlab:
http://www.scipy.org/NumPy_for_Matlab_Users