#!python
"""forest_floor: the scene for a painting. Half of it is above the ground and half below.

    forest_floor.py OUT.png [--seed N]

Above: a night wood, a moon, three morels. Below: soil, tree roots, and the fungus itself, which is the
pale threads. The mushrooms are only where it comes up for air. This picture is never shown as it is.
It is what brush.py looks at while it paints. Morel, Nov 2026.

What is true here and what isn't: morels do grow from a network of threads in the soil, and many fungi
trade with tree roots (water and minerals one way, sugar the other). Real threads are far too fine to
see like this, they don't glow, and how morels in particular deal with trees is still argued about.
"""
import argparse
import numpy as np
from PIL import Image, ImageFilter

ap = argparse.ArgumentParser(); ap.add_argument("out"); ap.add_argument("--seed", type=int, default=11)
a = ap.parse_args()
rng = np.random.default_rng(a.seed)
W, H = 720, 960
yy, xx = np.mgrid[0:H, 0:W].astype(np.float32)

def rgb(h): return np.array([int(h[i:i + 2], 16) for i in (0, 2, 4)], np.float32) / 255
def mix(a_, b_, t): return a_ + (b_ - a_) * np.clip(t, 0, 1)[..., None]
def smooth(lo, hi, v): t = np.clip((v - lo) / (hi - lo), 0, 1); return t * t * (3 - 2 * t)
def blur(layer, r):
    """gaussian blur of a single-channel float layer: one pass across, one pass down"""
    n = max(int(r * 3), 1)
    k = np.exp(-0.5 * (np.arange(-n, n + 1) / r) ** 2).astype(np.float32); k /= k.sum()
    out = layer.astype(np.float32)
    for axis in (0, 1):
        pad = np.pad(out, [(n, n) if ax == axis else (0, 0) for ax in (0, 1)], mode="edge")
        out = sum(w * (pad[i:i + out.shape[0]] if axis == 0 else pad[:, i:i + out.shape[1]]) for i, w in enumerate(k))
    return out
def noise(scale):
    small = rng.random((H // scale + 3, W // scale + 3)).astype(np.float32)
    return np.asarray(Image.fromarray(small, "F").resize((W, H), Image.BICUBIC), np.float32)
def stamp(layer, x, y, r, amount=1.0):
    """add a soft round dab to a single-channel layer"""
    x0, x1, y0, y1 = max(int(x - r - 2), 0), min(int(x + r + 3), W), max(int(y - r - 2), 0), min(int(y + r + 3), H)
    if x0 >= x1 or y0 >= y1: return
    d = np.hypot(xx[y0:y1, x0:x1] - x, yy[y0:y1, x0:x1] - y)
    layer[y0:y1, x0:x1] = np.maximum(layer[y0:y1, x0:x1], amount * np.clip(r + 0.5 - d, 0, 1))

ground_at = lambda x: 0.43 * H + 12 * np.sin(x * 0.012 + 1.0) + 6 * np.sin(x * 0.031)
gline = ground_at(xx)
below = smooth(-1.5, 1.5, yy - gline)                    # 0 in the air, 1 in the soil

# --- the air ---
sky = mix(rgb("0a1022"), rgb("1f4052"), (yy / gline) ** 1.6)
moon_x, moon_y = 0.71 * W, 0.12 * H
md = np.hypot(xx - moon_x, yy - moon_y)
sky += rgb("e9e2cf") * (np.exp(-(md / 120) ** 2) * 0.16)[..., None]
sky = mix(sky, rgb("efe8d2"), smooth(34, 31, md))
sky = mix(sky, rgb("cfc8b2"), smooth(9, 7, np.hypot(xx - moon_x + 9, yy - moon_y + 6)) * 0.5)   # one grey sea on the moon
stars = np.zeros((H, W), np.float32)
for _ in range(90):
    sx, sy = rng.uniform(0, W), rng.uniform(0, 0.36 * H)
    stamp(stars, sx, sy, rng.choice([0.6, 0.9, 1.3], p=[0.6, 0.3, 0.1]), rng.uniform(0.5, 1))
sky += rgb("e9e2cf") * (stars + blur(stars, 2.5) * 1.5)[..., None]
# far trees, then near ones
trees = np.zeros((H, W), np.float32)
far = np.zeros((H, W), np.float32)
for tx, tw in ((0.27, 9), (0.47, 12), (0.62, 8), (0.80, 10), (0.36, 6)):
    far = np.maximum(far, smooth(tw, tw - 2, np.abs(xx - tx * W - 6 * np.sin(yy * 0.01 + tx * 9))))
sky = mix(sky, rgb("10202e"), far * 0.85)
TRUNKS = ((0.10, 40), (0.555, 25), (0.93, 52))
for tx, tw in TRUNKS:
    flare = 1 + 1.6 * np.exp(-np.clip(gline - yy, 0, None) / 38)              # wider where it meets the ground
    trees = np.maximum(trees, smooth(tw * flare, tw * flare - 2.5, np.abs(xx - tx * W - 5 * np.sin(yy * 0.008 + tx * 5))))
bark = rgb("0b0e15") + (noise(5) - 0.5)[..., None] * 0.03 + rgb("1f4052") * (smooth(0.6, 1.0, trees) * 0)[..., None]
sky = mix(sky, bark, trees)
# branches leaving the top of the near trunks
br = np.zeros((H, W), np.float32)
for tx, tw in TRUNKS:
    for _ in range(4):
        x, y, ang, r = tx * W, rng.uniform(0.03, 0.25) * H, rng.choice([-1, 1]) * rng.uniform(0.5, 1.2), tw * 0.22
        while r > 1 and 0 < x < W and y > -10:
            stamp(br, x, y, r); x += np.sin(ang) * 2; y -= np.cos(ang) * 2; ang += rng.normal(0, 0.05); r *= 0.991
sky = mix(sky, rgb("0b0e15"), br)

# --- the soil ---
depth = np.clip((yy - gline) / (H - gline), 0, 1)
soil = mix(rgb("2e2016"), rgb("0e0a0b"), depth ** 0.8)
soil += ((noise(48) - 0.5) * 0.07 + (noise(9) - 0.5) * 0.035)[..., None] * rgb("ffd9a8")
band = np.sin(depth * 21 + noise(120) * 5)                                    # faint layers, the way soil settles
soil += (band * 0.012)[..., None]
for _ in range(26):                                                           # stones
    px, py, pr = rng.uniform(0, W), rng.uniform(0.5, 0.99) * H, rng.uniform(4, 13)
    d = np.hypot((xx - px) / 1.5, yy - py)
    soil = mix(soil, rgb("3b3632") + (py / H - 0.7) * -0.1, smooth(pr, pr - 2, d) * 0.8)
img = mix(sky, soil, below)
turf = smooth(9, 0, np.abs(yy - gline - 3))                                   # a dark green skin on top
img = mix(img, rgb("1c2e22"), turf * 0.9)

# --- roots ---
roots = np.zeros((H, W), np.float32)
root_tips = []
for tx, tw in TRUNKS:
    for k in range(6):
        x = tx * W + rng.uniform(-tw, tw); y = float(ground_at(x)) - 4
        ang = rng.uniform(-1.25, 1.25); r = tw * rng.uniform(0.22, 0.34)
        while r > 0.9 and y < H + 5 and -20 < x < W + 20:
            stamp(roots, x, y, r)
            x += np.sin(ang) * 1.6; y += np.cos(ang) * 1.6
            ang = np.clip(ang + rng.normal(0, 0.045), -1.45, 1.45); r *= 0.9962
            if r > 2.5 and rng.random() < 0.006:                              # a side root
                x2, y2, a2, r2 = x, y, ang + rng.choice([-1, 1]) * rng.uniform(0.5, 0.9), r * 0.6
                while r2 > 0.9 and y2 < H and 0 < x2 < W:
                    stamp(roots, x2, y2, r2); x2 += np.sin(a2) * 1.6; y2 += np.cos(a2) * 1.6
                    a2 = np.clip(a2 + rng.normal(0, 0.05), -1.5, 1.5); r2 *= 0.994
                root_tips.append((x2, y2))
        root_tips.append((x, y))
root_col = rgb("4a3020") + (noise(6) - 0.5)[..., None] * 0.05
img = mix(img, root_col * 0.55, blur(roots, 2.0) * below * 0.9)               # a dark edge first
img = mix(img, root_col, smooth(0.35, 0.9, blur(roots, 1.2)) * below)

# --- three morels ---
MORELS = ((0.335, 78), (0.43, 46), (0.715, 62))                               # where, how tall
bases = []
for mx, mh in MORELS:
    bx = mx * W; by = float(ground_at(bx)) + 3
    bases.append((bx, by))
    lean = (mx - 0.5) * 0.25
    sh, sw = mh * 0.42, mh * 0.17                                             # stem
    u = (by - yy) / sh
    cx = bx + lean * (by - yy)
    stem = smooth(sw * (1.05 - 0.2 * np.clip(u, 0, 1)), sw * (1.05 - 0.2 * np.clip(u, 0, 1)) - 2, np.abs(xx - cx)) * (u > -0.05) * (u < 1.1)
    shade = np.clip((xx - cx) / sw * 0.5 + 0.5, 0, 1)
    img = mix(img, mix(rgb("e4dbc0"), rgb("9a9079"), shade), stem)
    ch, cw = mh * 0.66, mh * 0.27                                             # cap: an egg, fatter low down
    cyc = by - sh - ch * 0.42
    v = (yy - cyc) / (ch * 0.55)
    capx = bx + lean * (by - cyc)
    width = cw * np.sqrt(np.clip(1 - v * v, 0, 1)) * (1 + 0.28 * v)
    cap = smooth(0, 2.5, width - np.abs(xx - capx)) * (np.abs(v) < 1)
    side = np.clip((xx - capx) / (cw * 1.2) * 0.5 + 0.5, 0, 1)
    capcol = mix(rgb("d9a651"), rgb("7d5526"), side * 0.9 + np.clip(v, -1, 1) * 0.12)
    cell = mh / 8.5                                                           # the pits: a honeycomb of dark hollows
    q = (yy - cyc) / (cell * 0.86); row = np.round(q)
    p = (xx - capx) / cell - 0.5 * (row % 2)
    pit = np.hypot((p - np.round(p)) * cell, (q - row) * cell * 0.86)
    capcol = mix(capcol, rgb("3d2912"), smooth(cell * 0.40, cell * 0.22, pit) * 0.85)
    img = mix(img, capcol, cap)
    img += rgb("e6c27a") * (np.exp(-(np.hypot(xx - capx, yy - cyc) / (mh * 0.9)) ** 2) * 0.05 * (1 - below))[..., None]

# --- the fungus: threads from under each morel, feeling their way down and towards the roots ---
threads = np.zeros((H, W), np.float32)
knots = np.zeros((H, W), np.float32)
def grow(x, y, ang, life, weight, depth_left):
    target = root_tips[rng.integers(len(root_tips))] if rng.random() < 0.7 else None
    while life > 0 and 2 < x < W - 2 and y < H - 2:
        if y > ground_at(x) + 2:
            stamp(threads, x, y, 0.75, weight)
        x += np.sin(ang) * 1.3; y += np.cos(ang) * 1.3; life -= 1
        ang += rng.normal(0, 0.16)
        if target is not None:                                                # lean a little towards a root
            want = np.arctan2(target[0] - x, target[1] - y)
            ang += 0.03 * np.arctan2(np.sin(want - ang), np.cos(want - ang))
            if np.hypot(target[0] - x, target[1] - y) < 6:
                stamp(knots, x, y, 2.2, 1.0); return
        ang = np.clip(ang, -2.0, 2.0)
        if y < ground_at(x) + 6: ang *= 0.9                                   # stay under the turf
        if depth_left and rng.random() < 0.0065:
            stamp(knots, x, y, 1.3, 0.7)
            grow(x, y, ang + rng.choice([-1, 1]) * rng.uniform(0.4, 1.1), life * rng.uniform(0.5, 0.9), weight * 0.78, depth_left - 1)
for bx, by in bases:
    for _ in range(6):
        grow(bx + rng.uniform(-7, 7), by + 4, rng.uniform(-1.5, 1.5), rng.integers(240, 560), 1.0, 3)
    stamp(knots, bx, by + 6, 4, 1.0)
glow = blur(threads, 0.6) * 0.9 + blur(threads, 5) * 0.55 + blur(knots, 1.2) * 1.0 + blur(knots, 9) * 1.1
img += rgb("a5e0d6") * (np.clip(glow, 0, 1.4) * below)[..., None] * 0.85
img = mix(img, rgb("eafaf2"), np.clip(blur(threads, 0.5) * 0.75 + blur(knots, 1.0), 0, 1) * below * 0.8)

img += rng.normal(0, 0.006, (H, W, 1)).astype(np.float32)
Image.fromarray((np.clip(img, 0, 1) * 255 + 0.5).astype(np.uint8)).save(a.out)
print("scene", W, "x", H, "threads cover %.1f%%" % (100 * (threads > 0.3).mean()))