How to draw Buddhabrot fractal?

Question

I'm trying to implement Buddhabrot fractal in Python. I read a lot of articles and posts but I think that I missunderstood something (just see the image). Someone can write a pseudocode?

My code is this:

from multiprocessing import Pool
from random import randrange
import matplotlib.pyplot as plt
import numpy as np
from math import ceil


maxiter       = 1000
points        = 1000
xmin, xmax    = -2, 1
ymin, ymax    = -2, 1
cores         = 4
width, height = 200, 200


maxn         = width * height
incrx, incry = abs(xmax - xmin) / width, abs(ymax - ymin) / height

def randomComplexGenerator():
  for i in range(points):
    n = randrange(maxn)
    yield complex(n // height * incrx, n % width * incry)

def buddhabrot(c):
  m, z, i = np.zeros((width, height)), c, 0
  while i < maxiter and abs(z) < 2:
    x, y     = ceil(z.real / incrx), ceil(z.imag / incry)
    m[x, y] += 1
    z        = z ** 2 + c
    i       += 1
  return m if i == maxiter else 0

if __name__ == '__main__':
  a = np.linspace(xmin, xmax, width)
  b = np.linspace(ymin, ymax, height)
  with Pool(cores) as p:
    ms = p.map(buddhabrot, (c for c in randomComplexGenerator()))
  res = 0
  for m in ms:
    res += m
  plt.axis('off')
  plt.imshow(res)
  plt.show()

The image generated with my code is this (lel):

Answer 1

After days, this is the code that I created, which seems to generate appropriately the fractal. Any performance suggestion is welcome.

from multiprocessing import Pool
from random import randrange
import matplotlib.pyplot as plt
import numpy as np


cores         = 4
maxiter       = 10000
points        = 1000000
width, height = 200, 200
rdom, idom    = (-2, 2), (-2, 2)

xdom, ydom    = (0, width - 1), (0, height - 1)

def randomComplex():
  r = np.interp(randrange(xdom[0], xdom[1]), xdom, rdom)
  i = np.interp(randrange(ydom[0], ydom[1]), ydom, idom)
  return (r, i)

def complex2pixel(c):
  x = int(np.interp(c[0], rdom, xdom))
  y = int(np.interp(c[1], idom, ydom))
  return (x, y)

def escapedPixels(c):
  pixels, z = {}, c
  for i in range(maxiter):
    z2 = (z[0] * z[0], z[1] * z[1]) 
    if z2[0] + z2[1] > 4: break
    p = complex2pixel(z)
    try: pixels[p] += 1
    except: pixels[p] = 1
    z = (z2[0] - z2[1] + c[0], 2 * z[0] * z[1] + c[1])
  return pixels if i < maxiter - 1 else {}

if __name__ == '__main__':
  with Pool(cores) as p:
    ds = p.map(escapedPixels, (randomComplex() for i in range(points)))
  m = np.zeros((width, height))
  for d in ds:
    for p in d:
      m[p] += d[p]
  plt.axis('off')
  plt.imshow(m)
  plt.show()