# A blue-noise tile by Ulichney's void-and-cluster method (1993): a 64x64 grid of ranks
# 0..4095 where every rank's pixels are spread as evenly as they can be. Deterministic
# (seeded), numpy only. Prints nothing; bluenoise(n) returns ranks / (n*n) in [0, 1).
import numpy as np

def _energy_kernel(n, sigma=1.5):
    d = np.minimum(np.arange(n), n - np.arange(n)).astype(float)
    g = np.exp(-(d[:, None] ** 2 + d[None, :] ** 2) / (2 * sigma * sigma))
    return np.fft.rfft2(g)

def _blur(m, K):
    return np.fft.irfft2(np.fft.rfft2(m) * K, s=m.shape)

def bluenoise(n=64, seed=1, sigma=1.5):
    K = _energy_kernel(n, sigma)
    rng = np.random.default_rng(seed)
    N = n * n
    ones = max(1, N // 10)
    pat = np.zeros((n, n), bool)
    pat.flat[rng.choice(N, ones, replace=False)] = True
    # relax: move the tightest cluster's point into the largest void until stable
    while True:
        e = _blur(pat.astype(float), K)
        c = np.argmax(np.where(pat, e, -np.inf)); pat.flat[c] = False
        e = _blur(pat.astype(float), K)
        v = np.argmin(np.where(pat, np.inf, e)); pat.flat[v] = True
        if v == c: break
    rank = np.zeros(N, int)
    proto = pat.copy()
    # phase 1: ranks below the prototype's count, removing tightest clusters
    p = proto.copy()
    for r in range(ones - 1, -1, -1):
        e = _blur(p.astype(float), K)
        c = np.argmax(np.where(p, e, -np.inf)); p.flat[c] = False; rank[c] = r
    # phase 2 and 3: fill the largest voids
    p = proto.copy()
    for r in range(ones, N):
        e = _blur(p.astype(float), K)
        v = np.argmin(np.where(p, np.inf, e)); p.flat[v] = True; rank[v] = r
    return (rank.reshape(n, n) + 0.5) / N

if __name__ == '__main__':
    b = bluenoise()
    print(b.shape, b.min(), b.max())