buddhabrot.py

#!/usr/bin/env python3
import numpy as np
from PIL import Image # Pillow fork of PIL
import time

width, height = 900, 900   # rotated 90 degrees
xmin, xmax = -1.6, 0.8
itsmax = 1500
qual = 5

yscale = (xmax - xmin) * height / width / 2
ymin, ymax = -yscale, yscale + 0.01
y, x = np.ogrid[ymin:ymax:1j*qual*height, xmin:xmax:1j*qual*width]

c = x + 1j*y
z = np.zeros_like(c)
np.warnings.filterwarnings('ignore', '(overflow|invalid)')
starttime = time.time()

for n in range(itsmax):
    mask = (z * z.conjugate()) < 4
    np.putmask(z, mask, z*z + c)

escaped = np.where(~mask)
cvals = c[escaped]
escvals = cvals.copy()

print("Stage1 complete:", time.time() - starttime, "seconds")
starttime = time.time()

y, x = np.ogrid[ymin:ymax:1j*height, xmin:xmax:1j*width]
surface = np.zeros_like(x * y, dtype=np.uint64)

for _ in range(itsmax):
    newvals = escvals ** 2 + cvals
    inside = np.where(abs(newvals < 10))
    escvals = newvals[inside]
    cvals = cvals[inside]
    xpoints = np.uint64((escvals.real - xmin) * (width - 2) / (xmax - xmin))
    ypoints = np.uint64((escvals.imag - ymin) * (height - 2) / (ymax - ymin))
    points = ((ypoints.clip(0, height-1), xpoints.clip(0, width-1)))
    surface[points] += 1

print("Stage2 complete:", time.time() - starttime, "seconds")
img = Image.fromarray(np.uint8(surface.clip(0,255)).T, mode="L")
img.save("buddha.png")
img.show()