Circular blur kernel optimizer

Disc kernels as sums of phased Gaussians (cos(c r²)·rw + sin(c r²)·iw)·e−b r² and, optionally, components with a repeated rate (multiplicity m up to 10) and extra r², r⁴, … terms. For rendering, the components are transformed together into separable channels, each a horizontal and a vertical complex 1-D pass, by a fixed linear map (reported with the kernel) chosen so that the passes stay about as small as the kernel and 16-bit intermediates suffice. Global search by differential evolution over decays and frequencies, with exact minimax weights solved by Remez exchange on CPU workers (or over all parameters on the GPU with jax-js-nonconsuming), then local Remez refinement or gradient fine-tuning; the blur demo is plain WebGPU shaders. After Olli Niemitalo: Circularly symmetric convolution and lens blur.

Blur demo WebGPU, using the current kernel

Shader code (WGSL)

The compute shader follows the intermediate storage and radius-interpolation settings. Its header describes the buffers, tap tables built from the separable channels under "Current kernel", and dispatch order.

Scene
Lens
Varying blur

How the varying blur works

Display
Rendering
Capture and bug reports

Kernel optimizer fits the disc

Starting…

Kernel
About multiplicities

Per component, 1 to 10: m coinciding rates in closed form, (a0 + a1 r² + … + am−1 r2(m−1))·eλr², costing m complex channels per pass (1 = single, 2 = merged pair). Example: 2, 2, 1. More in the README.

Target and error
About negative gain

Errors where the kernel is negative count "negative gain" times. Above 1, the optimizers trade ripple for fewer negative values, which look like cheap sharpening (as in the article's non-negative sets).

Differential evolution
Initialize population
About search ranges

Remez exchange refinement
How it works

Minimax on the dense grid (Dense grid samples per band, below), in double precision on the CPU. First solves the optimal weights for the current decays and frequencies exactly by an exchange (dual simplex) algorithm. Then each iteration tries a trust-region Gauss–Newton step whose linearized minimax problem is solved by exchange, and Newton steps on the equal-ripple (Remez) equations of the extremal points, keeping whichever lowers the ripple. Stops when no step improves. Usually converges in a few iterations from a kernel near a local optimum.

Gradient fine-tuning
How it works

Minimizes a log-sum-exp smooth maximum of |error| on the dense grid with Adam in parameter-relative units (learning rate defaults to 0.001/3N−2 and halves every third of the run; sharpness, relative to the current ripple, rises 10×). Keeps the best kernel by true maximum error.

2-D kernel normalized to its peak: positive values in gray levels, negative values in red, on the same brightness scale.

Current kernel

How the kernel is formed

Each component is separable on its own, but optimized components often have large weights, and their imaginary parts cancel: the kernel is only the real part of the sum. Passes that large would lose the kernel in 16-bit intermediates. So the blur instead uses channels from a joint eigendecomposition of the real kernel, whose passes are no larger than the kernel needs, at the same channel count. Each image below is normalized to its own peak (shown below it), in the colors of the 2-D kernel above, over the support |x|, |y| ≤ 1 + transition.

Components

Each component's own separable form: one complex channel √a0 eλx² × √a0 eλy² for a single component (the article's), or m channels x2j eλx² × Qj(y²) eλy² for multiplicity m, summed into the component u, and the components into the kernel K. Rows are real terms h(x) v(y) as in Channels; their sizes show any cancellation between components. The blur doesn't compute these: this is the math view, showing that each component (the sum of its terms) is circularly symmetric.

Channels

The blur's separable passes, transformed from the components as a whole. Each row is a horizontal 1-D kernel h(x) times a vertical one v(y); pairs of these real terms are the complex channels Re[h(x) v(y)] = hrevre − himvim, summed into K as the shader accumulates them. The shader never forms the two terms of a channel separately: its horizontal pass filters with complex h, and its vertical pass multiplies by complex v and keeps the real part, which is their sum. A channel alone is not circularly symmetric; only K is. Compare the peaks with those of the components' own forms.