import numpy as np from PIL import Image, ImageDraw from scipy.spatial import cKDTree from scipy.sparse import coo_matrix from scipy.sparse.csgraph import dijkstra import os SEED = 7 SIZE = 2000 GENS = 14 N_PTS = 34000 N_SINKS = 5000 # Flexoki, as given BG = (0x10, 0x0F, 0x0F) GREY = {"800": (0x40,0x3E,0x3C), "700": (0x57,0x56,0x53), "600": (0x6F,0x6E,0x69), "500": (0x87,0x85,0x80), "300": (0xB7,0xB5,0xAC), "200": (0xCE,0xCD,0xC3), "100": (0xE6,0xE4,0xD9), "50": (0xF2,0xF0,0xE5)} # hue order around the wheel: red, orange, yellow, green, cyan, blue, purple, magenta (600 tone, 400 tone) HUES = [((0xAF,0x30,0x29),(0xD1,0x4D,0x41)), ((0xBC,0x52,0x15),(0xDA,0x70,0x2C)), ((0xAD,0x83,0x01),(0xD0,0xA2,0x15)), ((0x66,0x80,0x0B),(0x87,0x9A,0x39)), ((0x24,0x83,0x7B),(0x3A,0xA9,0x9F)), ((0x20,0x5E,0xA6),(0x43,0x85,0xBE)), ((0x5E,0x40,0x9D),(0x8B,0x7E,0xC8)), ((0xA0,0x2F,0x6F),(0xCE,0x5D,0x97))] rng = np.random.default_rng(SEED) S = SIZE ROOT = np.array([0.70, 0.64]) # the crowded off-center junction VOID_C, VOID_R = np.array([0.24, 0.27]), 0.21 # large dark area, kept empty def in_void(p): return np.hypot(p[:,0]-VOID_C[0], p[:,1]-VOID_C[1]) < VOID_R # points, denser near the root and along spiral arms pts = rng.random((N_PTS*3, 2)) d = pts - ROOT r = np.hypot(d[:,0], d[:,1]); th = np.arctan2(d[:,1], d[:,0]) arm = 0.5 + 0.5*np.cos(2*(th - 3.2*np.log(r+0.05))) keep = rng.random(len(pts)) < (0.45 + 0.55*arm) * (0.6 + 0.4*np.exp(-r*2.0)) dv_ = np.hypot(pts[:,0]-VOID_C[0], pts[:,1]-VOID_C[1]) + 0.03*np.sin(9*np.arctan2(pts[:,1]-VOID_C[1], pts[:,0]-VOID_C[0])) keep &= rng.random(len(pts)) < np.clip((dv_-VOID_R)/0.06, 0, 1) pts = pts[keep][:N_PTS] pts = np.vstack([ROOT, pts]); n = len(pts) tree = cKDTree(pts) dist, nb = tree.query(pts, k=8) I = np.repeat(np.arange(n), 7); J = nb[:,1:].ravel(); L = dist[:,1:].ravel() A = np.minimum(I,J); B = np.maximum(I,J) key = np.unique(A*n+B); A = key//n; B = key%n EL = np.hypot(*(pts[A]-pts[B]).T) # weak sinks inside every region sink_idx = rng.choice(np.arange(1,n), N_SINKS, replace=False) sink_w = np.zeros(n); sink_w[sink_idx] = 0.4 + rng.random(N_SINKS) cum = np.zeros(len(A)) gen_flux = [] for g in range(GENS): cond = (0.02 + cum/ (cum.max()+1e-9))**1.1 if g else np.ones(len(A)) noise = np.exp(rng.normal(0, 0.55 if g else 0.0, len(A))) cost = EL*noise/cond M = coo_matrix((np.r_[cost,cost], (np.r_[A,B], np.r_[B,A])), shape=(n,n)).tocsr() D, P = dijkstra(M, indices=0, return_predecessors=True) order = np.argsort(-np.where(np.isfinite(D), D, -1)) flow = sink_w.copy() for v in order: p = P[v] if p >= 0: flow[p] += flow[v] # edge flux emap = {} f = np.zeros(len(A)) eid = dict(zip((A*n+B).tolist(), range(len(A)))) for v in range(1, n): p = P[v] if p < 0: continue a, b = (v, p) if v < p else (p, v) f[eid[a*n+b]] += flow[v] gen_flux.append(f) cum += f/ (f.max()+1e-9) total = np.sum(gen_flux, axis=0) tn = total/total.max() mid = (pts[A]+pts[B])/2 ang = (np.arctan2(mid[:,1]-ROOT[1], mid[:,0]-ROOT[0]) % (2*np.pi))/(2*np.pi) # Shared endpoint widths make neighboring limbs meet continuously. # The width law softly approaches a fine vein below the former cutoff. limb_edges = np.where(tn > 0.001)[0] raw_width = np.where(tn > 0.008, 4 + 95*tn**0.6, 0.0) node_width = np.zeros(n) np.maximum.at(node_width, A[limb_edges], raw_width[limb_edges]) np.maximum.at(node_width, B[limb_edges], raw_width[limb_edges]) # Short local relaxation spreads junction widths into their offshoots. # No random choices: routing, colors and subsequent jitter are unchanged. for _ in range(3): next_width = node_width.copy() decay = np.exp(-EL[limb_edges]*S/12.0) * 0.72 np.maximum.at(next_width, A[limb_edges], node_width[B[limb_edges]]*decay) np.maximum.at(next_width, B[limb_edges], node_width[A[limb_edges]]*decay) node_width = next_width img = Image.new("RGB", (S, S), BG) dr = ImageDraw.Draw(img) def px(p): return (p[0]*S, p[1]*S) for g in range(GENS): f = gen_flux[g]/gen_flux[g].max() t = g/(GENS-1) thresh = 0.0 if g == 0 else 0.0015 + 0.03*t**1.5 sel = np.where(f > thresh)[0] sel = sel[np.argsort(f[sel])] width_base = 4.0*(1-t)**1.3 + 1.0 spread = 4.0*(1-t)**1.2 + 0.4 # hue spread narrows with generations for e in sel: w = max(1, int(round(width_base*(0.4 + 1.8*np.sqrt(f[e]))))) hpos = ang[e]*7 + 0.6*g/GENS*8 * 0.4 + (rng.random()-0.5)*spread h = [4,5,6,7,0,1,2][int(np.floor(hpos)) % 7] if g < 2: col = HUES[h][0] # generation-0 twig fill uses full-intensity 600-tone hue arc elif g < GENS-4: col = HUES[h][0] if g < 5 and tn[e] < 0.3 else HUES[h][1] else: col = HUES[h][1] if t < 0.95 else GREY["100"] if g == GENS-1 and tn[e] > 0.5: col = GREY["50"] a = pts[A[e]]; b = pts[B[e]] dv = b-a; nv = np.array([-dv[1], dv[0]])/(np.hypot(*dv)+1e-9) off = nv*(rng.random()-0.5)*0.0035*(g/GENS+0.3) # slight per-generation offset dr.line([px(a+off), px(b+off)], fill=col, width=w) if w > 2: for q in (a+off, b+off): x, y = px(q); rr = w/2 dr.ellipse([x-rr, y-rr, x+rr, y+rr], fill=col) if g == 0: # primary limbs: thick, flat, saturated; colour set by angle, drifting with distance (no rng use) sel = limb_edges[np.maximum(node_width[A[limb_edges]], node_width[B[limb_edges]]) > 2.0] sel = sel[np.argsort(tn[sel])] rd = np.hypot(mid[:,0]-ROOT[0], mid[:,1]-ROOT[1]) for e in sel: w = raw_width[e] hpos = ang[e]*7 + rd[e]*2.2 h = [4,5,6,7,0,1,2][int(np.floor(hpos)) % 7] col = HUES[h][1] if tn[e] > 0.05 else HUES[h][0] a = pts[A[e]]; b = pts[B[e]] dv = b-a nv = np.array([-dv[1], dv[0]])/(np.hypot(*dv)+1e-9) ra = max(0.5, node_width[A[e]]/2) rb = max(0.5, node_width[B[e]]/2) dr.polygon([px(a+nv*ra/S), px(b+nv*rb/S), px(b-nv*rb/S), px(a-nv*ra/S)], fill=col) for q, rr in ((a, ra), (b, rb)): x, y = px(q) dr.ellipse([x-rr, y-rr, x+rr, y+rr], fill=col) img.save(os.path.join(os.path.dirname(os.path.abspath(__file__)), "render-1.png"))