From d939e5ad15ebc780c38b6d6b6588ab213dd08484 Mon Sep 17 00:00:00 2001 From: Mateusz Paprocki Date: Fri, 29 Apr 2022 23:36:44 +0200 Subject: [PATCH] Improvements to numpy_canvas_fractals example --- pyscriptjs/examples/fractals.py | 47 ++++ .../examples/numpy_canvas_fractals.html | 243 +++++++++++++++--- 2 files changed, 260 insertions(+), 30 deletions(-) diff --git a/pyscriptjs/examples/fractals.py b/pyscriptjs/examples/fractals.py index 8d4ad78e..65e0f31a 100644 --- a/pyscriptjs/examples/fractals.py +++ b/pyscriptjs/examples/fractals.py @@ -1,4 +1,6 @@ +from typing import Tuple import numpy as np +from numpy.polynomial import Polynomial def mandelbrot(width: int, height: int, *, x: float = -0.5, y: float = 0, zoom: int = 1, max_iterations: int = 100) -> np.array: @@ -60,3 +62,48 @@ def julia(width: int, height: int, *, div_time[m] = i return div_time + +Range = Tuple[float, float] + +def newton(width: int, height: int, *, + p: Polynomial, a: complex, xr: Range = (-2.5, 1), yr: Range = (-1, 1), max_iterations: int = 100) -> (np.array, np.array): + """ """ + # To make navigation easier we calculate these values + x_from, x_to = xr + y_from, y_to = yr + + # 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 + + # Compute the derivative + dp = p.deriv() + + # Compute roots + roots = p.roots() + epsilon = 1e-5 + + # Set the initial conditions + a = np.full(z.shape, a) + + # 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(a.shape, True, dtype=bool) + + # To keep track which root each point converged to + r = np.full(a.shape, 0, dtype=int) + + for i in range(max_iterations): + z[m] = z[m] - a[m]*p(z[m])/dp(z[m]) + + for j, root in enumerate(roots): + converged = (np.abs(z.real - root.real) < epsilon) & (np.abs(z.imag - root.imag) < epsilon) + m[converged] = False + r[converged] = j + 1 + + div_time[m] = i + + return div_time, r diff --git a/pyscriptjs/examples/numpy_canvas_fractals.html b/pyscriptjs/examples/numpy_canvas_fractals.html index cbee92fb..2247747d 100644 --- a/pyscriptjs/examples/numpy_canvas_fractals.html +++ b/pyscriptjs/examples/numpy_canvas_fractals.html @@ -1,6 +1,6 @@ - Visualization of Mandelbrot and Julia sets with NumPy and HTML5 canvas + Visualization of Mandelbrot, Julia and Newton sets with NumPy and HTML5 canvas @@ -17,6 +17,10 @@ animation: spin 1s ease-in-out infinite; } + canvas { + display: none; + } + @keyframes spin { to { transform: rotate(360deg); @@ -28,35 +32,68 @@ -
-
+
+
Mandelbrot set
-
+
+
-
+
Julia set
-
+
+ +
+
+
+
Newton set
+
+
p(z) =
+
a =
+
+ x = [ + + , + + ] +
+
+ y = [ + + , + + ] +
+
+
convergence
+
iterations
+
+
+
+
+
- numpy +- sympy - paths: - /palettes.py - /fractals.py -from pyodide import to_js +from pyodide import to_js, create_proxy import numpy as np +import sympy from palettes import Magma256 -from fractals import mandelbrot, julia +from fractals import mandelbrot, julia, newton from js import ( console, @@ -65,28 +102,20 @@ from js import ( ImageData, Uint8ClampedArray, CanvasRenderingContext2D as Context2d, + requestAnimationFrame, ) -def create_canvas(width: int, height: int, target: str) -> Context2d: - pixel_ratio = devicePixelRatio - - canvas = document.createElement("canvas") +def prepare_canvas(width: int, height: int, canvas: Element) -> Context2d: 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) + canvas.width = width + canvas.height = height 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: @@ -105,30 +134,184 @@ def draw_image(ctx: Context2d, image: np.array) -> None: 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") +width, height = 600, 600 + +async def draw_mandelbrot() -> None: + spinner = document.querySelector("#mandelbrot .loading") + canvas = document.querySelector("#mandelbrot canvas") + + spinner.style.display = "" + canvas.style.display = "none" + + ctx = prepare_canvas(width, height, canvas) console.log("Computing Mandelbrot set ...") console.time("mandelbrot") - array = mandelbrot(width, height) + iters = mandelbrot(width, height) console.timeEnd("mandelbrot") - image = color_map(array, Magma256) + image = color_map(iters, Magma256) draw_image(ctx, image) -def draw_julia(width: int, height: int) -> None: - ctx = create_canvas(width, height, "#julia") + spinner.style.display = "none" + canvas.style.display = "block" + +async def draw_julia() -> None: + spinner = document.querySelector("#julia .loading") + canvas = document.querySelector("#julia canvas") + + spinner.style.display = "" + canvas.style.display = "none" + + ctx = prepare_canvas(width, height, canvas) console.log("Computing Julia set ...") console.time("julia") - array = julia(width, height) + iters = julia(width, height) console.timeEnd("julia") - image = color_map(array, Magma256) + image = color_map(iters, Magma256) draw_image(ctx, image) -draw_mandelbrot(600, 600) -draw_julia(600, 600) + spinner.style.display = "none" + canvas.style.display = "block" + +def ranges(): + x0_in = document.querySelector("#x0") + x1_in = document.querySelector("#x1") + y0_in = document.querySelector("#y0") + y1_in = document.querySelector("#y1") + + xr = (float(x0_in.value), float(x1_in.value)) + yr = (float(y0_in.value), float(y1_in.value)) + + return xr, yr + +current_image = None +async def draw_newton() -> None: + spinner = document.querySelector("#newton .loading") + canvas = document.querySelector("#newton canvas") + + spinner.style.display = "" + canvas.style.display = "none" + + ctx = prepare_canvas(width, height, canvas) + + console.log("Computing Newton set ...") + + poly_in = document.querySelector("#poly") + coef_in = document.querySelector("#coef") + conv_in = document.querySelector("#conv") + iter_in = document.querySelector("#iter") + + xr, yr = ranges() + + # z**3 - 1 + # z**8 + 15*z**4 - 16 + # z**3 - 2*z + 2 + + expr = sympy.parse_expr(poly_in.value) + coeffs = [ complex(c) for c in reversed(sympy.Poly(expr, sympy.Symbol("z")).all_coeffs()) ] + poly = np.polynomial.Polynomial(coeffs) + + coef = complex(sympy.parse_expr(coef_in.value)) + + console.time("newton") + iters, roots = newton(width, height, p=poly, a=coef, xr=xr, yr=yr) + console.timeEnd("newton") + + if conv_in.checked: + n = poly.degree() + 1 + k = int(len(Magma256)/n) + + colors = Magma256[::k, :][:n] + colors[0, :] = [255, 0, 0] # red: no convergence + + image = color_map(roots, colors) + else: + image = color_map(iters, Magma256) + + global current_image + current_image = image + draw_image(ctx, image) + + spinner.style.display = "none" + canvas.style.display = "block" + +handler = create_proxy(lambda _event: draw_newton()) +document.querySelector("#newton fieldset").addEventListener("change", handler) + +canvas = document.querySelector("#newton canvas") + +is_selecting = False +init_sx, init_sy = None, None +sx, sy = None, None +async def mousemove(event): + global is_selecting + global init_sx + global init_sy + global sx + global sy + + def invert(sx, source_range, target_range): + source_start, source_end = source_range + target_start, target_end = target_range + factor = (target_end - target_start)/(source_end - source_start) + offset = -(factor * source_start) + target_start + return (sx - offset) / factor + + bds = canvas.getBoundingClientRect() + event_sx, event_sy = event.clientX - bds.x, event.clientY - bds.y + + ctx = canvas.getContext("2d") + + pressed = event.buttons == 1 + if is_selecting: + if not pressed: + xr, yr = ranges() + + x0 = invert(init_sx, xr, (0, width)) + x1 = invert(sx, xr, (0, width)) + y0 = invert(init_sy, yr, (0, height)) + y1 = invert(sy, yr, (0, height)) + + document.querySelector("#x0").value = x0 + document.querySelector("#x1").value = x1 + document.querySelector("#y0").value = y0 + document.querySelector("#y1").value = y1 + + is_selecting = False + init_sx, init_sy = None, None + sx, sy = init_sx, init_sy + + await draw_newton() + else: + ctx.save() + ctx.clearRect(0, 0, width, height) + draw_image(ctx, current_image) + sx, sy = event_sx, event_sy + ctx.beginPath() + ctx.rect(init_sx, init_sy, sx - init_sx, sy - init_sy) + ctx.fillStyle = "rgba(255, 255, 255, 0.4)" + ctx.strokeStyle = "rgba(255, 255, 255, 1.0)" + ctx.fill() + ctx.stroke() + ctx.restore() + else: + if pressed: + is_selecting = True + init_sx, init_sy = event_sx, event_sy + sx, sy = init_sx, init_sy + +canvas.addEventListener("mousemove", create_proxy(mousemove)) + +import asyncio + +_ = await asyncio.gather( + draw_mandelbrot(), + draw_julia(), + draw_newton(), +)