# The banding study: every dither candidate through Starlock's glow path, on dark test scenes, # at each painterly level and pixel ratio, measured against the float ideal. # python3 experiment.py -> results.json (all numbers) and a table on stdout # python3 experiment.py --quick -> fewer cases, for a fast look # python3 experiment.py --summary -> the entry's main table (out/summary.txt) import sys, json, itertools import numpy as np from multiprocessing import Pool import pipeline as P import dithers as D FRAME = 640 # render pixels on a side ORIGIN = (60., 1100.) # where the crop sits on the screen (gl_FragCoord of its corner): mid-phone PPD = {1: 45., 2: 64.} # pixels per degree: a 109 ppi desktop at 60 cm; a phone at 2x, 30 cm VIOLET, EMBER, WHITE = (.55, .3, 1.), (.75, .45, 1.), (1., .94, 1.) # --- scenes: lists of layers (kind, centre x, y, size w, h, colour, opacity, extra) in frame px --- def scenes(): c = FRAME / 2 pent = [(c + 120 * np.cos(2 * np.pi * k / 5 + .3), c + 120 * np.sin(2 * np.pi * k / 5 + .3)) for k in range(5)] return { 'glow': [('glow', c, c, 620, 620, VIOLET, .12, 0.)], 'faint': [('glow', c, c, 620, 620, VIOLET, .04, 0.)], 'stack': [('glow', x, y, 420 + 40 * k, 420 + 40 * k, (VIOLET, EMBER)[k % 2], .05 + .015 * k, 0.) for k, (x, y) in enumerate(pent)], 'concentric': [('glow', c, c, 380 + 60 * k, 380 + 60 * k, (VIOLET, EMBER)[k % 2], .05 + .015 * k, 0.) for k in range(5)], 'ring': [('ring', c, c, 620, 620, WHITE, .15, (.06, 0.)), ('ring', c, c, 420, 420, VIOLET, .3, (.1, .3))], } def layer_alpha(kind, u, v, extra): return P.glow_alpha(u, v, extra) if kind == 'glow' else P.ring_alpha(u, v, *extra) # --- candidates ------------------------------------------------------------------------------- # field: the noise; pdf: rpdf (+-0.5 LSB), tpdf (two fields, +-1 LSB, rail-guarded), tri (one field # remapped, rail-guarded); gate: Starlock's step(.002, a) cutoff; seeds: a different offset per layer; # css: noise evaluated per CSS pixel; model: how sin is computed (dithers.sin_hash); out: a dither on the final 8-bit write; target: u8 or f16. BASE = dict(field='sin', pdf='rpdf', gate=True, seeds=False, css=False, out=None, target='u8') CANDIDATES = { 'none': dict(BASE, field=None), 'C0 sin (now)': dict(BASE), 'C0 sin exact': dict(BASE, model='M1'), 'C1 sin, no gate': dict(BASE, gate=False), 'C2 sin, seeded': dict(BASE, seeds=True), 'C3 ign': dict(BASE, field='ign', gate=False, seeds=True), 'C3 ign shared': dict(BASE, field='ign', gate=False), 'C3 ign, gate': dict(BASE, field='ign', seeds=True), 'C1 sin, no gate, seeded': dict(BASE, gate=False, seeds=True), 'white, no gate, seeded': dict(BASE, field='white', gate=False, seeds=True), 'C4 blue': dict(BASE, field='blue', gate=False, seeds=True), 'C5 blue tpdf': dict(BASE, field='blue', pdf='tpdf', gate=False, seeds=True), 'C6 ign tri': dict(BASE, field='ign', pdf='tri', gate=False, seeds=True), 'C9 blue css': dict(BASE, field='blue', gate=False, seeds=True, css=True), 'C7 now+out': dict(BASE, out='blue'), 'C7 ign+out': dict(BASE, field='ign', gate=False, seeds=True, out='blue'), 'C7 blue+out': dict(BASE, field='blue', gate=False, seeds=True, out='blue'), 'C7 ign+ign out': dict(BASE, field='ign', gate=False, seeds=True, out='ign'), 'SHIP ign+ign-swap out': dict(BASE, field='ign', gate=False, seeds=True, out='ign-swap'), 'C7 ign+r2 out': dict(BASE, field='ign', gate=False, seeds=True, out='r2'), 'C8 f16+out': dict(BASE, field=None, target='f16', out='blue'), 'C8 f16': dict(BASE, field=None, target='f16'), } def noise(field, x, y, k, seeds, sin_model='M2'): if seeds: # integer offsets per layer: a different stretch of the same field x = x + 37 * k; y = y + 91 * k if field == 'sin': return D.sin_hash(x, y, sin_model) if field == 'ign': return D.ign(x, y) if field == 'blue': return D.blue(x, y) if field == 'white': return D.white(x, y, k) raise ValueError(field) def render(scene, cand, R, dpr): """The final 8-bit codes (H,W,3) for one case.""" sin_model = cand.get('model', 'M2') W = H = FRAME gx, gy = P.coords(W, H) gx = gx + np.float32(ORIGIN[0]); gy = gy + np.float32(ORIGIN[1]) nx, ny = (np.floor(gx / dpr) + .5, np.floor(gy / dpr) + .5) if cand['css'] else (gx, gy) t = np.zeros((H, W, 3)) for k, (kind, cx, cy, w, h, col, op, extra) in enumerate(scene): u, v, inside = P.quad_uv(W, H, cx, cy, w, h) a = layer_alpha(kind, u, v, extra) rgb = np.array(col)[None, None] * (a * op)[..., None] if cand['field'] and cand['target'] == 'u8': u1 = noise(cand['field'], nx, ny, k, cand['seeds'], sin_model) if cand['pdf'] == 'rpdf': d = u1 - .5 else: tp = (u1 - .5) + (noise(cand['field'], nx + 23, ny + 41, k + 7, True, sin_model) - .5) if cand['pdf'] == 'tpdf' else D.tri(u1) s = np.clip(255 * rgb.max(-1) - .5, 0, 1) # the rail guard: RPDF near black, TPDF above d = (1 - s) * (u1 - .5) + s * tp if cand['gate']: d = d * ((a * op) >= .002) rgb = rgb + (d / 255.)[..., None] if cand['target'] == 'u8': t = np.where(inside[..., None], P.to_unorm8(t / 255. + np.clip(rgb, 0, 1)), t) else: t = np.where(inside[..., None], (t + np.maximum(rgb, 0)).astype(np.float16).astype(np.float64), t) tex = t / 255. if cand['target'] == 'u8' else t disp = P.srgb(P.aces(P.kuwahara(tex, R, dpr))) if cand['out']: # the paint pass's own dither, in display codes (offset from the layers' noise) u = {'blue': lambda: D.blue(gx, gy, 64, 2), 'ign': lambda: D.ign(gx + 29, gy + 53), 'ign-swap': lambda: D.ign(gy, gx), 'r2': lambda: D.r2(gx, gy)}[cand['out']]() disp = disp + ((u - .5) / 255.)[..., None] return P.to_unorm8(disp) def ideal(scene, R, dpr): W = H = FRAME t = np.zeros((H, W, 3)) for kind, cx, cy, w, h, col, op, extra in scene: u, v, inside = P.quad_uv(W, H, cx, cy, w, h) t += np.where(inside[..., None], np.array(col)[None, None] * (layer_alpha(kind, u, v, extra) * op)[..., None], 0) return P.srgb(P.aces(P.kuwahara(t, R, dpr))) # --- metrics ---------------------------------------------------------------------------------- def luminance(disp): c = np.clip(disp, 0, 1) lin = np.where(c <= 0.04045, c / 12.92, ((c + 0.055) / 1.055) ** 2.4) return lin @ np.array([0.2126, 0.7152, 0.0722]) def lstar_of(Y): return np.where(Y > 216 / 24389, 116 * np.cbrt(Y) - 16, Y * 24389 / 27) def blur(img, sigma): r = int(np.ceil(4 * sigma)); k = np.exp(-np.arange(-r, r + 1) ** 2 / (2 * sigma * sigma)); k /= k.sum() p = np.pad(img, r, mode='reflect') p = np.apply_along_axis(lambda a: np.convolve(a, k, 'valid'), 0, p) return np.apply_along_axis(lambda a: np.convolve(a, k, 'valid'), 1, p) def flat_steps(Lcode): """Contour pixels: a step between two neighbours with a flat run of 8 on each side.""" n = 0 for a in (Lcode, Lcode.T): d = np.diff(a, axis=1) step = d != 0 flat = np.ones_like(step) for s in range(1, 8): flat[:, s:] &= d[:, :-s] == 0 flat[:, :-s] &= d[:, s:] == 0 n += (step & flat).sum() return n def metrics(out8, ideal_disp, dpr): B = 24 # leave out the frame's edge Yo, Yi = luminance(out8 / 255.), luminance(ideal_disp) sb, sv = .187 * PPD[dpr] / 2, .187 * PPD[dpr] / 8 # blurs that halve 2 and 8 cycles per degree # (the eye's optics average light, not lightness: blur luminance, then take L*) eb = lstar_of(blur(Yo, sb)) - lstar_of(blur(Yi, sb)) ev = lstar_of(blur(Yo, sv)) - lstar_of(blur(Yi, sv)) crop = lambda a: a[B:-B, B:-B] code = (out8 @ np.array([1, 256, 65536])).astype(np.int64) return dict( band=float(np.sqrt(np.mean(crop(eb) ** 2))), # low-frequency error: bands and blotches band99=float(np.percentile(np.abs(crop(eb)), 99)), grain=float(np.sqrt(np.mean(crop(ev - eb) ** 2))), # visible-scale noise left over steps=float(flat_steps(crop(code)) / (crop(code).size / 1e6)), # contour pixels per megapixel lift=float(crop(out8)[crop(ideal_disp.max(-1)) < .5 / 255].mean()) if (crop(ideal_disp.max(-1)) < .5 / 255).any() else 0., ) def case(args): sname, cname, R, dpr = args sc = scenes()[sname] return (sname, cname, R, dpr, metrics(render(sc, CANDIDATES[cname], R, dpr), IDEAL[(sname, R, dpr)], dpr)) IDEAL = {} def _init(keys): for sname, R, dpr in keys: IDEAL[(sname, R, dpr)] = ideal(scenes()[sname], R, dpr) SUMMARY = ['none', 'C0 sin (now)', 'C1 sin, no gate', 'C3 ign, gate', 'C3 ign', 'C4 blue', 'SHIP ign+ign-swap out', 'C8 f16+out'] def summary(): """The entry's main table: each scene at full paint and with paint off, pixel ratio 2.""" keys = [(s, R, 2) for s in scenes() for R in (5, 0)] jobs = [(s, c, R, d) for s, R, d in keys for c in SUMMARY] with Pool(4, initializer=_init, initargs=(keys,)) as pool: res = {(s, c, R): m for s, c, R, d, m in pool.map(case, jobs, chunksize=2)} print('band = low-frequency error in L* (lower is better); steps = contour pixels per megapixel') for R in (5, 0): print(f'\npaint {"full" if R else "off"} (radius {R}), pixel ratio 2') print(f'{"":24s}' + ''.join(f'{s:>18s}' for s in scenes())) for c in SUMMARY: print(f'{c:24s}' + ''.join(f"{res[(s, c, R)]['band']:9.3f} {res[(s, c, R)]['steps']:7.0f} " for s in scenes())) if __name__ == '__main__': if '--summary' in sys.argv: summary(); sys.exit() quick = '--quick' in sys.argv S = ['glow', 'stack'] if quick else list(scenes()) C = list(CANDIDATES) Rs, Ds = (0, 2, 5), (1, 2) keys = list(itertools.product(S, Rs, Ds)) jobs = list(itertools.product(S, C, Rs, Ds)) with Pool(4, initializer=_init, initargs=(keys,)) as pool: res = pool.map(case, jobs, chunksize=4) out = [dict(scene=s, cand=c, R=R, dpr=d, **m) for s, c, R, d, m in res] json.dump(out, open('results.json' if not quick else 'results-quick.json', 'w'), indent=0) for s in S: for d in Ds: print(f'\n{s}, pixel ratio {d}: band (L*), grain (L*), contour px/MP at paint off | lite | full') for c in C: row = [m for m in out if m['scene'] == s and m['cand'] == c and m['dpr'] == d] row.sort(key=lambda m: m['R']) print(f' {c:24s}' + ' | '.join(f"{m['band']:.3f} {m['grain']:.3f} {m['steps']:7.0f}" for m in row))