spin_render.py

1import time23import numpy as np4from fractions import Fraction5from math import atan2, gcd, sqrt6from pathlib import Path7from PIL import Image89HERE = Path(__file__).resolve().parent1011N = 5512ODDS = list(range(1, N + 1, 2))13L = len(ODDS)14GOLDEN = 137.50776415RESOLUTIONS = [256, 512, 1024, 2048]16LAYER_COUNTS = [4, 8, 14, 28]17FIG_R = 102418PEAK_COUNT = 2019PEAK_RADIUS = 0.0120ZOOM_R = 102421ZOOM_W = 0.022223def primes_up_to_count(k):24    out = []25    c = 226    while len(out) < k:27        if all(c % d for d in range(2, int(sqrt(c)) + 1)):28            out.append(c)29        c += 130    return out3132def gaussian_angle(n):33    best = None34    for b in range(1, int(sqrt(n)) + 1):35        a2 = n - b * b36        a = int(round(sqrt(a2)))37        if a * a == a2 and a >= b >= 1:38            best = (a, b)39    if best is None:40        return 0.041    return np.degrees(atan2(best[1], best[0]))4243def schedules():44    ps = primes_up_to_count(L - 1)45    return {46        "unspun": [0.0] * L,47        "degrees": [float(k + 1) for k in range(L)],48        "golden": [GOLDEN * (k + 1) for k in range(L)],49        "primes": [0.0] + [float(p) for p in ps],50        "random": list(np.random.default_rng(1).uniform(0.0, 360.0, L)),51        "gaussian": [gaussian_angle(n) for n in ODDS],52    }5354def parity(z):55    return z - np.floor(z) >= 0.55657def stack(angles, r, layers=L, dtype=np.float32, box=(0.0, 0.0, 1.0)):58    x0, y0, w = box59    step = (np.arange(r, dtype=dtype) + dtype(0.5)) * dtype(w / r)60    dx = step + dtype(x0 - 0.5)61    dy = step + dtype(y0 - 0.5)62    acc = np.zeros((r, r), dtype=np.uint8)63    for k in range(layers):64        n = ODDS[k]65        h = 0.5 * n66        t = np.float64(angles[k]) * np.pi / 180.067        cx = dtype(h * np.cos(t))68        sx = dtype(h * np.sin(t))69        c0 = dtype(h * 0.5)70        a = c0 + cx * dx[None, :] + sx * dy[:, None]71        b = c0 - sx * dx[None, :] + cx * dy[:, None]72        acc += parity(a) & parity(b)73    return (np.float32(layers) - acc) / np.float32(layers)7475def disc_mask(r):76    d = (np.arange(r, dtype=np.float64) + 0.5) / r - 0.577    return d[:, None] ** 2 + d[None, :] ** 2 <= 0.257879def maximum_filter3(f):80    p = np.pad(f, 1, mode="constant", constant_values=-np.inf)81    r = f.shape[0]82    m = f.copy()83    for dy in (0, 1, 2):84        for dx in (0, 1, 2):85            np.maximum(m, p[dy:dy + r, dx:dx + r], out=m)86    return m8788def plateau(field, val, iy, ix, wide):89    r = field.shape[0]90    y0, y1 = max(0, iy - wide), min(r, iy + wide + 1)91    x0, x1 = max(0, ix - wide), min(r, ix + wide + 1)92    w = field[y0:y1, x0:x1] == val93    seed = np.zeros_like(w)94    seed[iy - y0, ix - x0] = True95    while True:96        g = seed.copy()97        g[1:] |= seed[:-1]98        g[:-1] |= seed[1:]99        g[:, 1:] |= seed[:, :-1]100        g[:, :-1] |= seed[:, 1:]101        g &= w102        if g.sum() == seed.sum():103            return np.nonzero(seed)[0] + y0, np.nonzero(seed)[1] + x0104        seed = g105106def peaks(field, mask, count=PEAK_COUNT, radius=PEAK_RADIUS):107    r = field.shape[0]108    cand = (field >= maximum_filter3(field)) & mask109    ys, xs = np.nonzero(cand)110    vs = field[ys, xs]111    order = np.lexsort((xs, ys, -vs))112    ys, xs, vs = ys[order], xs[order], vs[order]113    py = (ys + 0.5) / r114    px = (xs + 0.5) / r115    alive = np.ones(len(vs), dtype=bool)116    wide = int(radius * r) + 1117    out = []118    for _ in range(count):119        idx = int(np.argmax(alive))120        if not alive[idx]:121            break122        cy, cx = plateau(field, vs[idx], ys[idx], xs[idx], wide)123        out.append((float((cx + 0.5).mean() / r), float((cy + 0.5).mean() / r), float(vs[idx]), float(px[idx]), float(py[idx]), len(cx)))124        alive &= ~((px - px[idx]) ** 2 + (py - py[idx]) ** 2 <= radius ** 2)125    return out126127def drift(a, b, r):128    pa = np.array([[p[0], p[1]] for p in a])129    pb = np.array([[p[0], p[1]] for p in b])130    dd = np.sqrt(((pa[:, None, :] - pb[None, :, :]) ** 2).sum(2)).min(1)131    return float(np.median(dd)), float(dd.max()), int((dd * r <= 1.0).sum())132133def exact_variance(layers):134    s = [2 * k + 1 for k in range(layers)]135    tot = Fraction(0)136    for m in s:137        for n in s:138            d = gcd(m, n)139            tot += Fraction((d * d - 1) * (2 * (m - 1) * (n - 1) + d * d - 1), 16 * m * m * n * n)140    return tot / (layers * layers)141142def diagonal_maximum():143    cuts = sorted({Fraction(k, n) for n in ODDS for k in range(1, n + 1)} | {Fraction(0), Fraction(1)})144    best = 0145    rows = []146    for i in range(len(cuts) - 1):147        u = (cuts[i] + cuts[i + 1]) / 2148        sel = frozenset(n for n in ODDS if (n * u).numerator // (n * u).denominator % 2 == 1)149        if len(sel) > best:150            best, rows = len(sel), []151        if len(sel) == best:152            rows.append((cuts[i], cuts[i + 1], sel))153    cells = [(a, b, c, d) for (a, b, s1) in rows for (c, d, s2) in rows if s1 == s2]154    area = sum(((b - a) * (d - c) for (a, b, c, d) in cells), Fraction(0))155    return best, rows, cells, area, len(cuts)156157def dark_layers(angles, x, y):158    out = []159    for k in range(L):160        n = ODDS[k]161        t = np.float64(angles[k]) * np.pi / 180.0162        u = 0.5 + (x - 0.5) * np.cos(t) + (y - 0.5) * np.sin(t)163        v = 0.5 - (x - 0.5) * np.sin(t) + (y - 0.5) * np.cos(t)164        if parity(0.5 * n * u) and parity(0.5 * n * v):165            out.append(n)166    return out167168def centre_value(field, r):169    i = r // 2170    return float(field[i - 1:i + 1, i - 1:i + 1].mean())171172def save_png(field, path, side):173    r = field.shape[0]174    b = r // side175    q = field.reshape(side, b, side, b).mean((1, 3))176    q = np.round(np.round(q * L) * (255.0 / L)).astype(np.uint8)177    Image.fromarray(q, mode="L").save(path, optimize=True)178    return path.stat().st_size179180def contact_sheet(fields, path, side):181    sheet = np.zeros((2 * side, 2 * side), dtype=np.uint8)182    for i, f in enumerate(fields):183        r = f.shape[0]184        b = r // side185        q = f.reshape(side, b, side, b).mean((1, 3))186        q = np.round(np.round(q * L) * (255.0 / L)).astype(np.uint8)187        sheet[side * (i // 2):side * (i // 2) + side, side * (i % 2):side * (i % 2) + side] = q188    Image.fromarray(sheet, mode="L").save(path, optimize=True)189    return path.stat().st_size190191def frac(v):192    return "%d/%d" % (round(v * L), L)193194def report(name, angles):195    print("schedule", name, "angles", " ".join("%.6g" % a for a in angles))196    tops = {}197    for r in RESOLUTIONS:198        f = stack(angles, r)199        msk = disc_mask(r)200        ink = 1.0 - f201        d = ink[msk]202        pk = peaks(ink, msk)203        tops[r] = pk204        top = pk[0][2]205        area = float((d == d.max()).sum()) / (r * r)206        print(207            " R %4d peak %s locus %.6g mean %.6f rms %.6f centre %s top" % (r, frac(top), area, float(f[msk].mean()), float(d.std()), frac(1.0 - centre_value(f, r))),208            " ".join("(%.5f,%.5f,%s,%dpx)" % (q[0], q[1], frac(q[2]), q[5]) for q in pk[:3]),209        )210    for r in RESOLUTIONS[:-1]:211        med, mx, hit = drift(tops[r], tops[2 * r], r)212        print(" drift %4d->%4d median %.6f (%.2f px) max %.6f (%.2f px) within 1 px %d/%d" % (r, 2 * r, med, med * r, mx, mx * r, hit, PEAK_COUNT))213    x, y, v, sx, sy, area = tops[RESOLUTIONS[-1]][0]214    z = stack(angles, ZOOM_R, box=(x - ZOOM_W / 2, y - ZOOM_W / 2, ZOOM_W))215    zi = 1.0 - z216    lev = float((zi == zi.max()).mean())217    print(" zoom window %g wide at effective R %d: peak %s, cell at the top peak %.6g" % (ZOOM_W, int(ZOOM_R / ZOOM_W), frac(float(zi.max())), lev * ZOOM_W * ZOOM_W))218    dk = dark_layers(angles, sx, sy)219    print(" top peak (%.5f, %.5f) at %s, plateau %d px at R %d, %d layers dark there" % (x, y, frac(v), area, RESOLUTIONS[-1], len(dk)), " ".join(str(n) for n in dk))220    return tops221222def main():223    t0 = time.perf_counter()224    sched = schedules()225    print("layers", L, "odd scales 1..%d" % N, "peak radius", PEAK_RADIUS, "peak count", PEAK_COUNT)226    print("centre cell inradius 1/%d in every layer at every angle, scales dark at the centre" % (2 * N), " ".join(str(n) for n in dark_layers([0.0] * L, 0.5, 0.5)))227    for name, angles in sched.items():228        report(name, angles)229230    c, rows, cells, area, cuts = diagonal_maximum()231    print("unspun ink maximum by exact scan of %d breakpoints on the diagonal: %d/%d on %d intervals" % (cuts, c, L, len(rows)), " ".join("[%s, %s)" % (a, b) for a, b, _ in rows))232    print("unspun maximum locus: %d cells, total area %s = %.6g, one cell %s = %.6g, scales" % (len(cells), area, float(area), (rows[0][1] - rows[0][0]) ** 2, float((rows[0][1] - rows[0][0]) ** 2)), " ".join(str(n) for n in sorted(rows[0][2])))233234    print("fade at R %d, c = rms * sqrt(layers)" % FIG_R)235    msk = disc_mask(FIG_R)236    for name, angles in sched.items():237        row = []238        for lc in LAYER_COUNTS:239            f = stack(angles, FIG_R, lc)240            row.append((float(f.std()) * sqrt(lc), float(f[msk].std()) * sqrt(lc)))241        print(242            " %-9s square" % name, " ".join("%.6f" % a for a, b in row),243            " disc", " ".join("%.6f" % b for a, b in row),244        )245    print(" exact   square", " ".join("%.6f" % sqrt(float(exact_variance(lc)) * lc) for lc in LAYER_COUNTS))246247    f32 = stack(sched["golden"], 512)248    f64 = stack(sched["golden"], 512, dtype=np.float64)249    print("float32 against float64 at R 512 golden: max abs diff %.6g, pixels differing %d of %d" % (float(np.abs(f32 - f64).max()), int((f32 != f64).sum()), 512 * 512))250251    figs = {}252    for name, fname in (("unspun", "stack-unspun.png"), ("degrees", "stack-degrees.png"), ("primes", "stack-primes.png"), ("gaussian", "stack-gaussian.png")):253        f = stack(sched[name], FIG_R)254        f = np.where(disc_mask(FIG_R), f, 1.0).astype(np.float32)255        figs[name] = f256        print("figure", fname, "bytes", save_png(f, HERE / fname, 512))257    print("figure stack-sheet.png bytes", contact_sheet(list(figs.values()), HERE / "stack-sheet.png", 256))258    print("seconds", round(time.perf_counter() - t0, 2))259260main()