#%% <- this is vscode stuff, this will run on the command line but plots might come out weird import numpy as np u32 = np.uint32 i32 = np.int32 def golden_ratio_sequence(i: u32) -> u32: return u32(i) * u32(2654435769) # 0.614... in 0.32 fixed point def reverse_bits32(x: u32) -> u32: x = u32(x) x = ((x >> u32(1)) & u32(0x55555555)) | ((x & u32(0x55555555)) << u32(1)) x = ((x >> u32(2)) & u32(0x33333333)) | ((x & u32(0x33333333)) << u32(2)) x = ((x >> u32(4)) & u32(0x0F0F0F0F)) | ((x & u32(0x0F0F0F0F)) << u32(4)) x = ((x >> u32(8)) & u32(0x00FF00FF)) | ((x & u32(0x00FF00FF)) << u32(8)) x = ( x >> u32(16) ) | ( x << u32(16)) return x def nested_uniform_scramble(x: u32) -> u32: x = reverse_bits32(x) x ^= x * u32(0x6c50b47c) x ^= x * u32(0xb82f1e52) x ^= x * u32(0xc7afe638) x ^= x * u32(0x8d22f6e6) return reverse_bits32(x) def okay_blue_noise(i: u32) -> u32: return golden_ratio_sequence(nested_uniform_scramble(i)) #%% import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = (9,2) def spectrum(seq): S = np.fft.rfft(2 * seq - 1.0) return np.abs(S) / len(seq) def plots(name: str, seq_u32: u32): seq = seq_u32 * np.ldexp(1, -32) figure, (histogram_axis, spectrum_axis) = plt.subplots(1, 2) histogram_axis.hist(seq, 128) histogram_axis.set_xlabel(f"{name} histogram: {len(seq)} points") dft = spectrum(seq) spectrum_axis.plot(dft) spectrum_axis.set_xlabel(f"{name} spectrum: {len(seq)} points") n = 3333 #%% from numpy.random import default_rng rng = default_rng() plots("PRNG", rng.integers(0, 0x1_0000_0000, size=n)) #%% i = np.arange(n) plots("golden ratio sequence", golden_ratio_sequence(i)) #%% plots("nested uniform scramble-shuffled grs", okay_blue_noise(i)) #%% def xorshift(x): x = u32(x) x ^= x << u32(13) x ^= x >> u32(17) x ^= x << u32(5) return x i = np.arange(16) print(i) print(xorshift(i) & 15) #%% print(xorshift(i + 64 + 16) & 15) #%% def xorshift_star(x): x = u32(x) x ^= x << u32(13) x ^= x >> u32(17) x ^= x << u32(5) x *= u32(0x9e02ad0d) return x print(xorshift_star(i + 64 + 16) & 15) #%% print(len(np.unique(xorshift_star(np.arange(0xffff)) & 0xffff)), 0xffff) print(len(np.unique(xorshift_star(np.arange(0x1ffff)) & 0x1ffff)), 0x1ffff) #%% def all_ones_below_high_bit(x: u32) -> u32: x = u32(x) x |= (x >> u32(16)) x |= (x >> u32(8)) x |= (x >> u32(4)) x |= (x >> u32(2)) x |= (x >> u32(1)) # note this last shift: we don't include the high bit return x >> u32(1) def unfolded_masked_xorshift(x: u32, cap_mask: u32) -> u32: mask = all_ones_below_high_bit(x & cap_mask) upper = x & ~mask lower = xorshift_star(x) & mask result = upper + lower return result def masked_xorshift(x: u32, bits: u32 = 8) -> u32: # all ones if (x & 0x100) == 0x100, all zeros otherwise sign_mask = i32(x << u32(31 - bits)) >> i32(31) sign_mask = u32(sign_mask) cap_mask = u32((1 << bits) - 1) return unfolded_masked_xorshift(x ^ sign_mask, cap_mask) ^ sign_mask i = np.arange(1024) plt.plot(i - masked_xorshift(i)) #%% def white_shuffle(i: u32) -> u32: s = i s = masked_xorshift(s) s = nested_uniform_scramble(s) return s def white(i: u32) -> u32: return golden_ratio_sequence(white_shuffle(i)) i = np.arange(3333) plots("white", white(i)) #%% plt.rcParams['image.cmap'] = 'gray' plt.rcParams['image.interpolation'] = 'none' px = 1/plt.rcParams['figure.dpi'] n = 512; i = np.arange(n*n) f, ax = plt.subplots(figsize=(n*px, n*px)); ax.axis('off') ax.imshow(white(i).reshape((n,n)) * np.ldexp(1, -32)) #%% def white_shuffle(i: u32) -> u32: s = i s = nested_uniform_scramble(s) s = masked_xorshift(s) s = nested_uniform_scramble(s) return s def white(i: u32) -> u32: return golden_ratio_sequence(white_shuffle(i)) f, ax = plt.subplots(figsize=(n*px, n*px)); ax.axis('off') ax.imshow(white(i).reshape((n,n)) * np.ldexp(1, -32)) #%% def kronecker_sequence(i: u32, a: u32) -> u32: return u32(i) * u32(a) def blue(i: u32) -> u32: s = white_shuffle(i >> 1) b = kronecker_sequence(s, 2654435770) # 0.31 fixed point golden ratio odd = u32(i & 1) b ^= (odd ^ u32(1)) - u32(1) # negate on odd indices b += odd b ^= b >> u32(6) # gray code round return b i = np.arange(3333) plots("blue", blue(i)) #%% def spectrum_2d(img): dft = np.abs(np.fft.fftshift(np.fft.fft2(img))) dft /= dft.shape[0] return np.clip(dft, 0, 1) def lookit(image): f, ax = plt.subplots(figsize=(2*image.shape[0]*px, image.shape[1]*px)) ax.axis('off') ax.imshow(np.hstack((image, spectrum_2d(image)))) n = 384 i = np.arange(n*n).reshape((n,n)) noise = blue(i) * np.ldexp(1, -32) lookit(noise) #%% def spiral(n: u32, lo: float, hi: float): x, y = np.meshgrid(np.linspace(lo, hi, n), np.linspace(lo, hi, n)) # two sqrts: one for the distance, two to adjust for spirals closer to 0 being tighter (think random sampling in a circle) xy = np.round(np.sqrt(np.sqrt(x*x + y*y)) * np.sqrt(n*n + n*n)) # rescaled to a reasonable range that makes debugging possible without a third eye angles = (np.arctan2(y, x) + np.pi) / (2 * np.pi) # sort by magnitudes then angles (I don't know why lexsort is little endian), then invert the sort spiral = np.lexsort((angles.flatten(), xy.flatten())).argsort() return spiral.reshape((n,n)) s = spiral(n, -1.0, 1.0) noise = blue(s) * np.ldexp(1, -32) lookit(noise) #%% s = spiral(n, 2.0, 4.0) noise = blue(s) * np.ldexp(1, -32) lookit(noise) #%% spiral(8, 2.0, 4.0) #%% def left_shift_2(x: u32) -> u32: x = (x ^ (x << 16)) & 0x0000ffff x = (x ^ (x << 8)) & 0x00ff00ff x = (x ^ (x << 4)) & 0x0f0f0f0f x = (x ^ (x << 2)) & 0x33333333 x = (x ^ (x << 1)) & 0x55555555 return u32(x) def z_order(x: u32, y: u32) -> u32: return left_shift_2(x) + (left_shift_2(y) << u32(1)) tile_bits = 6 tile_n = 1 << tile_bits tile_mask = tile_n - 1 tile_path = spiral(tile_n, 2.0, 4.0) def blue_2d(x: u32, y: u32) -> u32: x_lo = x & tile_mask y_lo = y & tile_mask x_hi = x >> tile_bits y_hi = y >> tile_bits tile = z_order(x_hi, y_hi) i = (tile << u32(2*tile_bits)) + tile_path[y_lo, x_lo] return blue(i) x, y = np.meshgrid(np.arange(n), np.arange(n)) noise = blue_2d(x, y) * np.ldexp(1, -32) lookit(noise) #%% plt.hist(noise.flatten(), 384) pass #%% s = (spiral(384, -0.08, 0.08) & 7) / 8 lookit(s) #%% s = u32(spiral(384, -0.08, 0.08)) s = masked_xorshift(s, 2) s ^= s >> 1 s = reverse_bits32(s) s = nested_uniform_scramble(s) lookit(s / s.max()) #%% f32 = np.float32 def uniform_to_triangle_dist(x): # From demofox @ https://www.shadertoy.com/view/4t2SDh x = (x + 0.5) % 1 orig = x * 2.0 - 1.0 nz = orig != 0 x[~nz] = -1 x[nz] = orig[nz] / np.sqrt(np.abs(orig[nz])) x = x - np.sign(orig) + 0.5 x = (x - 0.5) * 0.5 + 0.5 return x def quantize_dithered(img: f32, noise: f32, dither_bits: int, output_bits: int) -> f32: output_scale = 2**output_bits dither_scale = 2**dither_bits / output_scale img_plus_noise = img + dither_scale * (noise - 0.5) quantized = np.round(img_plus_noise * (output_scale - 1)) / (output_scale - 1) return np.clip(quantized, 0, 1) def dither_gradient(dither_bits, output_bits): n = 512 m = 100 gradient = np.tile(np.linspace(0,1,n), (m*5, 2)) vac_texture = np.concatenate((plt.imread("HDR_L_0.png").T, plt.imread("HDR_L_0.png").T)).T x, y = np.meshgrid(np.arange(n), np.arange(m)) random = rng.random((m,n)) lds_white = white(y*m + x) * np.ldexp(1, -32) lds_blue = blue_2d(x, y) * np.ldexp(1, -32) vac_blue = vac_texture[:m,:n] random_tri = 0.5*(rng.random((m,n)) + rng.random((m,n))) lds_white_tri = uniform_to_triangle_dist(lds_white) lds_blue_tri = uniform_to_triangle_dist(lds_blue) vac_blue_tri = uniform_to_triangle_dist(vac_texture[:m,:n]) noise = np.zeros((m*5, n*2)) noise[0*m:1*m//2, 0:2*n] = 0.5 noise[1*m:2*m, 0:1*n] = random noise[1*m:2*m, n:2*n] = random_tri noise[2*m:3*m, 0:1*n] = lds_white noise[2*m:3*m, n:2*n] = lds_white_tri noise[3*m:4*m, 0:1*n] = lds_blue noise[3*m:4*m, n:2*n] = lds_blue_tri noise[4*m:5*m, 0:1*n] = vac_blue noise[4*m:5*m, n:2*n] = vac_blue_tri image = quantize_dithered(gradient, noise, dither_bits, output_bits) image[1*m//2:1*m, 0:2*n] = gradient[1*m//2:1*m, 0:2*n] f, ax = plt.subplots(figsize=(image.shape[1]*px, image.shape[0]*px)); ax.axis('off') ax.imshow(image) dither_gradient(2,4) # %% def red(i: u32) -> u32: s = white_shuffle(i >> 1) r = kronecker_sequence(s, 2654435770) # 0.31 fixed point golden ratio r[(i & 1) == 1] ^= r[(i & 1) == 1] >> u32(6) return r def red_2d(x: u32, y: u32) -> u32: i = left_shift_2(x) + (left_shift_2(y) >> u32(1)) return red(i) noise = red_2d(x, y) * np.ldexp(1, -32) lookit(noise) #%% plt.hist(noise.flatten(), 384) pass #%% n = 64 border = 2 gap = border + 8 rows = 4 cols = 3 img = np.zeros((rows*(border-1) + n*rows + (rows-1)*gap, cols*(border-1) + 1 + n*cols + (cols-1)*gap)) mask = np.zeros_like(img) x = border y = border for i in range(rows): offset = rng.integers(0x1000)*n*n b = border x = b img[y-b:y+n+b, x-b:x+n+b] = 0.5 img[y:y+n, x:x+n] = f32(rng.random(size=(n,n))) < 0.5 mask[y-b:y+n+b, x-b:x+n+b] = 1 x += n + gap img[y-b:y+n+b, x-b:x+n+b] = 0.5 img[y:y+n, x:x+n] = (white(offset +np.arange(n*n)).reshape((n,n)) * np.ldexp(1, -32)) < 0.5 mask[y-b:y+n+b, x-b:x+n+b] = 1 x += n + gap img[y-b:y+n+b, x-b:x+n+b] = 0.5 bx, by = np.meshgrid(offset + np.arange(n), offset + np.arange(n)) img[y:y+n, x:x+n] = (blue_2d(bx, by) * np.ldexp(1, -32)) < 0.5 mask[y-b:y+n+b, x-b:x+n+b] = 1 y += n + gap img_rgba = np.zeros(img.shape + (4,)) img_rgba[...,0] = img img_rgba[...,1] = img img_rgba[...,2] = img img_rgba[...,3] = mask def scale_image(img, scale): return np.repeat(np.repeat(img, scale, axis=0), scale, axis=1) img_rgba = scale_image(img_rgba, 3) f, ax = plt.subplots(figsize=(img_rgba.shape[0]*px, img_rgba.shape[1]*px)); ax.axis('off') ax.imshow(img_rgba)