Tile Blur (WebGPU separable Gaussian tile-shared-memory blur)

Math essence

For each axis (x, y) apply 1D Gaussian convolution with kernel k_i = exp(-(i-µ)²/(2σ²)) / Σ k_j. Padded to keep spatial extent unchanged. Two passes replace O(K²) cost with O(2K). Step adds Gaussian perturbation plus two 'ink blot' disk features and oscillates sigma ±1.5 around the base.

Parameters

Reference

First 60 lines of reference.py. Full source: library_external/anadol_webgpu/tileBlur/reference.py

"""tileBlur — WebGPU separable tile-shared-memory Gaussian blur.

The classical GPU image-blur pattern:

  1. Horizontal pass: each workgroup (TILE x TILE) loads a tile of
     the input image plus a halo of border pixels into
     workgroup-shared memory, workgroupBarrier(), convolves each pixel
     with a 1D Gaussian kernel along the x-axis, writes to a ping-pong
     texture.
  2. Vertical pass: same algorithm along the y-axis reading the
     intermediate texture as input.

The 1D factorisation reduces an O(K²) kernel to O(2K) texture
accesses per pixel, which is exactly why "separable Gaussian blur"
is the canonical GPU tutorial.

Here we implement the algorithm in numpy using two passes of 1D
Gaussian convolution along rows then columns; the input image is a
synthesised textured pattern so the blur result is interesting.
Each step "tick" adds a small noise perturbation to the source so
the blur evolves over time.

Reference:
- https://www.mysimulator.uk/content/articles/webgpu-compute.html
  ("2D convolution using workgroup tiles ... separable Gaussian kernel
  in shared memory")
- https://lisyarus.github.io/blog/posts/compute-blur.html
  ("Compute shaders in graphics: Gaussian blur")
"""
from __future__ import annotations

from typing import ClassVar

import numpy as np


def _gauss_kernel_1d(sigma: float, radius: int | None = None) -> np.ndarray:
    if radius is None:
        radius = max(1, int(np.ceil(3.0 * sigma)))
    x = np.arange(-radius, radius + 1, dtype=np.float64)
    k = np.exp(-(x * x) / (2.0 * sigma * sigma))
    k /= k.sum()
    return k.astype(np.float64)


def _separable_gaussian_blur(img: np.ndarray, sigma_x: float, sigma_y: float) -> np.ndarray:
    """Apply a separable 2D Gaussian blur: 1D along rows then along cols."""
    H, W, C = img.shape
    out = img.astype(np.float32)

    if sigma_x > 0:
        kx = _gauss_kernel_1d(sigma_x)
        r = kx.size // 2
        # pad each row
        pad = np.zeros((H, r, C), dtype=np.float32)
        padded = np.concatenate([pad, out, pad], axis=1)
        tmp = np.zeros_like(out)
        for i in range(kx.size):
            tmp += kx[i] * padded[:, i:i + W, :]
        out = tmp