# The filters a dither is seen through, as tap lists [(dx, dy, weight)] with weights summing to 1.
import numpy as np

def kuwahara(R=5, stride=1):
    """Starlock's painterly pass (public/ideas/starlock/reveal3d.js), linearised: on a flat, dithered
    input every sector's variance is at most a quarter of a level squared, so its weight
    1 / (1 + (1000 var)^4) is within 2e-8 of 1 and the generalised Kuwahara is, to that precision, the
    plain mean of its 8 Gaussian-weighted sector means (the centre tap counts in all 8; the fresh check
    measured the difference at 2e-9 levels). Taps `stride` pixels apart."""
    r2 = R * R + .5; taps = []; w = np.zeros(8)
    for j in range(-R, R + 1):
        for i in range(-R, R + 1):
            d = i * i + j * j
            if d > r2: continue
            g = np.exp(-d / r2 * 2)
            if d < .5: taps.append((i, j, g, None)); w += g; continue
            k = int(np.floor((np.arctan2(j, i) + 3.14159265) / .78539816 + 1e-5) % 8)  # (the GPU's edge rule, see the 10-22 entry)
            w[k] += g; taps.append((i, j, g, k))
    out = []
    for i, j, g, k in taps:
        wt = sum(g / w[kk] for kk in range(8)) / 8 if k is None else g / w[k] / 8
        out.append((i * stride, j * stride, wt))
    return out

def box(n, stride=1):
    return [(i * stride, j * stride, 1 / n / n) for j in range(n) for i in range(n)]

def gaussian(sigma, stride=1):
    R = int(np.ceil(3 * sigma)); t = [(i, j, np.exp(-(i * i + j * j) / (2 * sigma * sigma))) for j in range(-R, R + 1) for i in range(-R, R + 1)]
    s = sum(w for _, _, w in t); return [(i * stride, j * stride, w / s) for i, j, w in t]

def lee_disc(r=2):
    """Equal weights on the lattice points with |i| + |j| <= r (the 13-point 'Lee' ball for r = 2)."""
    t = [(i, j) for j in range(-r, r + 1) for i in range(-r, r + 1) if abs(i) + abs(j) <= r]
    return [(i, j, 1 / len(t)) for i, j in t]

def response(taps, fx, fy):
    """H(xi) = sum_q h_q exp(2 pi i xi.q), evaluated exactly at frequencies (fx, fy) in cycles/pixel."""
    fx = np.asarray(fx, float); fy = np.asarray(fy, float); H = np.zeros(np.broadcast(fx, fy).shape, complex)
    for dx, dy, w in taps: H += w * np.exp(2j * np.pi * (fx * dx + fy * dy))
    return H

def grid_response(taps, N):
    """|H|^2 on the N x N grid of frequencies (j/N, i/N): row = y frequency, column = x frequency."""
    K = np.zeros((N, N))
    for dx, dy, w in taps: K[dy % N, dx % N] += w
    return np.abs(np.fft.fft2(K)) ** 2

def centre_plus_gaussian(c, sigma, stride=1):
    """A share c of the weight on the pixel itself, the rest a Gaussian blur: the shape of any filter that
    keeps the pixel (Kuwahara, bilateral, a TAA blend) in its simplest form."""
    g = gaussian(sigma, stride)
    return [(dx, dy, (1 - c) * w + (c if dx == 0 and dy == 0 else 0)) for dx, dy, w in g]