ldnoise.py

#%% <- 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)