diff --git a/pyscriptjs/examples/index.html b/pyscriptjs/examples/index.html
index 02d1eaa..6710acf 100644
--- a/pyscriptjs/examples/index.html
+++ b/pyscriptjs/examples/index.html
@@ -45,6 +45,8 @@
Minimal d3 demo demonstrating how to create a visualization
+
+ Visualization of Mandelbrot and Julia sets with NumPy and HTML5 canvas
A Python REPL (Read Eval Print Loop).
diff --git a/pyscriptjs/examples/numpy_canvas_fractals.html b/pyscriptjs/examples/numpy_canvas_fractals.html
new file mode 100644
index 0000000..c48653f
--- /dev/null
+++ b/pyscriptjs/examples/numpy_canvas_fractals.html
@@ -0,0 +1,236 @@
+
+
+ Visualization of Mandelbrot and Julia sets with NumPy and HTML5 canvas
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+- numpy
+
+
+
+from pyodide import to_js
+
+import numpy as np
+
+from js import (
+ console,
+ document,
+ devicePixelRatio,
+ ImageData,
+ Uint8ClampedArray,
+ CanvasRenderingContext2D as Context2d,
+)
+
+Magma256 = np.array([
+ [0x00, 0x00, 0x03], [0x00, 0x00, 0x04], [0x00, 0x00, 0x06], [0x01, 0x00, 0x07], [0x01, 0x01, 0x09], [0x01, 0x01, 0x0b],
+ [0x02, 0x02, 0x0d], [0x02, 0x02, 0x0f], [0x03, 0x03, 0x11], [0x04, 0x03, 0x13], [0x04, 0x04, 0x15], [0x05, 0x04, 0x17],
+ [0x06, 0x05, 0x19], [0x07, 0x05, 0x1b], [0x08, 0x06, 0x1d], [0x09, 0x07, 0x1f], [0x0a, 0x07, 0x22], [0x0b, 0x08, 0x24],
+ [0x0c, 0x09, 0x26], [0x0d, 0x0a, 0x28], [0x0e, 0x0a, 0x2a], [0x0f, 0x0b, 0x2c], [0x10, 0x0c, 0x2f], [0x11, 0x0c, 0x31],
+ [0x12, 0x0d, 0x33], [0x14, 0x0d, 0x35], [0x15, 0x0e, 0x38], [0x16, 0x0e, 0x3a], [0x17, 0x0f, 0x3c], [0x18, 0x0f, 0x3f],
+ [0x1a, 0x10, 0x41], [0x1b, 0x10, 0x44], [0x1c, 0x10, 0x46], [0x1e, 0x10, 0x49], [0x1f, 0x11, 0x4b], [0x20, 0x11, 0x4d],
+ [0x22, 0x11, 0x50], [0x23, 0x11, 0x52], [0x25, 0x11, 0x55], [0x26, 0x11, 0x57], [0x28, 0x11, 0x59], [0x2a, 0x11, 0x5c],
+ [0x2b, 0x11, 0x5e], [0x2d, 0x10, 0x60], [0x2f, 0x10, 0x62], [0x30, 0x10, 0x65], [0x32, 0x10, 0x67], [0x34, 0x10, 0x68],
+ [0x35, 0x0f, 0x6a], [0x37, 0x0f, 0x6c], [0x39, 0x0f, 0x6e], [0x3b, 0x0f, 0x6f], [0x3c, 0x0f, 0x71], [0x3e, 0x0f, 0x72],
+ [0x40, 0x0f, 0x73], [0x42, 0x0f, 0x74], [0x43, 0x0f, 0x75], [0x45, 0x0f, 0x76], [0x47, 0x0f, 0x77], [0x48, 0x10, 0x78],
+ [0x4a, 0x10, 0x79], [0x4b, 0x10, 0x79], [0x4d, 0x11, 0x7a], [0x4f, 0x11, 0x7b], [0x50, 0x12, 0x7b], [0x52, 0x12, 0x7c],
+ [0x53, 0x13, 0x7c], [0x55, 0x13, 0x7d], [0x57, 0x14, 0x7d], [0x58, 0x15, 0x7e], [0x5a, 0x15, 0x7e], [0x5b, 0x16, 0x7e],
+ [0x5d, 0x17, 0x7e], [0x5e, 0x17, 0x7f], [0x60, 0x18, 0x7f], [0x61, 0x18, 0x7f], [0x63, 0x19, 0x7f], [0x65, 0x1a, 0x80],
+ [0x66, 0x1a, 0x80], [0x68, 0x1b, 0x80], [0x69, 0x1c, 0x80], [0x6b, 0x1c, 0x80], [0x6c, 0x1d, 0x80], [0x6e, 0x1e, 0x81],
+ [0x6f, 0x1e, 0x81], [0x71, 0x1f, 0x81], [0x73, 0x1f, 0x81], [0x74, 0x20, 0x81], [0x76, 0x21, 0x81], [0x77, 0x21, 0x81],
+ [0x79, 0x22, 0x81], [0x7a, 0x22, 0x81], [0x7c, 0x23, 0x81], [0x7e, 0x24, 0x81], [0x7f, 0x24, 0x81], [0x81, 0x25, 0x81],
+ [0x82, 0x25, 0x81], [0x84, 0x26, 0x81], [0x85, 0x26, 0x81], [0x87, 0x27, 0x81], [0x89, 0x28, 0x81], [0x8a, 0x28, 0x81],
+ [0x8c, 0x29, 0x80], [0x8d, 0x29, 0x80], [0x8f, 0x2a, 0x80], [0x91, 0x2a, 0x80], [0x92, 0x2b, 0x80], [0x94, 0x2b, 0x80],
+ [0x95, 0x2c, 0x80], [0x97, 0x2c, 0x7f], [0x99, 0x2d, 0x7f], [0x9a, 0x2d, 0x7f], [0x9c, 0x2e, 0x7f], [0x9e, 0x2e, 0x7e],
+ [0x9f, 0x2f, 0x7e], [0xa1, 0x2f, 0x7e], [0xa3, 0x30, 0x7e], [0xa4, 0x30, 0x7d], [0xa6, 0x31, 0x7d], [0xa7, 0x31, 0x7d],
+ [0xa9, 0x32, 0x7c], [0xab, 0x33, 0x7c], [0xac, 0x33, 0x7b], [0xae, 0x34, 0x7b], [0xb0, 0x34, 0x7b], [0xb1, 0x35, 0x7a],
+ [0xb3, 0x35, 0x7a], [0xb5, 0x36, 0x79], [0xb6, 0x36, 0x79], [0xb8, 0x37, 0x78], [0xb9, 0x37, 0x78], [0xbb, 0x38, 0x77],
+ [0xbd, 0x39, 0x77], [0xbe, 0x39, 0x76], [0xc0, 0x3a, 0x75], [0xc2, 0x3a, 0x75], [0xc3, 0x3b, 0x74], [0xc5, 0x3c, 0x74],
+ [0xc6, 0x3c, 0x73], [0xc8, 0x3d, 0x72], [0xca, 0x3e, 0x72], [0xcb, 0x3e, 0x71], [0xcd, 0x3f, 0x70], [0xce, 0x40, 0x70],
+ [0xd0, 0x41, 0x6f], [0xd1, 0x42, 0x6e], [0xd3, 0x42, 0x6d], [0xd4, 0x43, 0x6d], [0xd6, 0x44, 0x6c], [0xd7, 0x45, 0x6b],
+ [0xd9, 0x46, 0x6a], [0xda, 0x47, 0x69], [0xdc, 0x48, 0x69], [0xdd, 0x49, 0x68], [0xde, 0x4a, 0x67], [0xe0, 0x4b, 0x66],
+ [0xe1, 0x4c, 0x66], [0xe2, 0x4d, 0x65], [0xe4, 0x4e, 0x64], [0xe5, 0x50, 0x63], [0xe6, 0x51, 0x62], [0xe7, 0x52, 0x62],
+ [0xe8, 0x54, 0x61], [0xea, 0x55, 0x60], [0xeb, 0x56, 0x60], [0xec, 0x58, 0x5f], [0xed, 0x59, 0x5f], [0xee, 0x5b, 0x5e],
+ [0xee, 0x5d, 0x5d], [0xef, 0x5e, 0x5d], [0xf0, 0x60, 0x5d], [0xf1, 0x61, 0x5c], [0xf2, 0x63, 0x5c], [0xf3, 0x65, 0x5c],
+ [0xf3, 0x67, 0x5b], [0xf4, 0x68, 0x5b], [0xf5, 0x6a, 0x5b], [0xf5, 0x6c, 0x5b], [0xf6, 0x6e, 0x5b], [0xf6, 0x70, 0x5b],
+ [0xf7, 0x71, 0x5b], [0xf7, 0x73, 0x5c], [0xf8, 0x75, 0x5c], [0xf8, 0x77, 0x5c], [0xf9, 0x79, 0x5c], [0xf9, 0x7b, 0x5d],
+ [0xf9, 0x7d, 0x5d], [0xfa, 0x7f, 0x5e], [0xfa, 0x80, 0x5e], [0xfa, 0x82, 0x5f], [0xfb, 0x84, 0x60], [0xfb, 0x86, 0x60],
+ [0xfb, 0x88, 0x61], [0xfb, 0x8a, 0x62], [0xfc, 0x8c, 0x63], [0xfc, 0x8e, 0x63], [0xfc, 0x90, 0x64], [0xfc, 0x92, 0x65],
+ [0xfc, 0x93, 0x66], [0xfd, 0x95, 0x67], [0xfd, 0x97, 0x68], [0xfd, 0x99, 0x69], [0xfd, 0x9b, 0x6a], [0xfd, 0x9d, 0x6b],
+ [0xfd, 0x9f, 0x6c], [0xfd, 0xa1, 0x6e], [0xfd, 0xa2, 0x6f], [0xfd, 0xa4, 0x70], [0xfe, 0xa6, 0x71], [0xfe, 0xa8, 0x73],
+ [0xfe, 0xaa, 0x74], [0xfe, 0xac, 0x75], [0xfe, 0xae, 0x76], [0xfe, 0xaf, 0x78], [0xfe, 0xb1, 0x79], [0xfe, 0xb3, 0x7b],
+ [0xfe, 0xb5, 0x7c], [0xfe, 0xb7, 0x7d], [0xfe, 0xb9, 0x7f], [0xfe, 0xbb, 0x80], [0xfe, 0xbc, 0x82], [0xfe, 0xbe, 0x83],
+ [0xfe, 0xc0, 0x85], [0xfe, 0xc2, 0x86], [0xfe, 0xc4, 0x88], [0xfe, 0xc6, 0x89], [0xfe, 0xc7, 0x8b], [0xfe, 0xc9, 0x8d],
+ [0xfe, 0xcb, 0x8e], [0xfd, 0xcd, 0x90], [0xfd, 0xcf, 0x92], [0xfd, 0xd1, 0x93], [0xfd, 0xd2, 0x95], [0xfd, 0xd4, 0x97],
+ [0xfd, 0xd6, 0x98], [0xfd, 0xd8, 0x9a], [0xfd, 0xda, 0x9c], [0xfd, 0xdc, 0x9d], [0xfd, 0xdd, 0x9f], [0xfd, 0xdf, 0xa1],
+ [0xfd, 0xe1, 0xa3], [0xfc, 0xe3, 0xa5], [0xfc, 0xe5, 0xa6], [0xfc, 0xe6, 0xa8], [0xfc, 0xe8, 0xaa], [0xfc, 0xea, 0xac],
+ [0xfc, 0xec, 0xae], [0xfc, 0xee, 0xb0], [0xfc, 0xf0, 0xb1], [0xfc, 0xf1, 0xb3], [0xfc, 0xf3, 0xb5], [0xfc, 0xf5, 0xb7],
+ [0xfb, 0xf7, 0xb9], [0xfb, 0xf9, 0xbb], [0xfb, 0xfa, 0xbd], [0xfb, 0xfc, 0xbf],
+], dtype="uint8")
+
+def mandelbrot(width: int, height: int, *,
+ x: float = -0.5, y: float = 0, zoom: int = 1, max_iterations: int = 100) -> np.array:
+ """
+ From https://www.learnpythonwithrune.org/numpy-compute-mandelbrot-set-by-vectorization/.
+ """
+ # To make navigation easier we calculate these values
+ x_width, y_height = 1.5, 1.5*height/width
+ x_from, x_to = x - x_width/zoom, x + x_width/zoom
+ y_from, y_to = y - y_height/zoom, y + y_height/zoom
+
+ # Here the actual algorithm starts
+ x = np.linspace(x_from, x_to, width).reshape((1, width))
+ y = np.linspace(y_from, y_to, height).reshape((height, 1))
+ c = x + 1j*y
+
+ # Initialize z to all zero
+ z = np.zeros(c.shape, dtype=np.complex128)
+
+ # To keep track in which iteration the point diverged
+ div_time = np.zeros(z.shape, dtype=int)
+
+ # To keep track on which points did not converge so far
+ m = np.full(c.shape, True, dtype=bool)
+ for i in range(max_iterations):
+ z[m] = z[m]**2 + c[m]
+ diverged = np.greater(np.abs(z), 2, out=np.full(c.shape, False), where=m) # Find diverging
+ div_time[diverged] = i # set the value of the diverged iteration number
+ m[np.abs(z) > 2] = False # to remember which have diverged
+
+ return div_time
+
+def julia(width: int, height: int, *,
+ c: complex = -0.4 + 0.6j, x: float = 0, y: float = 0, zoom: int = 1, max_iterations: int = 100) -> np.array:
+ """
+ From https://www.learnpythonwithrune.org/numpy-calculate-the-julia-set-with-vectorization/.
+ """
+ # To make navigation easier we calculate these values
+ x_width, y_height = 1.5, 1.5*height/width
+ x_from, x_to = x - x_width/zoom, x + x_width/zoom
+ y_from, y_to = y - y_height/zoom, y + y_height/zoom
+
+ # Here the actual algorithm starts
+ x = np.linspace(x_from, x_to, width).reshape((1, width))
+ y = np.linspace(y_from, y_to, height).reshape((height, 1))
+ z = x + 1j*y
+
+ # Initialize z to all zero
+ c = np.full(z.shape, c)
+
+ # To keep track in which iteration the point diverged
+ div_time = np.zeros(z.shape, dtype=int)
+
+ # To keep track on which points did not converge so far
+ m = np.full(c.shape, True, dtype=bool)
+ for i in range(max_iterations):
+ z[m] = z[m]**2 + c[m]
+ m[np.abs(z) > 2] = False
+ div_time[m] = i
+
+ return div_time
+
+def create_canvas(width: int, height: int, target: str) -> Context2d:
+ pixel_ratio = devicePixelRatio
+
+ canvas = document.createElement("canvas")
+ ctx = canvas.getContext("2d")
+
+ canvas.style.width = f"{width}px"
+ canvas.style.height = f"{height}px"
+
+ canvas.width = width*pixel_ratio
+ canvas.height = height*pixel_ratio
+
+ ctx.scale(pixel_ratio, pixel_ratio)
+ ctx.translate(0.5, 0.5)
+
+ ctx.clearRect(0, 0, width, height)
+
+ el = document.querySelector(target)
+ el.replaceChildren(canvas)
+
+ return ctx
+
+def color_map(array: np.array, palette: np.array) -> np.array:
+ size, _ = palette.shape
+ index = (array/array.max()*(size - 1)).round().astype("uint8")
+
+ width, height = array.shape
+ image = np.full((width, height, 4), 0xff, dtype=np.uint8)
+ image[:, :, :3] = palette[index]
+
+ return image
+
+def draw_image(ctx: Context2d, image: np.array) -> None:
+ data = Uint8ClampedArray.new(to_js(image.tobytes()))
+ width, height, _ = image.shape
+ image_data = ImageData.new(data, width, height)
+ ctx.putImageData(image_data, 0, 0)
+
+def draw_mandelbrot(width: int, height: int) -> None:
+ ctx = create_canvas(width, height, "#mandelbrot")
+
+ console.log("Computing Mandelbrot set ...")
+ console.time("mandelbrot")
+ array = mandelbrot(width, height)
+ console.timeEnd("mandelbrot")
+
+ image = color_map(array, Magma256)
+ draw_image(ctx, image)
+
+def draw_julia(width: int, height: int) -> None:
+ ctx = create_canvas(width, height, "#julia")
+
+ console.log("Computing Julia set ...")
+ console.time("julia")
+ array = julia(width, height)
+ console.timeEnd("julia")
+
+ image = color_map(array, Magma256)
+ draw_image(ctx, image)
+
+draw_mandelbrot(600, 600)
+draw_julia(600, 600)
+
+
+
+