import numpy as np
from PIL import Image
from scipy import ndimage as ndi
SEED = 7
N = 2000
rng = np.random.default_rng(SEED)
GROUND = (0x10, 0x0F, 0x0F)
INKS = [(0xD1,0x4D,0x41),(0xDA,0x70,0x2C),(0xD0,0xA2,0x15),(0x87,0x9A,0x39),
        (0x3A,0xA9,0x9F),(0x43,0x85,0xBE),(0x8B,0x7E,0xC8),(0xCE,0x5D,0x97),(0xFF,0xFC,0xF0)]
M = 110  # margin >5%
yy, xx = np.mgrid[0:N, 0:N].astype(np.float64)
px, py = xx.copy(), yy.copy()
# faults: lines through the frame fanning around a main diagonal; slip magnitudes heavy-tailed
# (a few huge slips, many tiny ones). Each slip shears everything on its far side along the line.
faults = []
base = np.deg2rad(-38)
mags = np.sort(rng.pareto(1.1, 14) + 0.15)[::-1]
mags = np.minimum(mags / mags[0], 1.0) * 400 + 14
for i, m in enumerate(mags):
    a = base + rng.normal(0, 0.12)
    c = rng.uniform(0.15, 0.85) * N if i else 0.42 * N
    d = np.array([np.cos(a), np.sin(a)])
    nrm = np.array([-d[1], d[0]])
    p0 = np.array([c, c * 0.9 + rng.uniform(-250, 250)])
    faults.append((p0, d, nrm, m * (1 if rng.random() < .5 else -1)))
gap = np.zeros((N, N), bool)
for p0, d, nrm, m in faults:
    s = (px - p0[0]) * nrm[0] + (py - p0[1]) * nrm[1]
    far = s > 0
    gap |= np.abs(s) < (8 + 9 * max(0.0, 1.0 - abs(m) / 420.0))
    px = np.where(far, px - d[0] * m, px)
    py = np.where(far, py - d[1] * m, py)
# stripes: bands across a fixed direction, widths drawn heavy-tailed, colours cycle
th = np.deg2rad(-72)
u = px * np.cos(th) + py * np.sin(th)
widths = np.minimum(rng.pareto(1.3, 80) + 0.2, 6) * 38 + 14
edges = np.cumsum(widths) - 2200
idx = np.searchsorted(edges, u)
pal = np.array([5,4,1,0,2,8,5,4,6])
lab = pal[(idx * 7 + idx // 9) % len(pal)]
# quiet area: one big slab, cut by the same faults, so its outline is torn too; the rest is ground
c0, c1 = np.cos(base), np.sin(base)
CX = CY = 1000.0
for _ in range(6):
    t = (px - CX) * c0 + (py - CY) * c1
    nn = -(px - CX) * c1 + (py - CY) * c0
    mask = (np.abs(t) < 780) & (nn > -280) & (nn < 360)
    ys, xs = np.nonzero(mask)
    CX += 1000 - xs.mean(); CY += 1000 - ys.mean()
mask &= (xx > M) & (xx < N - M) & (yy > M) & (yy < N - M)
mask &= ~gap
img = np.zeros((N, N), np.int16) - 1
img[mask] = lab[mask]
# enforce min mark width (>=8 px): per-colour opening, drop what is thinner
st = np.ones((9, 9), bool)
out = np.zeros((N, N), np.int16) - 1
for k in range(len(INKS)):
    m = ndi.binary_opening(img == k, structure=st)
    out[m] = k
# drop small detached fragments: any isolated piece under 600 px reads as a stray
# speck, not a fault offset; slab-scale pieces and gap-adjacent offsets are kept
for k in range(len(INKS)):
    lab, n = ndi.label(out == k)
    sizes = np.bincount(lab.ravel())
    out[np.isin(lab, [i for i in range(1, n + 1) if sizes[i] < 600])] = -1
# drop slivers adrift outside the mass: keep a small piece only when slab-scale
# mass flanks it on both cross-fault sides (fault-channel inlay); lone offsets
# floating past the torn edges read as strays, not as part of the gesture
st4 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
uni, un = ndi.label(out >= 0, structure=st4)
usizes = np.bincount(uni.ravel())
big = np.isin(uni, [i for i in range(1, un + 1) if usizes[i] >= 20000])
T = xx * c0 + yy * c1   # along-fault
S = -xx * c1 + yy * c0  # cross-fault
tb, sb = T[big], S[big]
for i in range(1, un + 1):
    if usizes[i] >= 20000:
        continue
    comp = uni == i
    sel = (tb >= T[comp].min() - 50) & (tb <= T[comp].max() + 50)
    sc = S[comp]
    if not (sel.any() and sb[sel].min() < sc.min() and sb[sel].max() > sc.max()):
        out[comp] = -1
rgb = np.zeros((N, N, 3), np.uint8); rgb[:] = GROUND
for k, c in enumerate(INKS):
    rgb[out == k] = c
Image.fromarray(rgb).save("render-1.png")