# Independent linearisation of Starlock's Kuwahara (reveal3d.js): on flat dithered input every sector
# weight a ~= 1, so out = (1/8) sum_k m_k, m_k = sum_{q in k} g_q c_q / sum_{q in k} g_q (centre in all 8).
# Returns a dict {(dx,dy): weight}. tie: how atan-boundary points are assigned ('up' = +eps, 'down' = -eps).
import math
import numpy as np
def kuw(R=5, tie='up'):
    r2 = R*R + .5; eps = 1e-7 if tie == 'up' else -1e-7
    sect = {}; centre = None
    for j in range(-R, R+1):
        for i in range(-R, R+1):
            d = i*i + j*j
            if d > r2: continue
            g = math.exp(-d/r2*2)
            if d < .5: centre = g; continue
            k = int(math.floor((math.atan2(j, i) + math.pi)/(math.pi/4) + eps)) % 8
            sect.setdefault(k, []).append((i, j, g))
    W = {k: sum(g for _, _, g in v) + centre for k, v in sect.items()}
    taps = {(0, 0): sum(centre/W[k] for k in range(8))/8}
    for k, v in sect.items():
        for i, j, g in v: taps[(i, j)] = taps.get((i, j), 0) + g/W[k]/8
    return taps
def H(taps, fx, fy):
    fx = np.asarray(fx, float); fy = np.asarray(fy, float); out = np.zeros(np.broadcast(fx, fy).shape, complex)
    for (dx, dy), w in taps.items(): out += w*np.exp(2j*np.pi*(fx*dx + fy*dy))
    return out
if __name__ == '__main__':
    t = kuw(); print(len(t), sum(t.values()))