# Independent re-derivation of the linearised Kuwahara taps and other filters.
import numpy as np
from fractions import Fraction
def sector(i, j):
    # exact: angle theta=atan2(j,i) in (-pi,pi]; k=floor((theta+pi)/(pi/4)) mod 8, boundary -> higher sector
    # use octant logic with integer comparisons. phi = theta+pi in (0, 2pi]
    import math
    th = math.atan2(j, i)
    # exact boundary detection: axes or diagonals
    on_b = (i == 0 or j == 0 or abs(i) == abs(j))
    k = (th + math.pi) / (math.pi / 4)
    if on_b: k = round(k)
    else: k = math.floor(k)
    return int(k) % 8
def kuw(R=5, stride=1):
    r2 = R*R + .5; g = {}; sec = {}; 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[(i,j)] = np.exp(-2*d/r2)
            if d == 0:
                W += g[(i,j)]; sec[(i,j)] = list(range(8))
            else:
                k = sector(i, j); W[k] += g[(i,j)]; sec[(i,j)] = [k]
    taps = {}
    for q, gq in g.items():
        taps[(q[0]*stride, q[1]*stride)] = sum(gq / W[k] for k in sec[q]) / 8
    return taps
def gauss_centre(c, s, stride=1):
    R = int(np.ceil(3*s)); t = {}
    for j in range(-R, R+1):
        for i in range(-R, R+1): t[(i*stride, j*stride)] = np.exp(-(i*i+j*j)/(2*s*s))
    S = sum(t.values()); t = {k: (1-c)*v/S for k, v in t.items()}; t[(0,0)] += c; return t
def box(n): return {(i,j): 1/n/n for i in range(n) for j in range(n)}
def H(t, fx, fy):
    fx = np.asarray(fx, float); fy = np.asarray(fy, float); out = 0
    for (dx, dy), w in t.items(): out = out + w*np.exp(-2j*np.pi*(fx*dx+fy*dy))
    return out
if __name__ == '__main__':
    import sys; import os; sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..'))
    import filters as F
    t = kuw(5); r = {(a,b): w for a,b,w in F.kuwahara(5)}
    print(len(t), len(r), sum(t.values()), max(abs(t[k]-r[k]) for k in t), set(t)==set(r))
    print('sum h^3/12 sqrt', np.sqrt(sum(w**3 for w in t.values())/12))
    for k in (2,):
        t2 = kuw(2); r2_ = {(a,b): w for a,b,w in F.kuwahara(2)}; print('R2', max(abs(t2[k]-r2_[k]) for k in t2))