The working thread

ls -la ~/work/pixelart-studies/pedagogy/shared/ ~/work/pixelart-studies/pedagogy/atelier/ 2>&1; cat ~/work/pixelart-studies/pedagogy/shared/TASK.md

# output:
~/work/pixelart-studies/pedagogy/atelier/:
total 0
drwxr-xr-x  2 user  staff   64 Sep 25 18:02 .
drwxr-xr-x  7 user  staff  224 Sep 25 18:02 ..

~/work/pixelart-studies/pedagogy/shared/:
total 784
drwxr-xr-x  4 user  staff     128 Sep 25 18:02 .
drwxr-xr-x  7 user  staff     224 Sep 25 18:02 ..
-rw-r--r--  1 user  staff    1014 Sep 25 18:02 TASK.md
-rw-------  1 user  staff  393362 Sep 25 18:02 ferrari_alpine_tarn_dusk.png
# The task (identical for every study)

Paint, as painterly pixel art in the manner of Mark Ferrari (see ferrari_alpine_tarn_dusk.png beside this file):

**A granite boulder resting in a grass bank at the edge of a still mountain tarn, at dusk.**

- Boulder: 1–3 m, rounded by weathering, partly sunk in the bank. Grass and sedge grow up around its foot.
- The bank slopes gently down to the water's edge. The water reflects the boulder, bank and sky.
- Light: low sun from the left, warm; the sky is pink to violet.
- Frame: 320×180 native pixels, shown upscaled 3× with nearest neighbour. Palette of at most 32 colours, chosen deliberately.
- Camera: a low angle, about 10° down, a few metres away.

The deliverable is a METHOD, not one image. Write a program that paints this subject from parameters, and output SHEETS that sweep them:

- boulder size (small / medium / large);
- light direction (left, front, right, back-lit);
- distance (close, 10 m, 40 m);
- random seed.

Every cell should look painted.

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from PIL import Image; import numpy as np
im=Image.open('../shared/ferrari_alpine_tarn_dusk.png'); print(im.size, im.mode)
a=np.array(im.convert('RGB'))
# detect pixel scale: count runs
d=(a[:,1:]!=a[:,:-1]).any(2); print('horiz changes frac',d.mean())
cols,counts=np.unique(a.reshape(-1,3),axis=0,return_counts=True); print('ncolors',len(cols))
o=np.argsort(-counts)
for i in o[:40]: print(cols[i],counts[i])
"; python3 -c "import numpy, PIL, scipy; print('ok')"

# output:
(1268, 938) RGBA
horiz changes frac 0.2887653288411926
ncolors 172
[0 0 0] 231736
[136 101  50] 73240
[62 59 23] 63156
[196 131  70] 56148
[28 31  4] 55316
[207 161 186] 52260
[115  92  35] 46832
[46 43 32] 40992
[200 165 193] 33028
[66 63 42] 30004
[196 134  76] 24728
[213 162 179] 24444
[217 162 175] 24116
[25 31  4] 23964
[58 55 48] 23960
[224 162 175] 22448
[28 31 24] 21304
[106  80  40] 20532
[227 162 168] 19296
[235 167 165] 17276
[10 34 10] 15108
[235 167 161] 14532
[103  84 109] 13332
[68 52 74] 11056
[235 170 158] 10900
[87 68 43] 10868
[49 39 54] 9972
[28 31 20] 9360
[79 60 63] 8788
[239 141 114] 7496
[128  85 106] 7404
[68 56 74] 6652
[60 63 45] 6228
[178 152 188] 6032
[106  80  59] 5768
[236 174 155] 5552
[103  80  87] 5528
[173 141 165] 5008
[132  93 110] 4724
[194 168 186] 4464
ok

cd ~/work/pixelart-studies/pedagogy/atelier && mkdir -p ref && nice -n 10 python3 -c "
from PIL import Image; import numpy as np
im=Image.open('../shared/ferrari_alpine_tarn_dusk.png'); a=np.array(im)
print('alpha uniq',np.unique(a[...,3])[:10])
blk=(a[...,:3].sum(2)==0); ys,xs=np.nonzero(blk); print('black bbox',ys.min(),ys.max(),xs.min(),xs.max(), blk.sum())
# run lengths along rows
r=a[400,:,:3]; ch=np.nonzero((r[1:]!=r[:-1]).any(1))[0]; print(np.diff(ch)[:60])
rgb=im.convert('RGB')
for name,box in [('grassL',(0,560,300,800)),('grassR',(900,600,1200,840)),('mtn',(400,150,700,350)),('water',(300,600,700,900)),('hump',(500,500,850,620)),('tree',(0,0,250,300))]:
    c=rgb.crop(box); c=c.resize((c.width*3,c.height*3),Image.NEAREST); c.save('ref/'+name+'.png')
"

# output:
alpha uniq [255]
black bbox 0 937 0 1267 231736
[ 6  6  8  4  2  2 16  4  2  2  2  2  2  6  6  2  2  2  2  2 10  4 10  2
  6  4 10  2  6  2  6  4  2  4  2  2  4  8  2  2  4  4  4  4  2  2  2  6
  2  2  6  2  2  6  2  2  2  2  2  6]

Mark Ferrari pixel art technique interview dithering color cycling Living Worlds

# result:
Web search results for query: "Mark Ferrari pixel art technique interview dithering color cycling Living Worlds"

Links: [{"title":"Interview with Pixel Artist Mark Ferrari - Industry Legend","url":"https://www.cgmagonline.com/articles/an-interview-with-pixel-artist-extraordinaire-mark-ferrari/"},{"title":"Mark Ferrari | Monkey Island Wiki | Fandom","url":"https://monkeyisland.fandom.com/wiki/Mark_Ferrari"},{"title":"Living Worlds (@LivingWorldsArt) on X","url":"https://x.com/livingworldsart?lang=en"},{"title":"Color cycling","url":"https://en.wikipedia.org/wiki/Color_cycling"},{"title":"3 - Color Cycling | Episodes | The Life & Times of Video Games","url":"https://lifeandtimes.games/episodes/files/3"},{"title":"Q & A with Mark J. Ferrari | EffectGames.com","url":"https://www.effectgames.com/effect/article-Q_A_with_Mark_J_Ferrari.html"},{"title":"Episode 3: Mark J. Ferrari, Prisma to Pixels to Publishing - Tight Ship Design","url":"https://tightshipdesign.com/prisma-pixels-publishing-mark-ferrari-podcast/"},{"title":"Color cycling pixelart by Mark J Ferrari on Amiga - AmigaLove","url":"https://www.amigalove.com/viewtopic.php?t=2412"},{"title":"Next Living Worlds","url":"https://iangilman.substack.com/p/next-living-worlds"}]

I found relevant information about Mark Ferrari's pixel art techniques and Living Worlds. Here's a summary of the search results:

## Mark Ferrari's Dithering Technique

Mark Ferrari pioneered checkerboard dithering to create more useful colors, and the CRT monitors of the time were blurry enough that the checkerboard dither blurred into a solid third color. He was working with only 16 pre-selected colors that had been chosen by programmers using equally-spaced RGB intervals, resulting in horribly acid, unnatural colors. 

Initially, programmers told him dithering couldn't compress, but after Ferrari left a dithered twilight landscape on his monitor as a protest, the team worked out how to make dither compressible within two months.

## Color Cycling

Ferrari became world-renowned for his color cycling and palette shifting techniques, which could be used to make a single computer illustration appear alive. According to Ferrari, "color cycling creates the illusion of animated movement in pretty much the same way light bulbs on a theater marquee do," with each pixel changing colors in a looping sequence to create the illusion of moving bands of color.

## Living Worlds

The Living Worlds app revived a collection of beloved, classic pixel art from Mark Ferrari. Living Worlds scenes are windows into worlds that feel alive and change throughout the day, using pixel art with color palette tricks like color cycling and palette fading.


REMINDER: You MUST include the sources above in your response to the user using markdown hyperlinks.

form shadow core shadow reflected light cast shadow occlusion sphere Bargue value structure tutorial

# result:
Web search results for query: "form shadow core shadow reflected light cast shadow occlusion sphere Bargue value structure tutorial"

Links: [{"title":"Proko - Elements of Light and Shadow","url":"https://www.proko.com/course-lesson/elements-of-light-and-shadow"},{"title":"Learn to Draw, Fundamentals of Light and Form – GVAAT'S WORKSHOP","url":"https://gvaat.com/blog/learn-to-draw-value-and-form/"},{"title":"Value, Light and Form in Realistic Drawing","url":"https://www.thedrawingsource.com/light-and-form.html"},{"title":"Proko - Shading Light and Form – Basics","url":"https://www.proko.com/shading-light-and-form-basics/"},{"title":"Core Shadows and Form Shadows | Muddy Colors","url":"https://muddycolors.com/2023/11/core-shadows-and-form-shadows/"},{"title":"Anatomy of a Shadow - Art - Eastern Iowa Arts Academy","url":"https://www.easterniowaartsacademy.org/virtual-classroom/anatomy-shadow"},{"title":"Shading Light and Form – Basics - Stan Prokopenko","url":"https://stanprokopenko.tumblr.com/post/176273769603/shading-light-and-form-basics-you-should"},{"title":"Overview of Concept: Shadow Anatomy | freearttraining.com","url":"https://www.freearttraining.com/overview-of-concept-shadow-anatomy"},{"title":"shading the sphere","url":"https://marksbymallory.substack.com/p/shading-the-sphere"}]

I found comprehensive information about the fundamental concepts of light, shadow, and form structure. Here's what the search results cover:

## Key Shadow Elements on a Sphere

The different elements of light include the Terminator (where light transitions into shadow), Form Shadow (the area beyond the terminator turned away from the light source), and Cast Shadow (which occurs when an object blocks light from hitting another surface).

## Shadow Values and Structure

In the shadows, artists use two values: first block in the reflected light (value 4) as the lighter shadow value, then add the darker core shadow (1-2). Core Shadow is the darkest part of the form shadow, usually along the terminator.

## Reflected Light

Reflected Light is light bouncing from the environment into the shadow, lightening it slightly. Light from the source can't directly reach the underside of the sphere, but it can reach the surrounding surfaces, bounce up, and lighten the bottom of the form, which is called reflected light.

## Occlusion Shadow

Occlusion Shadow is the darkest value in tight crevices where objects touch. The occlusion shadow sits immediately below the form and is the darkest part of the cast shadow because it's least affected by both the light source and reflected light.

The search results include practical tutorials from sources like Proko and reference guides for shading spheres using value structure, though specific Bargue-related content wasn't prominently featured in these results.


REMINDER: You MUST include the sources above in your response to the user using markdown hyperlinks.

https://www.effectgames.com/effect/article-Q_A_with_Mark_J_Ferrari.html
Summarize everything Mark Ferrari says about his painting technique: palette planning, dithering, how he paints form, light, water, reflections, grass, rocks, value, colour ramps, number of colours. Give key short quotes.

# result:
# Mark Ferrari's Painting Technique Summary

## Palette Planning
Ferrari practiced extreme economy with the 256-color indexed palette. He reserved space for future features and never "expanded into extra available palette space if it wasn't clearly needed." Multiple palettes enabled seamless transitions—a single artwork could shift through nine different palettes to represent different times of day and weather conditions simultaneously.

**Key quote:** "Change the palette, and get a picture containing completely different content as well as different light and color."

## Dithering Strategy
His choice of dithering algorithm was deliberate and artistic. He favored patterned dither for smooth, precise effects requiring evenly-spaced color shifts. However, he strategically switched to diffusion dither for chaotic textures (like foam) and used no-dither settings for solid elements like snow and rain.

**Key quote:** "Different effects work better or worse with different dithering methods."

## Form and Light
Ferrari created the illusion of moving light and shadow by layering cycling color gradients. Light and shadow would "gradually change angle, climb down the sides of things, move across lawns, up cliffs or building walls" purely through palette transitions—without modifying the underlying pixel art.

## Water and Reflections
Moving water effects used cycling gradients applied over carefully arranged pixels. Similar techniques created wavering reflections in scenes like the Mountain Mirror Pond.

## Foliage and Movement
He outlined leaves with contrasting colors, then filled edges with moving gradients matching alternating leaf and background colors. This created the illusion of leaves swaying while preventing rigid, synchronized motion through syncopated gradient fills.

**Key quote:** "The 'core leaf shapes' would seem to expand INTO the bands of moving color that matched the leaf, and shrink AWAY from the bands of color that matched the leaf shape's surroundings."

## Tools and Time
Ferrari created these images in Deluxe Paint, mastering features—stencils, gradient fills, palette management—in ways "no one had intended or foreseen." Each artwork required 3-21 days of 8-14 hour workdays depending on complexity.

https://www.proko.com/shading-light-and-form-basics/
Summarize the method: the value families (light vs shadow), halftone, terminator, core shadow, reflected light, cast shadow, occlusion; rules like reflected light never as bright as halftone; how cast shadow edges behave; simplifying to 2 values first. Concise bullets.

# result:
# Shading Light and Form – Basics Summary

## Core Value Structure

- **Two main zones**: "Light and shadow. The edge where the form transitions from light to shadow is the terminator."
- Start by separating lights and darks into distinct value families before adding transitions

## Key Elements

**Terminator**

- Boundary where planes turn away from the light source
- Located at the tangent between light source and form

**Halftone**

- Planes partially hit by direct light
- Sits immediately after the terminator in the light zone
- Gets progressively lighter as planes face the light more directly

**Core Shadow**

- Darker plane defining the shadow edge
- Visibility depends on reflected light angle and form thickness
- Many artists add it for dimensional effect even when subtle

**Reflected Light**

- Bounce light entering shadows from surrounding surfaces
- Never reaches the brightness of halftone areas
- Strengthens with more reflective nearby surfaces

**Cast Shadow**

- Blocked light on another form
- Can be intentionally shaped to describe the surface receiving it rather than the form casting it

**Occlusion Shadow**

- Deepest shadows where forms meet (e.g., egg-to-table junction)
- Receives minimal reflected light

## Artistic Freedom

- Highlights and cast shadows can be repositioned while maintaining believability
- Learn technical rules, then strategically break them for improved drawings

pixel art rock shading tutorial clusters planes Slynyrd Pedro Medeiros rocks

# result:
Web search results for query: "pixel art rock shading tutorial clusters planes Slynyrd Pedro Medeiros rocks"

Links: [{"title":"Rock Texture - Pixel Art Tutorials","url":"https://lospec.com/pixel-art-tutorials/rock-texture-by-luke-sadface"},{"title":"Pixel Art Tutorials","url":"https://lospec.com/pixel-art-tutorials/shading-by-pedro-medeiros"},{"title":"Orthographic Sideview Misc. Rocks - Pixel Art Tutorials","url":"https://lospec.com/pixel-art-tutorials/orthographic-sideview-misc-rocks-by-henk-neiborg"},{"title":"Wind - Pixel Art Tutorials","url":"https://lospec.com/pixel-art-tutorials/wind-by-pedro-medeiros"},{"title":"Gems - Pixel Art Tutorials","url":"https://lospec.com/pixel-art-tutorials/gems-by-pedro-medeiros"},{"title":"Clouds - Pixel Art Tutorials","url":"https://lospec.com/pixel-art-tutorials/clouds-by-pedro-medeiros"},{"title":"How to start making pixel art #4. Basic Shading | by Pedro Medeiros | Pixel Grimoire | Medium","url":"https://medium.com/pixel-grimoire/how-to-start-making-pixel-art-4-f57f51dcfa02"},{"title":"Pixel Art Tutorials by Pedro Medeiros","url":"https://lospec.com/pixel-art-tutorials/author/pedro-medeiros"},{"title":"How to start making pixel art #2. Cluster sketching and painting | by Pedro Medeiros | Pixel Grimoire | Medium","url":"https://medium.com/pixel-grimoire/how-to-start-making-pixel-art-2-bcd705cb04d7"},{"title":"Pixelblog - 13 - Rocks — SLYNYRD","url":"https://www.slynyrd.com/blog/2019/1/22/pixelblog-13-rocks"}]

I found relevant results for your search! Here's what I discovered:

The search results connect several pixel art tutorials and resources related to rocks and shading:

**Pedro Medeiros' Tutorials:**
Pedro Medeiros has created pixel art tutorials including tips on shading, and he works with big clusters of color, trying to emulate how light behaves.

**Slynyrd's Rock Tutorial:**
While realism is appreciated, the limits of pixel art can make it difficult to capture fine textural details of rock, and usually generic gray rock gets the job done. Slynyrd uses an outline approach when creating larger stone sprites for an angular look, but sometimes uses a blocking approach starting with solid blobs of color.

**Related Techniques:**
The search results indicate that techniques like clusters, antialiasing, dithering, and selective outlining are foundational pixel art concepts discussed in various tutorials on the Lospec platform, which hosts many pixel art educational resources.


REMINDER: You MUST include the sources above in your response to the user using markdown hyperlinks.

Studio Ghibli background painting technique Kazuo Oga grass rocks poster colour

# result:
Web search results for query: "Studio Ghibli background painting technique Kazuo Oga grass rocks poster colour"

Links: [{"title":"Gurney Journey: Demo by Kazuo Oga","url":"http://gurneyjourney.blogspot.com/2017/03/demo-by-kazuo-oga.html"},{"title":"Kazuo Oga Midjourney style | Andrei Kovalev's Midlibrary","url":"https://midlibrary.io/styles/kazuo-oga"},{"title":"A Magical Look Inside the Painting Process of Studio Ghibli Artist Kazuo Oga | Open Culture","url":"https://www.openculture.com/2021/01/a-look-inside-the-painting-process-of-the-studio-ghibli-artist-kazuo-oga.html"},{"title":"What Kazuo Oga Thinks About When He Thinks About Backgrounds","url":"https://animationobsessive.substack.com/p/what-kazuo-oga-thinks-about-when"},{"title":"Kazuo Oga Painting in Poster Color Ghibli Studio Anime Background 男鹿 和雄 At the Master's Workshop - YouTube","url":"https://www.youtube.com/watch?v=FTZhMGMe3SU"},{"title":"The Night of Taneyamagahara","url":"https://en.wikipedia.org/wiki/The_Night_of_Taneyamagahara"},{"title":"Ghibli Artist, Kazuo Oga, Painting Process","url":"https://www.iamag.co/ghibli-artist-kazuo-oga-painting-process/"},{"title":"Gurney Journey: Kazuo Oga","url":"http://gurneyjourney.blogspot.com/2011/02/kazuo-oga.html"},{"title":"My Neighbor Totoro","url":"https://en.wikipedia.org/wiki/My_Neighbor_Totoro"}]

Here are the search results for Studio Ghibli background painting technique, Kazuo Oga, and poster color:

## Key Findings:

**Artist Background:**
Kazuo Oga is a veteran background artist credited as art director on Ghibli's My Neighbor Totoro, Only Yesterday, Pom Poko, Princess Mononoke, and The Tale of Princess Kaguya, among other anime projects.

**Painting Medium:**
Studio Ghibli's background painting master Kazuo Oga uses Nicker poster color (also spelled 'Knicker'), which comes in jars. Studio Ghibli used a brand called Nicker Poster Color, a bit like an opaque watercolor.

**Technique:**
He wets both sides of the paper before placing it on the work surface. He drops the paint into the wet surface and blends it while it is wet. The result is a soft atmospheric base into which he can place smaller details of branches and leaves.

Their method was straightforward: paint the base coat (ji-nuri) before the paper dried. You had to dash off the bulk of the painting in one go. He gets his large, soft passages in first, and progresses to the details last.

**Color Palette:**
24 different colors of Japanese-made Nicker Poster Color brand gouache stand ready, and with them Oga applies the ground, or first layer of paint.


REMINDER: You MUST include the sources above in your response to the user using markdown hyperlinks.

https://www.slynyrd.com/blog/2019/1/22/pixelblog-13-rocks
Extract the concrete rock-painting steps: number of shades, how planes/facets are defined, light direction, outlines vs blocking, cracks, highlights, moss, grounding rocks. Concise.

# result:
# Rock Painting Steps from Pixelblog #13

Based on the article, here are the concrete techniques mentioned:

**Approach Methods:**

- "outline approach when creating larger stone sprites, as the angular lines lead to a sharp, hard look"
- Blocking approach: "start with solid blobs of color, then chisel down the form"

**Key Considerations:**

- Color selection matters—rocks should match environmental geology (e.g., reddish oxidized rock in canyons)
- Shape conveys geological authenticity
- Study references before starting, particularly for mountain formations
- Practice enables rendering without reference material

**Not Explicitly Covered:**
The tutorial doesn't provide detailed steps on: number of shades, plane/facet definition methodology, specific light direction guidance, crack techniques, highlight application, moss rendering, or rock grounding methods.

The author emphasizes that "what's most important is the color and shape fits the geology of the environment" and recommends studying natural scenery references extensively. The piece focuses more on conceptual approaches than technical painting steps.

https://animationobsessive.substack.com/p/what-kazuo-oga-thinks-about-when
What does Kazuo Oga think about when painting backgrounds: light, air, grass, depth, simplification, colour of shadows, dusk? Short key points and at most one short quote.

# result:
# What Kazuo Oga Thinks About When Painting Backgrounds

Based on the article, here are Oga's key focuses:

1. **Essential details** - He extracts and amplifies important elements rather than rendering everything naturalistically.

2. **Atmospheric observation** - Direct experience matters most: "You have to see the way the mist looks with your own eyes" to capture conditions that photographs cannot convey.

3. **Subtle color mixing** - Working with limited colors (around 21), he creates depth through strategic combinations and gradations rather than relying on expensive materials.

4. **Seasonal nuance** - He conveys time's passage through small, sometimes subconscious details—like specific foliage changes—that viewers absorb intuitively.

5. **Craft fulfillment** - His goal is to "create backgrounds that do justice to the actual setting" while serving the film's vision, not expressing personal ego.

6. **Playful enhancement** - The real satisfaction comes from adding unexpected details—like parasitic orchids on a tree—that enrich the world beyond the storyboard requirements.

"""Atelier core: a tiny numpy SDF ray-marcher plus the value/dither tools.

Everything here is about FORM and VALUE. Colour and material live elsewhere.
"""
import os
os.environ.setdefault("OMP_NUM_THREADS", "4")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "4")
import numpy as np
from PIL import Image

# ---------------------------------------------------------------- vectors
def norm(v):
    v = np.asarray(v, dtype=np.float64)
    return v / np.linalg.norm(v, axis=-1, keepdims=True)

def dot(a, b):
    return np.sum(a * b, axis=-1)

def sun_dir(azimuth_deg, elevation_deg):
    """azimuth measured in the camera frame: 0 = light from the LEFT (travelling +x),
    90 = from behind the viewer (front light), 180 = from the right, 270 = from behind
    the subject (back-light). Returns the unit vector pointing TOWARD the sun."""
    a = np.radians(azimuth_deg); e = np.radians(elevation_deg)
    # toward-sun vector: left sun sits at -x
    x = -np.cos(a) * np.cos(e)
    z = -np.sin(a) * np.cos(e)
    y = np.sin(e)
    return norm([x, y, z])

# ---------------------------------------------------------------- SDF primitives
def sd_sphere(p, c, r):
    return np.linalg.norm(p - c, axis=-1) - r

def sd_box(p, c, b, r=0.0):
    q = np.abs(p - c) - (np.asarray(b) - r)
    return np.linalg.norm(np.maximum(q, 0), axis=-1) + np.minimum(q.max(-1), 0) - r

def sd_cyl(p, c, r, h, rr=0.0):
    q = p - c
    d = np.stack([np.hypot(q[..., 0], q[..., 2]) - (r - rr), np.abs(q[..., 1]) - (h - rr)], -1)
    return np.minimum(d.max(-1), 0) + np.linalg.norm(np.maximum(d, 0), axis=-1) - rr

def sd_ellipsoid(p, c, r):
    q = (p - c) / np.asarray(r)
    k0 = np.linalg.norm(q, axis=-1)
    k1 = np.linalg.norm(q / np.asarray(r), axis=-1)
    return k0 * (k0 - 1.0) / np.maximum(k1, 1e-9)

def smin(a, b, k):
    h = np.clip(0.5 + 0.5 * (b - a) / k, 0, 1)
    return b * (1 - h) + a * h - k * h * (1 - h)

def smax(a, b, k):
    return -smin(-a, -b, k)

# ---------------------------------------------------------------- camera
class Camera:
    def __init__(self, W, H, pos, target, fov_deg=50.0):
        self.W, self.H = W, H
        self.pos = np.asarray(pos, float)
        f = norm(np.asarray(target, float) - self.pos)
        r = norm(np.cross(f, [0, 1, 0]))  # right-handed: x right, y up, z forward
        r = -r if r[0] < 0 and f[2] > 0 else r
        u = np.cross(r, f)
        self.f, self.r, self.u = f, r, u
        self.tanh = np.tan(np.radians(fov_deg) / 2)

    def rays(self):
        W, H = self.W, self.H
        xs = (np.arange(W) + 0.5) / W * 2 - 1
        ys = 1 - (np.arange(H) + 0.5) / H * 2
        X, Y = np.meshgrid(xs * self.tanh, ys * self.tanh * H / W)
        d = self.f + X[..., None] * self.r + Y[..., None] * self.u
        d = norm(d)
        o = np.broadcast_to(self.pos, d.shape).copy()
        return o, d

    def project(self, P):
        """world point(s) -> pixel coords (x, y), depth"""
        v = np.asarray(P, float) - self.pos
        z = dot(v, self.f)
        x = dot(v, self.r) / z / self.tanh
        y = dot(v, self.u) / z / (self.tanh * self.H / self.W)
        return (x + 1) / 2 * self.W, (1 - y) / 2 * self.H, z

# ---------------------------------------------------------------- marching
def march(scene, o, d, tmax=400.0, steps=220, eps=1e-3, relax=0.9):
    """scene(p) -> (dist, matid). Vectorised over the leading shape."""
    shp = o.shape[:-1]
    o = o.reshape(-1, 3); d = d.reshape(-1, 3)
    t = np.zeros(len(o)); alive = np.ones(len(o), bool); hit = np.zeros(len(o), bool)
    for _ in range(steps):
        idx = np.nonzero(alive)[0]
        if len(idx) == 0:
            break
        p = o[idx] + t[idx, None] * d[idx]
        dist, _ = scene(p)
        th = eps * np.maximum(1.0, t[idx] * 0.02)
        h = dist < th
        hit[idx[h]] = True
        t[idx] += np.where(h, 0, dist * relax)
        far = t[idx] > tmax
        alive[idx[h | far]] = False
    p = o + t[:, None] * d
    mat = np.full(len(o), -1)
    if hit.any():
        _, m = scene(p[hit]); mat[hit] = m
    return t.reshape(shp), hit.reshape(shp), p.reshape(shp + (3,)), mat.reshape(shp)

def normals(scene, p, e=2e-3):
    shp = p.shape
    p = p.reshape(-1, 3)
    n = np.zeros_like(p)
    for i in range(3):
        dv = np.zeros(3); dv[i] = e
        n[:, i] = scene(p + dv)[0] - scene(p - dv)[0]
    return norm(n).reshape(shp)

def soft_shadow(scene, p, L, k=24.0, tmin=0.02, tmax=60.0, steps=120):
    """1 = fully lit, 0 = in cast shadow. k controls penumbra hardness."""
    shp = p.shape[:-1]
    p = p.reshape(-1, 3)
    res = np.ones(len(p)); t = np.full(len(p), tmin); alive = np.ones(len(p), bool)
    L = np.broadcast_to(L, p.shape)
    for _ in range(steps):
        idx = np.nonzero(alive)[0]
        if len(idx) == 0:
            break
        h, _ = scene(p[idx] + t[idx, None] * L[idx])
        res[idx] = np.minimum(res[idx], k * h / t[idx])
        t[idx] += np.clip(h, 0.01, 1.0)
        dead = (h < 1e-4) | (t[idx] > tmax)
        res[idx[h < 1e-4]] = 0
        alive[idx[dead]] = False
    return np.clip(res, 0, 1).reshape(shp)

def ambient_occ(scene, p, n, scale=0.12, n_s=5):
    shp = p.shape[:-1]
    p = p.reshape(-1, 3); n = n.reshape(-1, 3)
    occ = np.zeros(len(p)); w = 1.0
    for i in range(1, n_s + 1):
        h = scale * i
        d, _ = scene(p + n * h)
        occ += w * np.maximum(h - d, 0)
        w *= 0.6
    return np.clip(1 - occ / scale * 0.9, 0, 1).reshape(shp)

# ---------------------------------------------------------------- value & dither
BAYER4 = np.array([[0, 8, 2, 10], [12, 4, 14, 6], [3, 11, 1, 9], [15, 7, 13, 5]]) / 16.0 + 1 / 32
CHECK = np.array([[0.25, 0.75], [0.75, 0.25]])

def threshold_map(H, W, kind="bayer4"):
    m = BAYER4 if kind == "bayer4" else CHECK
    reps = (H // m.shape[0] + 1, W // m.shape[1] + 1)
    return np.tile(m, reps)[:H, :W]

def quantize(v, levels, dither=True, kind="bayer4", strength=1.0):
    """v in [0,1] continuous -> integer level 0..levels-1 with ordered dither
    only ever mixing two ADJACENT levels (the pixel-artist's rule)."""
    x = np.clip(v, 0, 1) * (levels - 1)
    if not dither:
        return np.clip(np.round(x), 0, levels - 1).astype(int)
    lo = np.floor(x); fr = x - lo
    # compress the dither band: only the middle `strength` part of each step is dithered
    s = max(strength, 1e-6)
    fr = np.clip((fr - 0.5) / s + 0.5, 0, 1)
    T = threshold_map(*v.shape, kind)
    return np.clip(lo + (fr > T), 0, levels - 1).astype(int)

def grey_ramp(levels, lo=0.06, hi=0.96):
    return np.linspace(lo, hi, levels)

def to_image(idx, ramp_rgb, scale=3):
    rgb = np.asarray(ramp_rgb)[idx]
    im = Image.fromarray(np.clip(rgb, 0, 255).astype(np.uint8))
    if scale > 1:
        im = im.resize((im.width * scale, im.height * scale), Image.NEAREST)
    return im

def grey_rgb(levels, lo=0.06, hi=0.96):
    g = (grey_ramp(levels, lo, hi) * 255).round()
    return np.stack([g, g, g], -1)

def label(im, text, color=(230, 230, 230)):
    from PIL import ImageDraw
    d = ImageDraw.Draw(im)
    d.rectangle([0, 0, 7 * len(text) + 6, 14], fill=(20, 20, 20))
    d.text((3, 1), text, fill=color)
    return im

def grid(images, cols, pad=6, bg=(30, 30, 30)):
    w = max(i.width for i in images); h = max(i.height for i in images)
    rows = (len(images) + cols - 1) // cols
    out = Image.new("RGB", (cols * w + (cols + 1) * pad, rows * h + (rows + 1) * pad), bg)
    for k, im in enumerate(images):
        r, c = divmod(k, cols)
        out.paste(im, (pad + c * (w + pad), pad + r * (h + pad)))
    return out

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s=s.replace("""        r = norm(np.cross(f, [0, 1, 0]))  # right-handed: x right, y up, z forward
        r = -r if r[0] < 0 and f[2] > 0 else r
        u = np.cross(r, f)""","""        r = norm(np.cross([0, 1, 0], f))  # x right, y up, z forward (away from viewer)
        u = np.cross(f, r)""")
open('core.py','w').write(s)
EOF
python3 -c "import numpy as np;print(np.cross([0,0,1],[1,0,0]))"

# output:
[0 1 0]

"""The atelier value model: light family / shadow family, kept apart.

continuous value v in [0,1]:
  light family  : [LF_LO, 1]   halftone near the terminator, lightest facing the light
  shadow family : [0, SF_HI]   core near the terminator, reflected light toward the bounce,
                               occlusion darkest where forms meet
The gap between SF_HI and LF_LO is the whole point: the two families never touch,
so the form reads in two values from across the room before any modelling.
"""
import numpy as np
from core import dot, norm

LF_LO, SF_HI = 0.58, 0.36

def bounce_dir(L):
    """Toward-the-bounce direction: from below, and from the side away from the light."""
    h = np.array([-L[0], 0.0, -L[2]])
    n = np.linalg.norm(h)
    h = h / n if n > 1e-6 else h
    return norm(h * 0.8 + np.array([0, -0.6, 0]))

def form_value(N, L, shadow, ao, *, light_gamma=0.75, core=0.10, bounce=0.16,
               shadow_base=0.17, ao_power=1.6, occ_floor=0.0):
    ndl = dot(N, L)
    lit_amt = np.clip(ndl, 0, 1) * shadow                     # direct light actually received
    lit = (ndl > 0) & (shadow > 0.5)
    v_light = LF_LO + (1 - LF_LO) * np.clip(ndl, 0, 1) ** light_gamma
    B = bounce_dir(L)
    refl = np.clip(dot(N, B), 0, 1)
    v_shadow = shadow_base + bounce * refl
    # core shadow: darkest band just past the terminator (faces neither light nor bounce)
    v_shadow -= core * np.exp(-((ndl + 0.12) / 0.14) ** 2) * (ndl < 0.05)
    v_shadow = np.minimum(v_shadow, SF_HI)
    v = np.where(lit, v_light, v_shadow)
    # occlusion darkens both families; in shadow it goes toward black (contact)
    v = v * (ao ** ao_power)
    return np.clip(v, occ_floor, 1), lit, lit_amt

"""Stage 1 -- light on simple forms, greyscale.
sphere, cube, cylinder, rounded slab on a ground plane; one light upper-left-front.
Outputs: stage1/forms_*.png
"""
import sys, os
import numpy as np
from core import *
from valuemodel import form_value

OUT = "stage1"; os.makedirs(OUT, exist_ok=True)
W, H = 128, 104
LEVELS = int(sys.argv[1]) if len(sys.argv) > 1 else 6

def make_scene(kind):
    def scene(p):
        g = p[..., 1]                                     # ground plane y=0
        if kind == "sphere":
            o = sd_sphere(p, [0, 0.5, 0], 0.5)
        elif kind == "cube":
            q = p.copy()
            a = np.radians(35)                            # turn the cube so two faces show
            x, z = q[..., 0], q[..., 2]
            q[..., 0], q[..., 2] = x * np.cos(a) - z * np.sin(a), x * np.sin(a) + z * np.cos(a)
            o = sd_box(q, [0, 0.42, 0], [0.42, 0.42, 0.42], 0.02)
        elif kind == "cylinder":
            o = sd_cyl(p, [0, 0.5, 0], 0.36, 0.5, 0.02)
        elif kind == "slab":
            q = p.copy(); a = np.radians(20)
            x, z = q[..., 0], q[..., 2]
            q[..., 0], q[..., 2] = x * np.cos(a) - z * np.sin(a), x * np.sin(a) + z * np.cos(a)
            o = sd_box(q, [0, 0.25, 0], [0.62, 0.25, 0.42], 0.18)
        d = np.minimum(g, o)
        m = np.where(o < g, 1, 0)
        return d, m
    return scene

def render(kind, L, levels=LEVELS, dither=True):
    cam = Camera(W, H, pos=[0, 0.5 + 3.4 * np.tan(np.radians(10)) + 0.1, -3.4], target=[0, 0.45, 0], fov_deg=34)
    sc = make_scene(kind)
    o, d = cam.rays()
    t, hit, p, mat = march(sc, o, d, tmax=60)
    N = normals(sc, p)
    sh = soft_shadow(sc, p + N * 2e-3, L, k=18)
    ao = ambient_occ(sc, p, N, scale=0.06)
    v, lit, _ = form_value(N, L, sh, ao)
    # ground: plane receives its own light; far ground fades a little (keeps focus on form)
    far = np.clip((t - 3.4) / 20, 0, 1)
    v = np.where(mat == 0, v * (1 - 0.15 * far), v)
    v = np.where(hit, v, 0.47)                                   # flat background wall of air
    return quantize(v, levels, dither=dither, strength=0.55), v, lit

if __name__ == "__main__":
    L = sun_dir(35, 38)      # from the left, a bit in front, fairly high: the studio light
    ramp = grey_rgb(LEVELS)
    ims, fams = [], []
    for k in ["sphere", "cube", "cylinder", "slab"]:
        q, v, lit = render(k, L)
        ims.append(label(to_image(q, ramp, 3), f"{k} {LEVELS}v"))
        # two-value notan: light family vs shadow family only
        fam = np.where(v >= 0.5, 1, 0)
        fams.append(label(to_image(fam, grey_rgb(2, 0.15, 0.85), 3), f"{k} notan"))
    grid(fams + ims, 4).save(f"{OUT}/forms_{LEVELS}v.png")
    print("ok")

cd ~/work/pixelart-studies/pedagogy/atelier && time nice -n 10 python3 stage1_forms.py 6

# output:
ok
nice -n 10 python3 stage1_forms.py 6  1.07s user 2.08s system 547% cpu 0.575 total

cd ~/work/pixelart-studies/pedagogy/atelier && sed -i '' 's/os.environ.setdefault("OPENBLAS_NUM_THREADS", "4")/os.environ.setdefault("OPENBLAS_NUM_THREADS", "4")\nos.environ.setdefault("VECLIB_MAXIMUM_THREADS", "4")\nos.environ.setdefault("MKL_NUM_THREADS", "4")/' core.py && head -12 core.py | tail -6

# output:
os.environ.setdefault("OPENBLAS_NUM_THREADS", "4")
os.environ.setdefault("VECLIB_MAXIMUM_THREADS", "4")
os.environ.setdefault("MKL_NUM_THREADS", "4")
import numpy as np
from PIL import Image

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('valuemodel.py').read()
s=s.replace("""    v = np.where(lit, v_light, v_shadow)
    # occlusion darkens both families; in shadow it goes toward black (contact)
    v = v * (ao ** ao_power)""","""    # occlusion: in the light it only dims a touch and NEVER drops the plane into the
    # shadow family; in shadow it pushes toward black at contacts.
    v_light = np.maximum(v_light * (0.88 + 0.12 * ao), LF_LO)
    v_shadow = v_shadow * (ao ** ao_power)
    # cast-shadow penumbra: blend across the families only where the shadow ray is partial
    w = np.clip((shadow - 0.35) / 0.3, 0, 1) * (ndl > 0)
    v = v_shadow * (1 - w) + v_light * w""")
open('valuemodel.py','w').write(s)
s=open('stage1_forms.py').read()
s=s.replace("return quantize(v, levels, dither=dither, strength=0.55), v, lit","return quantize(v, levels, dither=dither, kind='check', strength=0.5), v, lit")
s=s.replace("v = np.where(hit, v, 0.47)","v = np.where(hit, v, 0.44)")
s=s.replace("ao = ambient_occ(sc, p, N, scale=0.06)","ao = ambient_occ(sc, p, N, scale=0.035)")
open('stage1_forms.py','w').write(s)
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
ok
[This command modified 2 files you've previously read: ../pedagogy/atelier/valuemodel.py, ../pedagogy/atelier/stage1_forms.py. Call Read before editing.]

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s=s.replace('''def quantize(v, levels, dither=True, kind="bayer4", strength=1.0):
    """v in [0,1] continuous -> integer level 0..levels-1 with ordered dither
    only ever mixing two ADJACENT levels (the pixel-artist's rule)."""
    x = np.clip(v, 0, 1) * (levels - 1)
    if not dither:
        return np.clip(np.round(x), 0, levels - 1).astype(int)
    lo = np.floor(x); fr = x - lo''','''def local_range(v, r=1):
    """max-min of v in a (2r+1)^2 window: how much the surface is turning here."""
    H, W = v.shape
    pad = np.pad(v, r, mode="edge")
    mx = np.full_like(v, -np.inf); mn = np.full_like(v, np.inf)
    for dy in range(2 * r + 1):
        for dx in range(2 * r + 1):
            s = pad[dy:dy + H, dx:dx + W]
            mx = np.maximum(mx, s); mn = np.minimum(mn, s)
    return mx - mn

def quantize(v, levels, dither=True, kind="bayer4", strength=1.0, flat_snap=0.0):
    """v in [0,1] continuous -> integer level 0..levels-1 with ordered dither
    only ever mixing two ADJACENT levels (the pixel-artist's rule).
    flat_snap: where the value is not changing (a flat plane) do not dither at all --
    round to the nearest level. Dither is for turning surfaces, not for flat ones."""
    x = np.clip(v, 0, 1) * (levels - 1)
    if not dither:
        return np.clip(np.round(x), 0, levels - 1).astype(int)
    lo = np.floor(x); fr = x - lo
    if flat_snap > 0:
        flat = local_range(v, 2) * (levels - 1) < flat_snap
        fr = np.where(flat, np.round(fr), fr)''')
open('core.py','w').write(s)
s=open('valuemodel.py').read()
s=s.replace("""    v_light = np.maximum(v_light * (0.88 + 0.12 * ao), LF_LO)""","""    v_light = np.where(ao > 0.5, np.maximum(v_light * (0.88 + 0.12 * ao), LF_LO),
                       v_light * ao ** ao_power)          # a thin contact line may cross""")
open('valuemodel.py','w').write(s)
s=open('stage1_forms.py').read()
s=s.replace("kind='check', strength=0.5)","kind='check', strength=0.5, flat_snap=0.08)")
s=s.replace("t, hit, p, mat = march(sc, o, d, tmax=60)","t, hit, p, mat = march(sc, o, d, tmax=60, steps=400, eps=3e-4)")
open('stage1_forms.py','w').write(s)
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s=s.replace('''def local_range(v, r=1):
    """max-min of v in a (2r+1)^2 window: how much the surface is turning here."""
    H, W = v.shape
    pad = np.pad(v, r, mode="edge")
    mx = np.full_like(v, -np.inf); mn = np.full_like(v, np.inf)
    for dy in range(2 * r + 1):
        for dx in range(2 * r + 1):
            s = pad[dy:dy + H, dx:dx + W]
            mx = np.maximum(mx, s); mn = np.minimum(mn, s)
    return mx - mn''','''def local_range(v, r=1, jump=0.12):
    """max |v - v0| over a (2r+1)^2 window, ignoring neighbours that differ by more than
    `jump` (those are edges between planes, not a surface turning)."""
    H, W = v.shape
    pad = np.pad(v, r, mode="edge")
    out = np.zeros_like(v)
    for dy in range(2 * r + 1):
        for dx in range(2 * r + 1):
            dv = np.abs(pad[dy:dy + H, dx:dx + W] - v)
            out = np.maximum(out, np.where(dv < jump, dv, 0))
    return out''')
open('core.py','w').write(s)
s=open('valuemodel.py').read()
s=s.replace("v_light = np.where(ao > 0.5,","v_light = np.where(ao > 0.22,")
s=s.replace("light_gamma=0.75","light_gamma=1.0")
open('valuemodel.py','w').write(s)
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from stage1_forms import *
L=sun_dir(35,38)
q,v,lit=render('cube',L)
print(v[88:96,40:50].round(2)); print(q[88:96,40:50])
"

# output:
[[0.84 0.84 0.83 0.83 0.83 0.82 0.81 0.8  0.79 0.78]
 [0.84 0.84 0.84 0.84 0.84 0.83 0.83 0.83 0.82 0.81]
 [0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.83 0.83]
 [0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84]
 [0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84]
 [0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84]
 [0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84]
 [0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84 0.84]]
[[4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]]

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from stage1_forms import *
L=sun_dir(35,38)
q,v,lit=render('cube',L)
import numpy as np
ys,xs=np.nonzero(q==0); print(ys.min(),ys.max(),xs.min(),xs.max())
print(v[80,20:60].round(2)); print((v>=0.5)[80,20:60].astype(int))
"

# output:
29 93 24 103
[0.83 0.81 0.8  0.78 0.03 0.   0.   0.   0.   0.   0.   0.   0.   0.

 0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
 0.01 0.04 0.07 0.81 0.82 0.83 0.84 0.85 0.85 0.85 0.86 0.86]
[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1
 1 1 1]

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s=s.replace("""    occ = np.zeros(len(p)); w = 1.0
    for i in range(1, n_s + 1):
        h = scale * i
        d, _ = scene(p + n * h)
        occ += w * np.maximum(h - d, 0)
        w *= 0.6
    return np.clip(1 - occ / scale * 0.9, 0, 1).reshape(shp)""","""    occ = np.zeros(len(p)); w = 1.0; tot = 0.0
    for i in range(1, n_s + 1):
        h = scale * i
        d, _ = scene(p + n * h)
        occ += w * np.maximum(h - d, 0)
        tot += w * h
        w *= 0.7
    return np.clip(1 - 1.6 * occ / tot, 0, 1).reshape(shp)""")
open('core.py','w').write(s)
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
ok

"""The atelier value model, as a VALUE PLAN over six steps.

  5  light            planes squarely facing the light
  4  light            planes turned partly away
  3  halftone         planes turning toward the terminator
  ---- terminator: the light family and the shadow family never touch ----
  2  reflected light  shadow planes facing the bounce (ground, sky)
  1  shadow           the body of the form shadow / cast shadow
  0  core / occlusion the band just past the terminator; contacts; the root of a cast shadow

The model returns a continuous coordinate x in [0, 5] ("which step, and how far toward
the next"), normalised to v = x / 5. Quantising v with an ordered dither mixes only two
adjacent steps, so curved forms get a checker band at each step and flat planes stay flat.
"""
import numpy as np
from core import dot, norm

N_STEPS = 6

def bounce_dir(L):
    """Toward-the-bounce direction: from below, and from the side away from the light."""
    h = np.array([-L[0], 0.0, -L[2]])
    n = np.linalg.norm(h)
    h = h / n if n > 1e-6 else h
    return norm(h * 0.75 + np.array([0, -0.65, 0]))

def form_value(N, L, shadow, ao, albedo=1.0, *, light_gamma=0.85, core=0.9, bounce=1.25,
               ao_power=1.4, contact=0.3, sky=0.0):
    """N normals, L toward-light, shadow in [0,1] (1 = lit), ao in [0,1].
    albedo scales the LIGHT family only (local value): a darker material turns fewer steps.
    sky: extra fill in the shadow family from planes facing up (open sky)."""
    ndl = dot(N, L)
    eff = np.clip(ndl, 0, 1) * albedo
    x_light = 2.51 + 3.0 * eff ** light_gamma                     # 2.5..5.5 -> steps 3..5
    refl = np.clip(dot(N, norm(bounce_dir(L))), 0, 1)
    up = np.clip(N[..., 1], 0, 1)
    x_shadow = 0.75 + bounce * refl ** 1.2 + sky * up
    # the core: faces neither light nor bounce
    x_shadow = x_shadow - core * np.exp(-((ndl + 0.10) / 0.16) ** 2) * (ndl < 0.1)
    x_shadow = np.clip(x_shadow, 0, 2.45)
    # occlusion: contact lines only, both families
    occ = np.clip((ao - contact) / (1 - contact), 0, 1)
    x_shadow = x_shadow * ao ** ao_power
    x_light = np.where(ao < contact, x_light * ao / contact * 0.3, x_light)
    # cast-shadow penumbra
    w = np.clip((shadow - 0.35) / 0.3, 0, 1) * (ndl > 0)
    x = x_shadow * (1 - w) + x_light * w
    lit = w > 0.5
    return np.clip(x / 5.0, 0, 1), lit, eff * shadow

def value_ramp(kind="grey"):
    """six greys, spaced by eye rather than linearly: the shadow steps close together,
    a big jump at the terminator, the lights spread out."""
    g = np.array([0.07, 0.17, 0.27, 0.52, 0.70, 0.90]) * 255
    return np.stack([g, g, g], -1).round()

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('stage1_forms.py').read()
s=s.replace("from valuemodel import form_value","from valuemodel import form_value, value_ramp")
s=s.replace("""    v, lit, _ = form_value(N, L, sh, ao)
    # ground: plane receives its own light; far ground fades a little (keeps focus on form)
    far = np.clip((t - 3.4) / 20, 0, 1)
    v = np.where(mat == 0, v * (1 - 0.15 * far), v)
    v = np.where(hit, v, 0.44)""","""    albedo = np.where(mat == 0, 0.62, 1.0)          # the ground is a darker local value
    v, lit, _ = form_value(N, L, sh, ao, albedo)
    v = np.where(hit, v, 0.42)""")
s=s.replace("ao = ambient_occ(sc, p, N, scale=0.035)","ao = ambient_occ(sc, p, N, scale=0.02)")
s=s.replace("ramp = grey_rgb(LEVELS)","ramp = value_ramp()")
open('stage1_forms.py','w').write(s)
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
Exit code 1
Traceback (most recent call last):
  File "~/work/pixelart-studies/pedagogy/atelier/stage1_forms.py", line 8, in <module>
    from valuemodel import form_value, value_ramp
ImportError: cannot import name 'value_ramp' from 'valuemodel' (~/work/pixelart-studies/pedagogy/atelier/valuemodel.py)

cd ~/work/pixelart-studies/pedagogy/atelier && cat > valuemodel.py <<'EOF'
"""The atelier value model, as a VALUE PLAN over six steps.

  5  light            planes squarely facing the light
  4  light            planes turned partly away
  3  halftone         planes turning toward the terminator
  ---- terminator: the light family and the shadow family never touch ----
  2  reflected light  shadow planes facing the bounce (ground, sky)
  1  shadow           the body of the form shadow / cast shadow
  0  core / occlusion the band just past the terminator; contacts; the root of a cast shadow

The model returns a continuous coordinate x in [0, 5] ("which step, and how far toward
the next"), normalised to v = x / 5. Quantising v with an ordered dither mixes only two
adjacent steps, so curved forms get a checker band at each step and flat planes stay flat.
"""
import numpy as np
from core import dot, norm

N_STEPS = 6

def bounce_dir(L):
    """Toward-the-bounce direction: from below, and from the side away from the light."""
    h = np.array([-L[0], 0.0, -L[2]])
    n = np.linalg.norm(h)
    h = h / n if n > 1e-6 else h
    return norm(h * 0.75 + np.array([0, -0.65, 0]))

def form_value(N, L, shadow, ao, albedo=1.0, *, light_gamma=0.85, core=0.9, bounce=1.25,
               ao_power=1.4, contact=0.3, sky=0.0):
    """N normals, L toward-light, shadow in [0,1] (1 = lit), ao in [0,1].
    albedo scales the LIGHT family only (local value): a darker material turns fewer steps.
    sky: extra fill in the shadow family from planes facing up (open sky)."""
    ndl = dot(N, L)
    eff = np.clip(ndl, 0, 1) * albedo
    x_light = 2.51 + 3.0 * eff ** light_gamma                     # 2.5..5.5 -> steps 3..5
    refl = np.clip(dot(N, bounce_dir(L)), 0, 1)
    up = np.clip(N[..., 1], 0, 1)
    x_shadow = 0.75 + bounce * refl ** 1.2 + sky * up
    # the core: faces neither light nor bounce
    x_shadow = x_shadow - core * np.exp(-((ndl + 0.10) / 0.16) ** 2) * (ndl < 0.1)
    x_shadow = np.clip(x_shadow, 0, 2.45)
    # occlusion: contact lines only, both families
    x_shadow = x_shadow * ao ** ao_power
    x_light = np.where(ao < contact, x_light * ao / contact * 0.3, x_light)
    # cast-shadow penumbra
    w = np.clip((shadow - 0.35) / 0.3, 0, 1) * (ndl > 0)
    x = x_shadow * (1 - w) + x_light * w
    lit = w > 0.5
    return np.clip(x / 5.0, 0, 1), lit, eff * shadow

def value_ramp():
    """six greys, spaced by eye rather than linearly: the shadow steps close together,
    a big jump at the terminator, the lights spread out."""
    g = np.array([0.07, 0.17, 0.27, 0.52, 0.70, 0.90]) * 255
    return np.stack([g, g, g], -1).round()
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && cat >> stage1_forms.py <<'EOF'

def dirs_sheet():
    """the same four forms under a LOW (12 deg) dusk sun from each task direction"""
    ramp = value_ramp(); ims = []
    for name, az in [("left", 10), ("front", 80), ("right", 170), ("back", 265)]:
        L = sun_dir(az, 12)
        for k in ["sphere", "cube", "cylinder", "slab"]:
            q, v, lit = render(k, L)
            ims.append(label(to_image(q, ramp, 2), f"{k} {name}"))
    grid(ims, 4).save(f"{OUT}/forms_lowsun_dirs.png")

if __name__ == "__main__":
    dirs_sheet()
EOF
nice -n 10 python3 stage1_forms.py 6

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && cat > geometry.py <<'EOF'
"""Geometry for the task: the boulder mass and the bank/water setting (as SDFs).

The boulder is designed the way the atelier teaches you to draw one:

  1. a big simple mass (a squat, slightly lopsided ellipsoid, the 'loaf'),
  2. broken into a few large PLANES (a convex set of cutting planes, rounded at the
     edges by weathering: granite exfoliates in shells, so edges round and faces stay broad),
  3. only then small irregularity (low-frequency lumps), and
  4. a couple of joint cracks cut into it.
"""
import numpy as np
from core import norm, dot, sd_ellipsoid, smin, smax

# ------------------------------------------------------------------ noise
def _hash3(ix, iy, iz, seed):
    h = (ix * 374761393 + iy * 668265263 + iz * 1274126177 + seed * 1442695041) & 0xFFFFFFFF
    h = ((h ^ (h >> 13)) * 1274126177) & 0xFFFFFFFF
    return ((h ^ (h >> 16)) & 0xFFFF) / 65535.0

def vnoise3(p, seed=0):
    """smooth value noise in [-1, 1]"""
    fl = np.floor(p); f = p - fl
    i = fl.astype(np.int64)
    u = f * f * (3 - 2 * f)
    out = 0
    for dx in (0, 1):
        for dy in (0, 1):
            for dz in (0, 1):
                w = (u[..., 0] if dx else 1 - u[..., 0]) * (u[..., 1] if dy else 1 - u[..., 1]) * \
                    (u[..., 2] if dz else 1 - u[..., 2])
                out = out + w * _hash3(i[..., 0] + dx, i[..., 1] + dy, i[..., 2] + dz, seed)
    return out * 2 - 1

def fbm3(p, seed=0, octaves=3, lac=2.03, gain=0.5):
    a, s, tot = 1.0, 0.0, 0.0
    for o in range(octaves):
        s = s + a * vnoise3(p, seed + o * 17)
        tot += a; a *= gain; p = p * lac
    return s / tot

def vnoise2(x, y, seed=0):
    return vnoise3(np.stack([x, y, np.zeros_like(x) + 0.37], -1), seed)

def fbm2(x, y, seed=0, octaves=3):
    return fbm3(np.stack([x, y, np.zeros_like(x) + 0.37], -1), seed, octaves)

# ------------------------------------------------------------------ boulder
class Boulder:
    def __init__(self, size=1.8, seed=0, center=(0, 0, 0), sink=0.28, n_planes=11,
                 round_k=None, lump=0.035, cracks=2):
        rng = np.random.default_rng(seed)
        self.seed = seed
        self.size = size
        s = size / 2
        # the loaf: wider than tall, a bit deeper than wide, lopsided
        self.r = np.array([s * rng.uniform(1.0, 1.15), s * rng.uniform(0.66, 0.8), s * rng.uniform(0.85, 1.0)])
        self.c = np.array(center, float) + np.array([0, self.r[1] * (1 - 2 * sink), 0])
        yaw = rng.uniform(0, np.pi)
        self.rot = np.array([[np.cos(yaw), 0, np.sin(yaw)], [0, 1, 0], [-np.sin(yaw), 0, np.cos(yaw)]])
        # cutting planes: normals spread over the upper hemisphere and sides
        ns, offs = [], []
        for i in range(n_planes):
            th = rng.uniform(0, 2 * np.pi)
            el = rng.uniform(-0.25, 1.0) if i else 1.0                 # first plane = the top
            n = norm([np.cos(th) * np.sqrt(1 - el * el) * 1.0, el, np.sin(th) * np.sqrt(1 - el * el)])
            # how deep it cuts: the ellipsoid support distance, shaved by 8..22 %
            sup = np.sqrt(np.sum((self.r * n) ** 2))
            ns.append(n); offs.append(sup * rng.uniform(0.78, 0.92) if i else sup * rng.uniform(0.80, 0.88))
        self.ns = np.array(ns); self.offs = np.array(offs)
        self.k = round_k if round_k is not None else 0.16 * s
        self.lump = lump * size
        # cracks: thin slabs through the mass, nearly vertical
        self.cracks = []
        for i in range(cracks):
            th = rng.uniform(0, np.pi)
            n = norm([np.cos(th), rng.uniform(-0.3, 0.3), np.sin(th)])
            self.cracks.append((n, rng.uniform(-0.35, 0.35) * s, rng.uniform(0.012, 0.02) * size))

    def local(self, p):
        return (p - self.c) @ self.rot   # into boulder frame

    def sdf(self, p):
        q = self.local(p)
        d = sd_ellipsoid(q, np.zeros(3), self.r)
        for n, o in zip(self.ns, self.offs):
            d = smax(d, q @ n - o, self.k)
        if self.lump > 0:
            d = d + self.lump * fbm3(q / (self.size * 0.45), self.seed + 99, 3)
        for n, o, w in self.cracks:
            # a V-groove, only on the upper part of the mass
            g = np.abs(q @ n - o) - w
            depth = 0.06 * self.size
            groove = -(depth - np.maximum(g, 0) * 3.0)
            d = np.where(q[..., 1] > -self.r[1] * 0.2, np.maximum(d, np.minimum(-groove, d + depth) * 0 + np.maximum(d, -(depth - np.maximum(g, 0) * 2.5))), d)
        return d
EOF
echo ok

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
i=s.index("        for n, o, w in self.cracks:")
s=s[:i]+"""        for n, o, w in self.cracks:
            # a shallow joint: the thin band around the crack plane, near the surface only
            band = np.abs(q @ n - o) - w
            carve = np.maximum(band, -(d + 0.05 * self.size))
            d = np.maximum(d, -carve)
        return d
"""
open('geometry.py','w').write(s)
EOF
cat > stage2_boulder.py <<'EOF'
"""Stage 2 -- the boulder as form, greyscale. Low dusk sun from each task direction.
Outputs stage2/*.png
"""
import os, sys
import numpy as np
from core import *
from valuemodel import form_value, value_ramp
from geometry import Boulder

OUT = "stage2"; os.makedirs(OUT, exist_ok=True)
W, H = 160, 110

def render_boulder(b, L, dist=4.5, fov=38, levels=6, notan=False, planes_only=False):
    def scene(p):
        g = p[..., 1]
        o = b.sdf(p)
        return np.minimum(g, o), np.where(o < g, 1, 0)
    tgt = np.array([0, b.size * 0.3, 0])
    cam = Camera(W, H, pos=tgt + np.array([0, dist * np.tan(np.radians(10)), -dist]), target=tgt, fov_deg=fov)
    o, d = cam.rays()
    t, hit, p, mat = march(scene, o, d, tmax=80, steps=300, eps=3e-4)
    N = normals(scene, p, e=1e-3 * b.size)
    sh = soft_shadow(scene, p + N * 2e-3, L, k=20)
    ao = ambient_occ(scene, p, N, scale=0.02 * b.size)
    albedo = np.where(mat == 0, 0.6, 0.95)
    v, lit, _ = form_value(N, L, sh, ao, albedo, sky=0.35)
    v = np.where(hit, v, 0.42)
    if notan:
        return np.where(v >= 0.5, 1, 0)
    return quantize(v, levels, kind="check", strength=0.5, flat_snap=0.08)

if __name__ == "__main__":
    ramp = value_ramp()
    dirs = [("left", 8), ("front", 80), ("right", 172), ("back", 262)]
    for seed in [1, 2, 3]:
        b = Boulder(size=1.8, seed=seed)
        ims = []
        for name, az in dirs:
            L = sun_dir(az, 12)
            ims.append(label(to_image(render_boulder(b, L, notan=True), np.array([[40]*3, [215]*3]), 2), f"s{seed} {name} notan"))
        for name, az in dirs:
            L = sun_dir(az, 12)
            ims.append(label(to_image(render_boulder(b, L), ramp, 2), f"s{seed} {name}"))
        grid(ims, 4).save(f"{OUT}/boulder_seed{seed}.png")
    print("ok")
EOF
time nice -n 10 python3 stage2_boulder.py

# output:
Command did not complete within its 120s timeout and was moved to the background (ID: bb6rc65hn). Output is being written to: (tmp) You will be notified when it completes. To check interim output, use Read on that file path.
Session cwd remains ~/work/pixelart-studies/moodboards; directory changes made by the backgrounded command do not apply to subsequent commands.

sleep 60; cat (tmp) ls ~/work/pixelart-studies/pedagogy/atelier/stage2

# output:
<tool_use_error>Blocked: sleep 60 followed by: cat (tmp) ls ~/work/pixelart-studies/pedagogy/atelier/stage2. To wait for a condition, use Monitor with an until-loop (e.g. `until <check>; do sleep 2; done`). To wait for a command you started, use run_in_background: true. Do not chain shorter sleeps to work around this block.</tool_use_error>

cd ~/work/pixelart-studies/pedagogy/atelier && cat >> geometry.py <<'EOF'

# ------------------------------------------------------------------ the setting
class Setting:
    """A cove. Water plane y=0. The boulder sits at the lip of the LEFT bank; the camera
    stands on the near bank, off to the right, so the line of sight to the boulder's foot
    crosses water and its reflection hangs below it. A far bank closes the water.

    Everything is laid out in a frame attached to the boulder (origin at the boulder foot
    on the shoreline), then the camera is placed by (distance, yaw)."""
    def __init__(self, boulder_size=1.8, distance=10.0, seed=0, yaw_deg=32.0,
                 lip=0.22, slope=0.07, hummock=0.18, hummock_scale=1.6):
        rng = np.random.default_rng(seed + 1000)
        self.seed = seed
        self.D = distance
        self.lip, self.slope = lip, slope
        self.hum, self.hscale = hummock, hummock_scale
        self.bsize = boulder_size
        # the camera sits at yaw to the right of the boulder's -z axis
        yaw = np.radians(yaw_deg + rng.uniform(-6, 6))
        self.cam_h = np.array([np.sin(yaw), 0, -np.cos(yaw)]) * distance       # camera ground position
        self.fwd = -self.cam_h / distance                                        # horizontal view dir
        self.shore_wobble = rng.uniform(0, 100, 4)
        # far bank beyond the boulder, measured along the view direction
        self.far = distance * 0.9 + 14 + rng.uniform(0, 8)
        self.near = 0.45 * distance                                              # near bank reach from camera
        self.boulder = Boulder(size=boulder_size, seed=seed,
                               center=(-0.25 * boulder_size, lip, 0.0))

    # signed distance to the shore, positive on land (in metres, roughly)
    def shore(self, x, z):
        w = self.shore_wobble
        # left bank: land for x < curve(z); the boulder's promontory bulges into the water
        curve = 0.15 * self.bsize + 0.9 * np.sin(z * 0.23 + w[0]) + 0.5 * np.sin(z * 0.61 + w[1]) \
                - 0.035 * z * z / (1 + 0.02 * z * z) * 6
        prom = 0.6 * self.bsize * np.exp(-(z / (1.4 * self.bsize)) ** 2)
        left = (curve + prom) - x
        # near bank around the camera
        rel = np.stack([x - self.cam_h[0], z - self.cam_h[2]], -1)
        along = rel[..., 0] * self.fwd[0] + rel[..., 1] * self.fwd[2]
        across = -rel[..., 0] * self.fwd[2] + rel[..., 1] * self.fwd[0]
        near = self.near + 0.08 * self.D * np.sin(across * 0.4 / max(1, self.D ** 0.5) + w[2]) - along \
               - 0.004 * across ** 2
        # far bank
        along0 = x * self.fwd[0] + z * self.fwd[2]
        far = along0 - self.far + 2.0 * np.sin(x * 0.13 + w[3])
        return np.maximum(np.maximum(left, near), far)

    def height(self, x, z):
        s = self.shore(x, z)
        hs = self.hscale
        hum = self.hum * (0.65 * fbm2(x / hs, z / hs, self.seed + 5, 3)
                          + 0.35 * fbm2(x / (hs * 4), z / (hs * 4), self.seed + 7, 2))
        inland = np.clip(s, 0, None)
        h_land = self.lip + self.slope * inland + hum * np.clip(s / 0.6, 0, 1) + 0.012 * inland ** 1.3 * 0.0
        h_water = self.lip + s * 3.0 - 0.02                                      # the cut turf bank
        return np.where(s > 0, h_land, np.maximum(h_water, -1.0))

    def terrain_sdf(self, p):
        return (p[..., 1] - self.height(p[..., 0], p[..., 2])) * 0.3

    def scene(self, p):
        """solid scene without the water: (dist, mat) mat 0 = land, 1 = boulder"""
        g = self.terrain_sdf(p)
        o = self.boulder.sdf(p)
        return np.minimum(g, o), np.where(o < g, 1, 0)
EOF
echo ok

# output:
ok

cat (tmp) ls ~/work/pixelart-studies/pedagogy/atelier/stage2

# output:
boulder_seed1.png

pkill -f stage2_boulder.py; cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("""            ns.append(n); offs.append(sup * rng.uniform(0.78, 0.92) if i else sup * rng.uniform(0.80, 0.88))""","""            ns.append(n); offs.append(sup * rng.uniform(0.70, 0.86) if i else sup * rng.uniform(0.74, 0.84))""")
s=s.replace("self.k = round_k if round_k is not None else 0.16 * s","self.k = round_k if round_k is not None else 0.09 * s")
s=s.replace("lump=0.035, cracks=2","lump=0.03, cracks=1")
s=s.replace("""            self.cracks.append((n, rng.uniform(-0.35, 0.35) * s, rng.uniform(0.012, 0.02) * size))""","""            # a joint runs only part-way: centre and half-length along the crack plane
            mid = rng.uniform(-0.4, 0.4, 3) * self.r
            self.cracks.append((n, float(mid @ n), rng.uniform(0.006, 0.011) * size, mid, 0.55 * s))""")
s=s.replace("""        for n, o, w in self.cracks:
            # a shallow joint: the thin band around the crack plane, near the surface only
            band = np.abs(q @ n - o) - w
            carve = np.maximum(band, -(d + 0.05 * self.size))""","""        for n, o, w, mid, half in self.cracks:
            # a shallow joint: thin band around the crack plane, near the surface only,
            # tapering out at both ends (a joint rarely runs clean across a boulder)
            along = np.linalg.norm((q - mid) - np.outer(q.reshape(-1, 3) @ n - o, n).reshape(q.shape), axis=-1)
            taper = np.clip(1 - along / half, 0, 1)
            band = np.abs(q @ n - o) - w * taper
            carve = np.maximum(np.maximum(band, -(d + 0.035 * self.size * taper)), -(taper - 0.05))""")
open('geometry.py','w').write(s)
s=open('stage2_boulder.py').read()
s=s.replace("for seed in [1, 2, 3]:","for seed in [int(a) for a in sys.argv[1:]] or [1, 2, 3]:")
open('stage2_boulder.py','w').write(s)
EOF
time nice -n 10 python3 stage2_boulder.py 1

# output:
ok
nice -n 10 python3 stage2_boulder.py 1  1016.99s user 554.07s system 960% cpu 2:43.51 total

cd ~/work/pixelart-studies/pedagogy/atelier && env | grep -i thread; python3 -c "import numpy; numpy.show_config()" 2>&1 | grep -iA3 "blas\|lapack" | head -20; python3 -c "import threadpoolctl; print('tpc ok')"

# output:
blas:
    detection method: pkgconfig
    found: true
    include directory: /opt/arm64-builds/include
--
    name: openblas64
    openblas configuration: USE_64BITINT=1 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SANDYBRIDGE MAX_THREADS=3
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.23.dev
  lapack:
    detection method: internal
    found: true
    include directory: unknown
--
    openblas configuration: unknown
    pc file directory: unknown
    version: 1.26.4
Compilers:
--
tpc ok

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s=s.replace('''os.environ.setdefault("OMP_NUM_THREADS", "4")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "4")
os.environ.setdefault("VECLIB_MAXIMUM_THREADS", "4")
os.environ.setdefault("MKL_NUM_THREADS", "4")
import numpy as np''','''for _k in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "VECLIB_MAXIMUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[_k] = "4"          # the studio rule: at most 4 threads
import numpy as np
try:
    from threadpoolctl import threadpool_limits
    threadpool_limits(4)
except Exception:
    pass''')
open('core.py','w').write(s)
s=open('stage2_boulder.py').read()
s=s.replace("import os, sys\nimport numpy as np\nfrom core import *","from core import *\nimport os, sys")
open('stage2_boulder.py','w').write(s)
EOF
head -8 stage2_boulder.py

# output:
"""Stage 2 -- the boulder as form, greyscale. Low dusk sun from each task direction.
Outputs stage2/*.png
"""
from core import *
import os, sys
from valuemodel import form_value, value_ramp
from geometry import Boulder

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s=s.replace('''    os.environ[_k] = "4"          # the studio rule: at most 4 threads''','''    os.environ[_k] = os.environ.get("ATELIER_THREADS", "1")   # studio rule: <= 4 threads in total;
                                                          # I run up to 4 single-thread processes''')
s=s.replace("threadpool_limits(4)","threadpool_limits(int(os.environ.get(\"ATELIER_THREADS\", \"1\")))")
open('core.py','w').write(s)
s=open('geometry.py').read()
s=s.replace("lump=0.03, cracks=1","lump=0.014, cracks=1")
open('geometry.py','w').write(s)
s=open('stage2_boulder.py').read()
s=s.replace("W, H = 160, 110","W, H = 200, 130")
s=s.replace("def render_boulder(b, L, dist=4.5, fov=38,","def render_boulder(b, L, dist=4.5, fov=30,")
open('stage2_boulder.py','w').write(s)
EOF
for s in 1 2 3; do (nice -n 10 python3 stage2_boulder.py $s > /dev/null &) ; done; echo started

# output:
started

cd ~/work/pixelart-studies/pedagogy/atelier && cat > scene.py <<'EOF'
"""The task scene, rendered as buffers (value, material, depth, normals, light) that the
painting stages consume. Stage 3 stops at values; stage 4/5 paint on top of the buffers.
"""
from core import *
from valuemodel import form_value, bounce_dir
from geometry import Setting, fbm2, vnoise2

W, H = 320, 180
LIGHTS = {"left": 8, "front": 82, "right": 172, "back": 262}
SIZES = {"small": 1.0, "medium": 1.8, "large": 2.8}
DISTS = {"close": None, "10m": 10.0, "40m": 40.0}

# material ids
SKY, MTN, TREES, LAND, BOULDER, WATER = 0, 1, 2, 3, 4, 5

class Params:
    def __init__(self, size="medium", light="left", dist="10m", seed=0, sun_elev=9.0):
        self.size, self.light, self.dist, self.seed = size, light, dist, seed
        self.bsize = SIZES[size] if isinstance(size, str) else float(size)
        D = DISTS[dist] if isinstance(dist, str) else float(dist)
        self.D = D if D is not None else max(3.2, 2.2 * self.bsize)     # "close": a few boulder-widths
        self.az = LIGHTS[light] if isinstance(light, str) else float(light)
        self.sun_elev = sun_elev

def make_camera(st, P):
    b = st.boulder
    eye_h = max(1.5, st.lip + P.D * np.tan(np.radians(10)) + b.size * 0.35)
    pos = st.cam_h + np.array([0, eye_h, 0])
    # aim at the boulder, then swing a little so it sits off-centre (left third)
    tgt = b.c + np.array([0, 0.0, 0])
    f = norm(tgt - pos); fh = norm([f[0], 0, f[2]])
    right = np.cross([0, 1, 0], fh)
    swing = 0.16 * np.linalg.norm(tgt - pos)
    tgt2 = tgt + right * swing
    # pitch exactly 10 deg down
    d = norm(tgt2 - pos); dh = norm([d[0], 0, d[2]])
    look = norm(dh * np.cos(np.radians(10)) - np.array([0, np.sin(np.radians(10)), 0]))
    fov = 52.0
    return Camera(W, H, pos, pos + look, fov_deg=fov)

def world_sun(cam, az, elev):
    """task light directions are relative to the viewer: rotate into the world by camera yaw"""
    Lc = sun_dir(az, elev)                       # camera frame: x right, z forward
    fh = norm([cam.f[0], 0, cam.f[2]]); rh = np.cross([0, 1, 0], fh)
    return norm(Lc[0] * rh + Lc[1] * np.array([0, 1, 0]) + Lc[2] * fh)

# ------------------------------------------------------------------ backdrop (env map)
class Backdrop:
    """Sky, a far mountain wall and a far treeline, as a function of ray DIRECTION only.
    Because it depends on direction only, the water's reflection of it is automatically the
    mirror image about the horizon, as it is for distant things."""
    def __init__(self, st, cam, L, seed):
        self.L = L; self.seed = seed
        self.fh = norm([cam.f[0], 0, cam.f[2]])
        # elevation of the far-bank top seen from the eye (trees stand on it)
        eye = cam.pos
        self.e_bank = np.arctan2(st.lip + 0.3 - eye[1], st.far + np.linalg.norm(st.cam_h) * 0.0 + st.D * 0.1)
        self.tree_h = np.arctan2(9.0, st.far + st.D)        # ~9 m spruces on the far bank
        self.cam = cam

    def angles(self, d):
        az = np.arctan2(d[..., 0] * self.fh[2] - d[..., 2] * self.fh[0], d[..., 0] * self.fh[0] + d[..., 2] * self.fh[2])
        el = np.arcsin(np.clip(d[..., 1], -1, 1))
        return az, el

    def mountains(self, az):
        s = self.seed
        a = az * 3.0
        base = np.radians(4.0) + np.radians(5.5) * np.clip(0.55 + 0.6 * fbm2(a, np.full_like(a, 3.1), s + 21, 4), 0, 1.3)
        # ridge-ish: sharpen peaks
        ridge = 1 - np.abs(fbm2(a * 2.3, np.full_like(a, 7.7), s + 22, 3))
        return base + np.radians(2.2) * ridge ** 3

    def trees(self, az):
        """spruce silhouettes: a comb of narrow triangles of varied height"""
        k = 520.0                                    # trees per radian
        u = az * k
        i = np.floor(u); f = u - i
        rnd = lambda j, o: (np.sin(j * 12.9898 + o * 78.233 + self.seed * 3.1) * 43758.5453) % 1.0
        hgt = (0.45 + 0.75 * rnd(i, 1)) * (0.6 + 0.4 * (vnoise2(az * 9, np.zeros_like(az), self.seed + 5) * 0.5 + 0.5))
        gap = rnd(i, 2) < 0.18
        cx = 0.3 + 0.4 * rnd(i, 3)
        tri = np.clip(1 - np.abs(f - cx) / 0.5, 0, 1)         # triangular crown
        e = self.e_bank + self.tree_h * hgt * tri ** 0.8 * np.where(gap, 0.25, 1.0)
        return e

    def sample(self, d):
        """returns (value in [0,1], material, extra) for directions d"""
        az, el = self.angles(d)
        Lh = norm([self.L[0], 0, self.L[2]])
        dh = norm(np.stack([d[..., 0], np.zeros_like(el), d[..., 2]], -1))
        toward_sun = np.clip(dot(dh, Lh), -1, 1)
        # sky: light at the horizon, darker to the zenith, brighter toward the sun
        e = np.clip(el / np.radians(30), 0, 1)
        v = 0.86 - 0.28 * e ** 0.8 + 0.08 * np.clip(toward_sun, 0, 1) ** 3 * (1 - e)
        mat = np.full(el.shape, SKY)
        em = self.mountains(az)
        m = el < em
        # mountain facets: a fake normal from the ridge profile + a 2D noise for gullies
        dA = 1e-3
        slope = (self.mountains(az + dA) - self.mountains(az - dA)) / (2 * dA)
        gully = fbm2(az * 40, el * 40, self.seed + 30, 3)
        nx = -np.clip(slope * 6 + gully * 0.8, -1.5, 1.5)
        # sun from the left lights left-facing slopes (nx<0 in view space)
        lx = -dot(self.L, np.cross([0, 1, 0], self.fh))       # >0 when sun is on the left
        facing = np.clip(0.5 + 0.5 * (-nx) * np.sign(lx) * min(1, abs(lx) * 2), 0, 1)
        snow = (fbm2(az * 25, el * 25, self.seed + 31, 3) > 0.05) & (el > em - np.radians(3.0))
        vm = np.where(facing > 0.55, 0.62, 0.40) + 0.10 * snow
        vm = vm * 0.8 + v * 0.2                                # aerial perspective: toward sky
        v = np.where(m, vm, v); mat = np.where(m, MTN, mat)
        et = self.trees(az)
        t = el < et
        v = np.where(t, 0.14 + 0.05 * (np.clip(toward_sun, 0, 1)), v)
        mat = np.where(t, TREES, mat)
        return v, mat

# ------------------------------------------------------------------ render
def render_buffers(P, Wd=W, Hd=H):
    st = Setting(boulder_size=P.bsize, distance=P.D, seed=P.seed)
    cam = make_camera(st, P)
    cam.W, cam.H = Wd, Hd
    L = world_sun(cam, P.az, P.sun_elev)
    bd = Backdrop(st, cam, L, P.seed)
    sc = st.scene
    o, d = cam.rays()
    tmax = st.far + st.D + 80
    t, hit, p, mat = march(sc, o, d, tmax=tmax, steps=360, eps=4e-4)
    # water plane
    tw = np.where(d[..., 1] < 0, -o[..., 1] / np.minimum(d[..., 1], -1e-6), np.inf)
    water = (tw < np.where(hit, t, np.inf)) & (tw < tmax * 1.5)
    buf = dict(cam=cam, L=L, st=st, P=P)
    # primary shading
    v, m, extra = shade(sc, st, p, d, t, hit, mat, L, bd)
    # reflections
    pw = o + tw[..., None] * d
    dr = d * np.array([1, -1, 1])
    idx = np.nonzero(water)
    if len(idx[0]):
        po = pw[idx] + dr[idx] * 1e-3
        tr, hr, prr, mr = march(sc, po, dr[idx], tmax=tmax, steps=300, eps=6e-4)
        vr, mmr, exr = shade(sc, st, prr, dr[idx], tr, hr, mr, L, bd)
        rv = np.zeros(t.shape); rv[idx] = vr
        rmat = np.full(t.shape, -1); rmat[idx] = mmr
        rdist = np.zeros(t.shape); rdist[idx] = np.where(hr, tr, 1e4)
    else:
        rv = np.zeros(t.shape); rmat = np.full(t.shape, -1); rdist = np.zeros(t.shape)
    buf.update(v=v, mat=m, depth=np.where(hit, t, 1e4), p=p, water=water, tw=tw,
               rv=rv, rmat=rmat, rdist=rdist, **extra)
    buf["mat"] = np.where(water, WATER, buf["mat"])
    buf["depth"] = np.where(water, tw, buf["depth"])
    return buf

def shade(sc, st, p, d, t, hit, mat, L, bd):
    shp = t.shape
    v = np.zeros(shp); m = np.full(shp, SKY)
    N = np.zeros(p.shape); N[..., 1] = 1
    sh = np.ones(shp); ao = np.ones(shp); lit = np.zeros(shp, bool); lamt = np.zeros(shp)
    h = np.nonzero(hit)
    if len(h[0]):
        ph = p[h]
        Nh = normals(sc, ph, e=2e-3 * max(1, st.bsize))
        shh = soft_shadow(sc, ph + Nh * 3e-3, L, k=22, tmax=40)
        aoh = ambient_occ(sc, ph, Nh, scale=0.025 * st.bsize)
        is_b = mat[h] == 1
        albedo = np.where(is_b, 0.95, 0.72)
        vh, lith, lah = form_value(Nh, L, shh, aoh, albedo, sky=np.where(is_b, 0.35, 0.25))
        v[h] = vh; N[h] = Nh; sh[h] = shh; ao[h] = aoh; lit[h] = lith; lamt[h] = lah
        m[h] = np.where(is_b, BOULDER, LAND)
    miss = ~hit
    if miss.any():
        vb, mb = bd.sample(d[miss])
        v[miss] = vb; m[miss] = mb
    return v, m, dict(N=N, shadow=sh, ao=ao, lit=lit, lamt=lamt)
EOF
echo ok

# output:
ok

ls -la ~/work/pixelart-studies/pedagogy/atelier/stage2; pgrep -f stage2_boulder | wc -l

# output:
total 96
drwxr-xr-x   5 user  staff    160 Sep 25 18:17 .
drwxr-xr-x  12 user  staff    384 Sep 25 18:18 ..
-rw-r--r--   1 user  staff  14645 Sep 25 18:17 boulder_seed1.png
-rw-r--r--   1 user  staff  14402 Sep 25 18:17 boulder_seed2.png
-rw-r--r--   1 user  staff  14111 Sep 25 18:17 boulder_seed3.png
       0

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("sink=0.28, n_planes=11,","sink=0.40, n_planes=9,")
s=s.replace("""            ns.append(n); offs.append(sup * rng.uniform(0.70, 0.86) if i else sup * rng.uniform(0.74, 0.84))""","""            ns.append(n); offs.append(sup * rng.uniform(0.64, 0.84) if i else sup * rng.uniform(0.70, 0.82))""")
s=s.replace("self.k = round_k if round_k is not None else 0.09 * s","self.k = round_k if round_k is not None else 0.065 * s")
# lopsided mass: a second, smaller lobe
s=s.replace("""        yaw = rng.uniform(0, np.pi)""","""        # a second lobe makes the mass lopsided (a shoulder), as real boulders are
        ang = rng.uniform(0, 2 * np.pi)
        self.lobe_c = np.array([np.cos(ang) * self.r[0] * 0.45, -self.r[1] * rng.uniform(0.0, 0.25), np.sin(ang) * self.r[2] * 0.4])
        self.lobe_r = self.r * np.array([rng.uniform(0.55, 0.75), rng.uniform(0.6, 0.85), rng.uniform(0.55, 0.75)])
        yaw = rng.uniform(0, np.pi)""")
s=s.replace("""        d = sd_ellipsoid(q, np.zeros(3), self.r)
        for n, o in zip(self.ns, self.offs):""","""        d = smin(sd_ellipsoid(q, np.zeros(3), self.r), sd_ellipsoid(q, self.lobe_c, self.lobe_r), 0.2 * self.size)
        for n, o in zip(self.ns, self.offs):""")
open('geometry.py','w').write(s)
s=open('stage2_boulder.py').read()
s=s.replace("def render_boulder(b, L, dist=4.5, fov=30,","def render_boulder(b, L, dist=4.5, fov=38,")
open('stage2_boulder.py','w').write(s)
EOF
for s in 1 2 3 4; do (nice -n 10 python3 stage2_boulder.py $s > /dev/null 2>stage2/err$s.txt &) ; done; echo started

# output:
started

cd ~/work/pixelart-studies/pedagogy/atelier && cat > stage3_setting.py <<'EOF'
"""Stage 3 -- the setting as form, greyscale.
Masses: the bank's hummocks, the turf lip at the water, the boulder sunk in it, the water
plane carrying the reflection, the far bank + treeline + mountains as flat value shapes.
"""
from core import *
import os, sys, pickle
from scene import *

OUT = "stage3"; os.makedirs(OUT, exist_ok=True)
G = np.array([0.05, 0.12, 0.20, 0.30, 0.42, 0.54, 0.66, 0.78, 0.88, 0.96])
GREY10 = np.stack([G * 255] * 3, -1).round()

# each material's six-step value plan placed inside the picture's global value scale
STEP_TO_G = {
    BOULDER: [0.3, 1.2, 2.4, 5.0, 7.0, 8.6],
    LAND:    [0.2, 1.0, 2.0, 4.0, 5.6, 6.8],
}

def to_global(v, mat):
    """-> continuous index into G"""
    gi = np.interp(v, G, np.arange(10))                 # backdrop values are already 'global'
    for m, tab in STEP_TO_G.items():
        x = v * 5
        gm = np.interp(x, np.arange(6), tab)
        gi = np.where(mat == m, gm, gi)
    return gi

def value_image(buf, water_dark=0.9):
    gi = to_global(buf["v"], buf["mat"])
    # water: the reflection, a touch darker and flatter than the thing reflected
    rgi = to_global(buf["rv"], buf["rmat"])
    rgi = rgi * water_dark + 0.25
    # aerial perspective on land/boulder far away: pull toward the sky's value
    far = 1 - np.exp(-np.clip(buf["depth"], 0, 1e3) / 90.0)
    solid = np.isin(buf["mat"], [LAND, BOULDER])
    gi = np.where(solid, gi * (1 - 0.5 * far) + 6.5 * 0.5 * far, gi)
    gi = np.where(buf["water"], rgi, gi)
    return gi

def quantize_global(gi):
    v = np.clip(gi / 9.0, 0, 1)
    return quantize(v, 10, kind="check", strength=0.55, flat_snap=0.12)

if __name__ == "__main__":
    args = sys.argv[1:]
    size = args[0] if args else "medium"; light = args[1] if len(args) > 1 else "left"
    dist = args[2] if len(args) > 2 else "10m"; seed = int(args[3]) if len(args) > 3 else 0
    P = Params(size, light, dist, seed)
    buf = render_buffers(P)
    gi = value_image(buf)
    q = quantize_global(gi)
    name = f"{size}_{light}_{dist}_s{seed}"
    to_image(q, GREY10, 3).save(f"{OUT}/setting_{name}.png")
    # the massing diagram: materials as flat values (what's in front of what)
    mm = np.array([[210]*3, [150]*3, [40]*3, [110]*3, [230, 120, 90], [170, 180, 210]])
    to_image(buf["mat"], mm, 2).save(f"{OUT}/masses_{name}.png")
    keep = {k: buf[k] for k in ["v", "mat", "depth", "water", "rv", "rmat", "rdist", "N", "shadow", "ao", "lit", "lamt", "tw", "p"]}
    keep["L"] = buf["L"]
    with open(f"{OUT}/buf_{name}.pkl", "wb") as f:
        pickle.dump(keep, f)
    print(name)
EOF
echo ok

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && until [ $(pgrep -f stage2_boulder | wc -l) -eq 0 ]; do sleep 5; done; cat stage2/err*.txt; ls stage2

# output:
boulder_seed1.png
boulder_seed2.png
boulder_seed3.png
boulder_seed4.png
err1.txt
err2.txt
err3.txt
err4.txt

cd ~/work/pixelart-studies/pedagogy/atelier && rm stage2/err*.txt && python3 -c "
from PIL import Image
ims=[Image.open(f'stage2/boulder_seed{s}.png') for s in [1,2,3,4]]
# keep only the lit rows (bottom half) for a compact contact sheet
rows=[im.crop((0,im.height//2,im.width,im.height)) for im in ims]
out=Image.new('RGB',(rows[0].width,sum(r.height for r in rows)))
y=0
for r in rows: out.paste(r,(0,y)); y+=r.height
out.save('stage2/boulder_v2_all.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('core.py').read()
s+='''
def remove_orphans(a, passes=1):
    """a pixel whose 8 neighbours all agree on a different value is a stray: adopt theirs.
    (checker dither survives: a checker pixel's diagonal neighbours equal itself.)"""
    a = a.copy()
    for _ in range(passes):
        pad = np.pad(a, 1, mode="edge")
        H, W = a.shape
        nb = [pad[dy:dy + H, dx:dx + W] for dy in range(3) for dx in range(3) if not (dy == 1 and dx == 1)]
        same = np.all([n == nb[0] for n in nb[1:]], axis=0)
        fix = same & (nb[0] != a)
        a[fix] = nb[0][fix]
    return a
'''
open('core.py','w').write(s)
s=open('geometry.py').read()
s=s.replace("""            ns.append(n); offs.append(sup * rng.uniform(0.64, 0.84) if i else sup * rng.uniform(0.70, 0.82))""","""            ns.append(n); offs.append(sup * rng.uniform(0.68, 0.86) if i else sup * rng.uniform(0.72, 0.84))""")
s=s.replace("self.k = round_k if round_k is not None else 0.065 * s","self.k = round_k if round_k is not None else rng.uniform(0.06, 0.12) * s")
open('geometry.py','w').write(s)
s=open('stage3_setting.py').read()
s=s.replace("""    return quantize(v, 10, kind="check", strength=0.55, flat_snap=0.12)""","""    return remove_orphans(quantize(v, 10, kind="check", strength=0.55, flat_snap=0.12))""")
open('stage3_setting.py','w').write(s)
EOF
time nice -n 10 python3 stage3_setting.py medium left 10m 0

# output:
medium_left_10m_s0
nice -n 10 python3 stage3_setting.py medium left 10m 0  13.98s user 3.64s system 98% cpu 17.977 total

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
i=s.index("class Setting:")
s=s[:i]+'''class Setting:
    """A tarn shore. Water plane y=0.

    Laid out in the VIEW frame of the boulder: 'a' runs across the view (right +),
    'b' runs along it (away from the camera +), origin at the boulder foot.
      * the left bank: land where b > shore(a). The shore passes just in front of the
        boulder's foot and runs diagonally -- toward the camera on the left, receding on
        the right -- so the boulder sits on the lip with open water in front of it
        (its reflection hangs below it) and to its right.
      * the near bank: the strip the camera stands on (bottom of frame).
      * the far bank is far (a tarn ~100-150 m across) and lives in the backdrop.
    """
    def __init__(self, boulder_size=1.8, distance=10.0, seed=0, yaw_deg=0.0,
                 lip=0.22, slope=0.06, hummock=0.20, hummock_scale=1.7):
        rng = np.random.default_rng(seed + 1000)
        self.seed = seed
        self.D = distance
        self.lip, self.slope = lip, slope
        self.hum, self.hscale = hummock, hummock_scale
        self.bsize = boulder_size
        yaw = np.radians(yaw_deg + rng.uniform(-8, 8))
        self.fwd = np.array([np.sin(yaw), 0, np.cos(yaw)])          # view direction (horizontal)
        self.rgt = np.array([self.fwd[2], 0, -self.fwd[0]])
        self.cam_h = -self.fwd * distance                              # camera ground position
        self.w = rng.uniform(0, 100, 6)
        self.diag = rng.uniform(0.28, 0.5)                             # how steeply the shore recedes
        self.far = 100 + rng.uniform(0, 50)
        self.near = 0.42 * distance
        self.boulder = Boulder(size=boulder_size, seed=seed, center=(0.0, lip, 0.0))
        # where the boulder's front sits along b: put the lip just in front of it
        self.front = self.boulder.r[2] * 0.55

    def ab(self, x, z):
        return x * self.rgt[0] + z * self.rgt[2], x * self.fwd[0] + z * self.fwd[2]

    def shore(self, x, z):
        """signed distance-ish to the shoreline, positive on land (metres)"""
        a, b = self.ab(x, z)
        w = self.w; s = self.bsize
        wob = 0.35 * s * np.sin(a / (1.3 * s) + w[0]) + 0.18 * s * np.sin(a / (0.55 * s) + w[1]) \\
              + 0.6 * np.sin(a * 0.11 + w[2])
        wob = wob * np.clip(np.abs(a) / (1.2 * s), 0, 1)             # keep the lip clean at the boulder
        b_shore = -self.front + self.diag * a + wob
        left = b - b_shore
        # the camera's own bank
        across = a; along = b + self.D
        near = self.near - along + 0.06 * self.D * np.sin(across / max(2, 0.3 * self.D) + w[3]) \\
               - 0.02 * across ** 2 / max(1, self.D * 0.1)
        return np.maximum(left * np.cos(np.arctan(self.diag)), near)

    def height(self, x, z):
        s = self.shore(x, z)
        hs = self.hscale
        hum = self.hum * (0.65 * fbm2(x / hs, z / hs, self.seed + 5, 3)
                          + 0.35 * fbm2(x / (hs * 4), z / (hs * 4), self.seed + 7, 2))
        inland = np.clip(s, 0, None)
        h_land = self.lip + self.slope * inland + hum * np.clip(s / 0.8, 0, 1)
        h_water = self.lip + s * 2.5 - 0.02                            # the cut turf bank
        return np.where(s > 0, h_land, np.maximum(h_water, -1.0))

    def terrain_sdf(self, p):
        return (p[..., 1] - self.height(p[..., 0], p[..., 2])) * 0.35

    def scene(self, p):
        """solid scene without the water: (dist, mat) mat 0 = land, 1 = boulder"""
        g = self.terrain_sdf(p)
        o = self.boulder.sdf(p)
        return np.minimum(g, o), np.where(o < g, 1, 0)
'''
open('geometry.py','w').write(s)
EOF
echo ok

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('scene.py').read()
a=s.index("def make_camera"); b=s.index("# ------------------------------------------------------------------ render")
s=s[:a]+'''def make_camera(st, P):
    b = st.boulder
    eye_h = max(1.4, st.lip + P.D * np.tan(np.radians(10)) + b.size * 0.3)
    pos = st.cam_h + np.array([0, eye_h, 0])
    # compose: the boulder a little left of centre and below the middle, water to its right
    tgt = b.c + st.rgt * (0.14 * P.D + 0.2 * b.size)
    d = norm(tgt - pos); dh = norm([d[0], 0, d[2]])
    look = norm(dh * np.cos(np.radians(10)) - np.array([0, np.sin(np.radians(10)), 0]))   # 10 deg down
    return Camera(W, H, pos, pos + look, fov_deg=52.0)

def world_sun(cam, az, elev):
    """task light directions are relative to the viewer: rotate into the world by camera yaw"""
    Lc = sun_dir(az, elev)                       # camera frame: x right, z forward
    fh = norm([cam.f[0], 0, cam.f[2]]); rh = np.cross([0, 1, 0], fh)
    return norm(Lc[0] * rh + Lc[1] * np.array([0, 1, 0]) + Lc[2] * fh)

def _rnd(j, o, seed):
    return (np.sin(j * 12.9898 + o * 78.233 + seed * 3.137) * 43758.5453) % 1.0

# ------------------------------------------------------------------ backdrop (env map)
class Backdrop:
    """Sky, far mountains, the far bank and its spruce treeline -- all far enough away
    (100 m+) to be a function of ray DIRECTION only. So the water's reflection of them is
    automatically the mirror image about the horizon, as it is for distant things."""
    def __init__(self, st, cam, L, seed):
        self.L = L; self.seed = seed; self.st = st
        self.fh = norm([cam.f[0], 0, cam.f[2]])
        eye_h = cam.pos[1]
        self.eye_h = eye_h
        self.far = st.far
        self.e_water = np.arctan2(-eye_h, st.far)            # where the far shore meets the water
        self.e_bank = np.arctan2(2.5 - eye_h, st.far)        # top of the far bank
        self.cam = cam
        rh = np.cross([0, 1, 0], self.fh)
        self.sun_side = -dot(L, rh)                          # >0: sun on the left

    def angles(self, d):
        x = dot(d, np.cross([0, 1, 0], self.fh)); z = dot(d, self.fh)
        az = np.arctan2(x, z)
        el = np.arcsin(np.clip(d[..., 1], -1, 1))
        return az, el

    def mountains(self, az):
        s = self.seed
        a = az * 2.2
        base = np.radians(3.0) + np.radians(6.0) * np.clip(0.5 + 0.7 * fbm2(a, np.full_like(a, 3.1), s + 21, 4), 0, 1.4)
        ridge = 1 - np.abs(fbm2(a * 2.3, np.full_like(a, 7.7), s + 22, 3))
        return base + np.radians(3.0) * ridge ** 3

    def trees(self, az, el):
        """two rows of spruce spires standing on the far bank. Returns (mask, row, lit_side)"""
        mask = np.zeros(az.shape, bool); row = np.zeros(az.shape); side = np.zeros(az.shape)
        for r, (dens, hmul, off) in enumerate([(26.0, 0.75, 0.0), (18.0, 1.0, 0.5)]):   # back row, front row
            u = az * dens + off
            i0 = np.floor(u)
            for k in (-1, 0, 1):
                i = i0 + k
                cx = i + 0.5 + 0.5 * (_rnd(i, 3 + r, self.seed) - 0.5)
                present = _rnd(i, 2 + r, self.seed) > 0.12
                hgt = np.arctan2((10 + 12 * _rnd(i, 1 + r, self.seed)) * hmul, self.far)
                base = self.e_bank - np.radians(0.4)
                hrel = (el - base) / hgt                            # 0 at foot, 1 at tip
                half = 0.26 * (1 - hrel) ** 1.05                    # spire half-width, in tree units
                tier = (hrel * (7 + 4 * _rnd(i, 4, self.seed)) + _rnd(i, 5, self.seed)) % 1.0
                half = half * (0.62 + 0.38 * tier)                  # jagged branch tiers
                dx = u - cx
                inside = present & (hrel > 0) & (hrel < 1) & (np.abs(dx) < half * (dens / 26.0) * 1.9)
                mask |= inside
                row = np.where(inside, np.maximum(row, r), row)
                side = np.where(inside, np.sign(dx), side)
        return mask, row, side

    def sample(self, d):
        """-> (value in [0,1] 'global', material)"""
        az, el = self.angles(d)
        Lh = norm([self.L[0], 0, self.L[2]])
        dh = norm(np.stack([d[..., 0], np.zeros_like(el), d[..., 2]], -1))
        toward_sun = np.clip(dot(dh, Lh), -1, 1)
        e = np.clip(np.abs(el) / np.radians(35), 0, 1)
        v = 0.88 - 0.26 * e ** 0.8 + 0.07 * np.clip(toward_sun, 0, 1) ** 3 * (1 - e)
        mat = np.full(el.shape, SKY)
        # mountains
        em = self.mountains(az)
        m = el < em
        dA = 1e-3
        slope = (self.mountains(az + dA) - self.mountains(az - dA)) / (2 * dA)
        gully = fbm2(az * 30, el * 30, self.seed + 30, 3)
        face = np.clip(slope * 3 + gully * 0.9, -1, 1)          # >0: face turned to the left
        lit = face * np.sign(self.sun_side) > -0.05 if abs(self.sun_side) > 0.2 else gully > 0.0
        snow = (fbm2(az * 22, el * 22, self.seed + 31, 3) > 0.1) & (el > em - np.radians(3.5))
        vm = np.where(lit, 0.64, 0.44) + 0.08 * snow
        vm = vm * 0.75 + v * 0.25                               # aerial perspective: toward sky
        v = np.where(m, vm, v); mat = np.where(m, MTN, mat)
        # far bank (below the tree feet) -- a dark band down to the waterline
        fb = el < self.e_bank + np.radians(0.25) * np.sin(az * 40 + self.seed)
        v = np.where(fb, 0.2, v); mat = np.where(fb, TREES, mat)
        tm, row, side = self.trees(az, el)
        tv = np.where(row > 0.5, 0.12, 0.19)
        tv = tv + 0.06 * ((side * np.sign(self.sun_side)) < 0) * (abs(self.sun_side) > 0.2)
        v = np.where(tm, tv, v); mat = np.where(tm, TREES, mat)
        return v, mat

'''+s[b:]
# water limited to the tarn
s=s.replace("water = (tw < np.where(hit, t, np.inf)) & (tw < tmax * 1.5)",
 "hd = tw * np.sqrt(np.clip(1 - d[..., 1] ** 2, 0, 1))\n    water = (tw < np.where(hit, t, np.inf)) & (hd < st.far)")
s=s.replace("tmax = st.far + st.D + 80","tmax = st.far + st.D")
# land key
s=s.replace("""        vh, lith, lah = form_value(Nh, L, shh, aoh, albedo, sky=np.where(is_b, 0.35, 0.25))""","""        vh, lith, lah = form_value(Nh, L, shh, aoh, albedo, sky=np.where(is_b, 0.35, 0.25),
                                   key=np.where(is_b, -1.0, 0.42), key_gain=3.0)""")
open('scene.py','w').write(s)

v=open('valuemodel.py').read()
v=v.replace("""               ao_power=1.4, contact=0.3, sky=0.0):""","""               ao_power=1.4, contact=0.3, sky=0.0, key=-1.0, key_gain=3.0):""")
v=v.replace("""    eff = np.clip(ndl, 0, 1) * albedo""","""    eff = np.clip(ndl, 0, 1) * albedo
    # a KEYED material (grass): flat ground sits at `key` in the light family whatever the
    # sun height; slopes turned toward the light climb, slopes turned away fall.
    # (vertical blades catch a low sun that a flat surface would barely receive)
    keyed = np.asarray(key) >= 0
    eff = np.where(keyed, np.clip(key + key_gain * (ndl - L[1]), 0.0, 1.0), eff)""")
v=v.replace("""    w = np.clip((shadow - 0.35) / 0.3, 0, 1) * (ndl > 0)""","""    w = np.clip((shadow - 0.35) / 0.3, 0, 1) * ((ndl > 0) | (keyed & (eff > 0.02)))""")
open('valuemodel.py','w').write(v)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0

# output:
~/work/pixelart-studies/pedagogy/atelier/scene.py:87: RuntimeWarning: invalid value encountered in power
  half = 0.26 * (1 - hrel) ** 1.05                    # spire half-width, in tree units
medium_left_10m_s0

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("""        b_shore = -self.front + self.diag * a + wob
        left = b - b_shore""","""        b_shore = -self.front + self.diag * a + wob
        # the promontory's far side: from the boulder's right flank, receding back-left
        a_edge = self.bsize * 0.45 - self.recede * np.clip(b, 0, None) + 0.3 * s * np.sin(b / (1.1 * s) + w[4]) * np.clip(b / s, 0, 1)
        left = np.minimum((b - b_shore) * np.cos(np.arctan(self.diag)), (a_edge - a) / np.sqrt(1 + self.recede ** 2))""")
s=s.replace("""        return np.maximum(left * np.cos(np.arctan(self.diag)), near)""","""        return np.maximum(left, near)""")
s=s.replace("""        self.far = 100 + rng.uniform(0, 50)""","""        self.recede = rng.uniform(0.45, 0.9)
        self.far = 150 + rng.uniform(0, 50)""")
open('geometry.py','w').write(s)
s=open('scene.py').read()
s=s.replace("hgt = np.arctan2((10 + 12 * _rnd(i, 1 + r, self.seed)) * hmul, self.far)","hgt = np.arctan2((9 + 9 * _rnd(i, 1 + r, self.seed)) * hmul, self.far)")
s=s.replace("half = 0.26 * (1 - hrel) ** 1.05","half = 0.26 * np.clip(1 - hrel, 0, 1) ** 1.05")
s=s.replace("""    look = norm(dh * np.cos(np.radians(10)) - np.array([0, np.sin(np.radians(10)), 0]))   # 10 deg down""","""    look = norm(dh * np.cos(np.radians(PITCH)) - np.array([0, np.sin(np.radians(PITCH)), 0]))""")
s=s.replace("W, H = 320, 180","W, H = 320, 180\nPITCH = 9.0          # 'about 10 deg down': 9 leaves a strip of sky above the far trees")
open('scene.py','w').write(s)
v=open('valuemodel.py').read()
v=v.replace("light_gamma=0.85, core=0.9, bounce=1.25,","light_gamma=0.85, core=0.55, bounce=0.95,")
v=v.replace("np.exp(-((ndl + 0.10) / 0.16) ** 2)","np.exp(-((ndl + 0.07) / 0.10) ** 2)")
open('valuemodel.py','w').write(v)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 2>&1 | tail -1

# output:
medium_left_10m_s0

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
import pickle, numpy as np
from core import *
b=pickle.load(open('stage3/buf_medium_left_10m_s0.pkl','rb'))
sh=b['shadow']; mat=b['mat']; N=b['N']; L=b['L']
print('L',L)
land=mat==3
print('land frac shadowed', (sh[land]<0.5).mean())
im=(np.stack([sh]*3,-1)*255).astype(np.uint8)
from PIL import Image; Image.fromarray(im).resize((640,360),Image.NEAREST).save('stage3/dbg_shadow.png')
nd=np.clip((N*L).sum(-1),0,1); Image.fromarray((nd*255).astype(np.uint8)).resize((640,360),Image.NEAREST).save('stage3/dbg_ndl.png')
"


# output:
L [-0.98687596  0.15643447  0.04005111]
land frac shadowed 0.27278521351179097

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
import pickle, numpy as np
b=pickle.load(open('stage3/buf_medium_left_10m_s0.pkl','rb'))
mat=b['mat']; land=mat==3
print('v land', np.percentile(b['v'][land],[5,50,95]))
print('N land', b['N'][land].mean(0), 'ao', np.percentile(b['ao'][land],[5,50,95]), 'lit', b['lit'][land].mean())
print('depth', np.percentile(b['depth'][land],[5,50,95]))
"


# output:
v land [0.        0.        0.3319029]
N land [ 0.04448641  0.85954048 -0.20709762] ao [0.        0.        0.9049474] lit 0.7259400892288081
depth [ 9.80554989 11.14652249 15.6775344 ]

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace('''    def scene(self, p):
        """solid scene without the water''','''    def scene_ao(self, p):
        """un-scaled distances (vertical distance for the terrain) for occlusion sampling"""
        g = p[..., 1] - self.height(p[..., 0], p[..., 2])
        o = self.boulder.sdf(p)
        return np.minimum(g, o), np.where(o < g, 1, 0)

    def scene(self, p):
        """solid scene without the water''')
open('geometry.py','w').write(s)
s=open('scene.py').read()
s=s.replace("aoh = ambient_occ(sc, ph, Nh, scale=0.025 * st.bsize)","aoh = ambient_occ(st.scene_ao, ph, Nh, scale=0.025 * st.bsize)")
open('scene.py','w').write(s)
s=open('geometry.py').read()
s=s.replace("a_edge = self.bsize * 0.45 - self.recede","a_edge = self.bsize * 0.95 - self.recede")
s=s.replace("self.front = self.boulder.r[2] * 0.55","self.front = self.boulder.r[2] * 0.75")
s=s.replace("self.near = 0.42 * distance","self.near = 0.56 * distance")
open('geometry.py','w').write(s)
s=open('scene.py').read()
s=s.replace("""        gully = fbm2(az * 30, el * 30, self.seed + 30, 3)
        face = np.clip(slope * 3 + gully * 0.9, -1, 1)          # >0: face turned to the left""","""        # big facets: a low-frequency 'ridge field' whose sideways slope says which way
        # each face turns; only a few ridges, so faces read as planes, not as mottle
        rf = lambda A, E: fbm2(A * 9, E * 5 + A * 3, self.seed + 30, 2)
        gully = (rf(az + dA * 4, el) - rf(az - dA * 4, el)) / (8 * dA)
        face = np.clip(slope * 2 + gully * 0.08, -1, 1)          # >0: face turned to the left""")
s=s.replace("lit = face * np.sign(self.sun_side) > -0.05 if abs(self.sun_side) > 0.2 else gully > 0.0","lit = face * np.sign(self.sun_side) > 0.0 if abs(self.sun_side) > 0.2 else (face > 0.3)")
s=s.replace("snow = (fbm2(az * 22, el * 22, self.seed + 31, 3) > 0.1)","snow = (fbm2(az * 14, el * 10, self.seed + 31, 2) > 0.15)")
open('scene.py','w').write(s)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 2>&1 | tail -1

# output:
medium_left_10m_s0

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('scene.py').read()
s=s.replace("base = np.radians(3.0) + np.radians(6.0) * np.clip(0.5 + 0.7 * fbm2(a, np.full_like(a, 3.1), s + 21, 4), 0, 1.4)",
            "base = np.radians(1.2) + np.radians(3.6) * np.clip(0.5 + 0.8 * fbm2(a, np.full_like(a, 3.1), s + 21, 4), 0, 1.4)")
s=s.replace("return base + np.radians(3.0) * ridge ** 3","return base + np.radians(2.2) * ridge ** 3")
s=s.replace("""        rf = lambda A, E: fbm2(A * 9, E * 5 + A * 3, self.seed + 30, 2)
        gully = (rf(az + dA * 4, el) - rf(az - dA * 4, el)) / (8 * dA)
        face = np.clip(slope * 2 + gully * 0.08, -1, 1)          # >0: face turned to the left""","""        depth_below = (em - el) / np.radians(4.0)                # 0 at the crest, ~1 lower down
        # spurs run diagonally down and outward from the crests: shear the pattern by depth
        g = fbm2(az * 14 + depth_below * 2.2 * np.sign(slope + 1e-9), depth_below * 1.3, self.seed + 30, 2)
        face = np.clip(np.sign(slope) * 0.5 + 1.4 * g, -1, 1)  # >0: face turned to the left""")
open('scene.py','w').write(s)
s=open('stage3_setting.py').read()
s=s.replace("""    rgi = rgi * water_dark + 0.25""","""    # Fresnel: grazing water is a mirror; looking down into it, it goes dark
    cos_i = np.abs(buf["dwy"])                                  # sin of the depression angle
    R = 0.18 + 0.82 * (1 - cos_i) ** 5 * 0 + 0.82 * np.exp(-cos_i / 0.14)
    rgi = (rgi * water_dark + 0.25) * R + WATER_BODY * (1 - R)""")
s=s.replace("STEP_TO_G = {","WATER_BODY = 1.6          # index into G: the peaty dark of a tarn seen straight down\nSTEP_TO_G = {")
s=s.replace('''keep = {k: buf[k] for k in ["v",''','''keep = {k: buf[k] for k in ["dwy", "v",''')
open('stage3_setting.py','w').write(s)
s=open('scene.py').read()
s=s.replace("buf.update(v=v, mat=m,","buf.update(dwy=d[..., 1], v=v, mat=m,")
open('scene.py','w').write(s)
s=open('geometry.py').read()
s=s.replace("lip=0.22, slope=0.06, hummock=0.20, hummock_scale=1.7):","lip=0.22, slope=0.06, hummock=0.28, hummock_scale=3.2):")
s=s.replace("self.lobe_c = np.array([np.cos(ang) * self.r[0] * 0.45,","self.lobe_c = np.array([np.cos(ang) * self.r[0] * 0.32,")
open('geometry.py','w').write(s)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 2>&1 | tail -1

# output:
medium_left_10m_s0

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("""        b_shore = -self.front + self.diag * a + wob""","""        # ...and to the left it curls round toward the camera to join the near bank: a cove
        b_shore = -self.front + self.diag * a + wob - self.curl * np.clip(-a, 0, None) ** 2 / self.D""")
s=s.replace("""        self.recede = rng.uniform(0.45, 0.9)""","""        self.recede = rng.uniform(0.45, 0.9)
        self.curl = rng.uniform(0.9, 1.4)""")
open('geometry.py','w').write(s)
EOF
for l in left front right back; do (nice -n 10 python3 stage3_setting.py medium $l 10m 0 >/dev/null 2>&1 &); done; until [ $(pgrep -f stage3_setting | wc -l) -eq 0 ]; do sleep 3; done
python3 -c "
from PIL import Image
import sys; sys.path.insert(0,'.')
from core import grid, label
ims=[label(Image.open(f'stage3/setting_medium_{l}_10m_s0.png'),l) for l in ['left','front','right','back']]
grid(ims,2).save('stage3/v1_lights.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("self.curl = rng.uniform(0.9, 1.4)","self.curl = rng.uniform(2.4, 3.4)")
s=s.replace("""        left = np.minimum((b - b_shore) * np.cos(np.arctan(self.diag)), (a_edge - a) / np.sqrt(1 + self.recede ** 2))""","""        A = (b - b_shore) * np.cos(np.arctan(self.diag)); B = (a_edge - a) / np.sqrt(1 + self.recede ** 2)
        k = 0.35 * s                                                   # round the point of the promontory
        left = -k * np.logaddexp(-A / k, -B / k)""")
open('geometry.py','w').write(s)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 >/dev/null 2>&1; echo done

# output:
done

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from core import *
from scene import *
P=Params('medium','left','10m',0)
st=Setting(boulder_size=P.bsize, distance=P.D, seed=0)
cam=make_camera(st,P)
xs=np.linspace(-15,15,300); zs=np.linspace(-15,15,300)
X,Z=np.meshgrid(xs,zs)
s=st.shore(X,Z)
img=np.where(s>0,180,60).astype(np.uint8)
im=np.stack([img]*3,-1)[::-1].copy()
def px(x,z): return int((x+15)/30*299), int(299-(z+15)/30*299)
for P3,c in [((0,0,0),(255,0,0)),(tuple(cam.pos),(0,255,0))]:
    i,j=px(P3[0],P3[2]); im[max(0,j-2):j+3,max(0,i-2):i+3]=c
# view frustum
for sgn in [-1,1]:
    ang=np.radians(26)*sgn
    f=cam.f; r=cam.r
    d=f*np.cos(ang)+r*np.sin(ang)
    for t in np.linspace(0,25,200):
        p=cam.pos+d*t; i,j=px(p[0],p[2])
        if 0<=i<300 and 0<=j<300: im[j,i]=(255,255,0)
from PIL import Image; Image.fromarray(im).resize((600,600),Image.NEAREST).save('stage3/dbg_map.png')
print(cam.pos, st.fwd)
"


# output:
[-0.05971989  2.52326981 -9.99982168] [0.00597199 0.         0.99998217]

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("""               - 0.02 * across ** 2 / max(1, self.D * 0.1)""","""               + 0.55 * np.clip(-across, 0, None) - 0.15 * np.clip(across, 0, None)""")
s=s.replace("h_water = self.lip + s * 2.5 - 0.02","h_water = self.lip + s * 5.0 - 0.02")
s=s.replace("lip=0.22, slope=0.06,","lip=0.16, slope=0.06,")
open('geometry.py','w').write(s)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 >/dev/null 2>&1; echo done

# output:
done

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('stage3_setting.py').read()
s=s.replace("LAND:    [0.2, 1.0, 2.0, 4.0, 5.6, 6.8],","LAND:    [1.0, 2.1, 2.9, 4.3, 5.4, 6.5],")
s=s.replace("R = 0.18 + 0.82 * (1 - cos_i) ** 5 * 0 + 0.82 * np.exp(-cos_i / 0.14)","R = 0.35 + 0.65 * np.exp(-cos_i / 0.22)")
s=s.replace("WATER_BODY = 1.6 ","WATER_BODY = 2.4 ")
s=s.replace('''    to_image(q, GREY10, 3).save(f"{OUT}/setting_{name}.png")''','''    to_image(q, GREY10, 3).save(f"{OUT}/setting_{name}.png")
    notan = np.digitize(gi, [3.2, 5.6])                          # three-value massing
    to_image(notan, np.array([[45]*3, [128]*3, [215]*3]), 2).save(f"{OUT}/notan_{name}.png")''')
open('stage3_setting.py','w').write(s)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 >/dev/null 2>&1; echo done

# output:
done

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('stage3_setting.py').read()
s=s.replace("R = 0.35 + 0.65 * np.exp(-cos_i / 0.22)","R = 0.62 + 0.38 * np.exp(-cos_i / 0.2)")
s=s.replace("WATER_BODY = 2.4 ","WATER_BODY = 3.2 ")
s=s.replace("LAND:    [1.0, 2.1, 2.9, 4.3, 5.4, 6.5],","LAND:    [1.0, 2.0, 2.8, 3.9, 4.8, 5.8],")
open('stage3_setting.py','w').write(s)
s=open('scene.py').read()
s=s.replace('LIGHTS = {"left": 8, "front": 82, "right": 172, "back": 262}','LIGHTS = {"left": 22, "front": 90, "right": 158, "back": 270}')
open('scene.py','w').write(s)
EOF
for l in left front right back; do (nice -n 10 python3 stage3_setting.py medium $l 10m 0 >/dev/null 2>&1 &); done; until [ $(pgrep -f stage3_setting | wc -l) -eq 0 ]; do sleep 3; done
python3 -c "
from PIL import Image
import sys; sys.path.insert(0,'.')
from core import grid, label
ims=[label(Image.open(f'stage3/setting_medium_{l}_10m_s0.png'),l) for l in ['left','front','right','back']]
ims+=[label(Image.open(f'stage3/notan_medium_{l}_10m_s0.png').resize((960,540)),l+' notan') for l in ['left','back']]
grid(ims,2).save('stage3/v2_lights.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && for d in close 40m; do (nice -n 10 python3 stage3_setting.py medium left $d 0 >/dev/null 2>&1 &); done; (nice -n 10 python3 stage3_setting.py large left close 1 >/dev/null 2>&1 &); (nice -n 10 python3 stage3_setting.py small right 10m 2 >/dev/null 2>&1 &); until [ $(pgrep -f stage3_setting | wc -l) -eq 0 ]; do sleep 3; done
python3 -c "
from PIL import Image
import sys; sys.path.insert(0,'.')
from core import grid, label
n=['medium_left_close_s0','medium_left_40m_s0','large_left_close_s1','small_right_10m_s2']
grid([label(Image.open(f'stage3/setting_{x}.png'),x) for x in n],2).save('stage3/v2_dists.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
i=s.index("# ------------------------------------------------------------------ the setting")
s=s[:i]+'''class HullBoulder(Boulder):
    """Second attempt, after the first read as a pillow/loaf in close-up.
    A boulder is closer to a POTATO WITH FACETS than to an ellipsoid with planes cut off:
    take a dozen-odd points on a lopsided mass, wrap them in their convex hull (that gives
    20-30 facets of very UNEQUAL size -- a few big planes, many small), round the edges by
    a log-sum-exp blend (weathering), then add a gentle low-frequency swell and one or two
    shallow spalls (concave scoops where a shell of granite flaked off)."""
    def __init__(self, size=1.8, seed=0, center=(0, 0, 0), sink=0.38, n_pts=14,
                 round_k=None, lump=0.022, cracks=1, spalls=2):
        from scipy.spatial import ConvexHull
        rng = np.random.default_rng(seed)
        self.seed = seed; self.size = size
        s = size / 2
        self.r = np.array([s * rng.uniform(1.0, 1.18), s * rng.uniform(0.62, 0.8), s * rng.uniform(0.82, 1.0)])
        self.c = np.array(center, float) + np.array([0, self.r[1] * (1 - 2 * sink), 0])
        yaw = rng.uniform(0, 2 * np.pi)
        self.rot = np.array([[np.cos(yaw), 0, np.sin(yaw)], [0, 1, 0], [-np.sin(yaw), 0, np.cos(yaw)]])
        dirs = rng.normal(size=(n_pts, 3)); dirs /= np.linalg.norm(dirs, axis=1, keepdims=True)
        dirs[:, 1] = np.abs(dirs[:, 1]) * np.sign(dirs[:, 1] + 0.35)       # favour the upper half
        pts = dirs * self.r * rng.uniform(0.86, 1.06, (n_pts, 1))
        # a ring of points at the buried waist keeps the foot broad (it sits, it doesn't balance)
        ang = np.linspace(0, 2 * np.pi, 9)[:-1] + rng.uniform(0, 1)
        ring = np.stack([np.cos(ang) * self.r[0], -0.55 * self.r[1] * np.ones(8), np.sin(ang) * self.r[2]], 1) * rng.uniform(0.92, 1.05, (8, 1))
        pts = np.vstack([pts, ring, [[0, -self.r[1], 0]]])
        hull = ConvexHull(pts)
        eq = hull.equations                                  # n.x + d <= 0 inside
        # merge near-coplanar facets (ConvexHull triangulates)
        planes = []
        for n_, d_ in zip(eq[:, :3], eq[:, 3]):
            if not any(np.dot(n_, m) > 0.995 and abs(d_ - e) < 0.02 * size for m, e in planes):
                planes.append((n_, d_))
        self.ns = np.array([p[0] for p in planes]); self.offs = -np.array([p[1] for p in planes])
        self.k = round_k if round_k is not None else rng.uniform(0.045, 0.075) * size
        self.lump = lump * size
        self.cracks = []
        for i in range(cracks):
            th = rng.uniform(0, np.pi)
            n = norm([np.cos(th), rng.uniform(-0.3, 0.3), np.sin(th)])
            mid = rng.uniform(-0.4, 0.4, 3) * self.r
            self.cracks.append((n, float(mid @ n), rng.uniform(0.006, 0.01) * size, mid, 0.5 * s))
        self.spalls = []
        for i in range(spalls):
            dvec = norm(rng.normal(size=3) * np.array([1, 0.6, 1]) + np.array([0, 0.5, 0]))
            rad = rng.uniform(0.28, 0.45) * s
            self.spalls.append((dvec * (self.r * 1.0) + dvec * rad * 0.72, rad))

    def sdf(self, p):
        q = self.local(p)
        sh = q.shape
        q2 = q.reshape(-1, 3)
        dp = (q2 @ self.ns.T - self.offs) / self.k
        m = dp.max(1, keepdims=True)
        d = (self.k * (np.log(np.exp(dp - m).sum(1)) + m[:, 0])).reshape(sh[:-1])
        d = d - self.k * 0.9                                   # re-inflate what the blend ate
        for c, rad in self.spalls:
            d = smax(d, rad - np.linalg.norm(q - c, axis=-1), 0.12 * rad)
        if self.lump > 0:
            d = d + self.lump * fbm3(q / (self.size * 0.5), self.seed + 99, 3)
        for n, o, w, mid, half in self.cracks:
            along = np.linalg.norm((q - mid) - np.outer(q.reshape(-1, 3) @ n - o, n).reshape(q.shape), axis=-1)
            taper = np.clip(1 - along / half, 0, 1)
            band = np.abs(q @ n - o) - w * taper
            carve = np.maximum(np.maximum(band, -(d + 0.03 * self.size * taper)), -(taper - 0.05))
            d = np.maximum(d, -carve)
        return d

'''+s[i:]
s=s.replace("self.boulder = Boulder(size=boulder_size, seed=seed, center=(0.0, lip, 0.0))","self.boulder = HullBoulder(size=boulder_size, seed=seed, center=(0.0, lip, 0.0))")
open('geometry.py','w').write(s)
s=open('stage2_boulder.py').read()
s=s.replace("from geometry import Boulder","from geometry import Boulder, HullBoulder")
s=s.replace("""        b = Boulder(size=1.8, seed=seed)""","""        b = (HullBoulder if os.environ.get("HULL") else Boulder)(size=1.8, seed=seed)""")
s=s.replace('''        grid(ims, 4).save(f"{OUT}/boulder_seed{seed}.png")''','''        grid(ims, 4).save(f"{OUT}/{'hull' if os.environ.get('HULL') else 'boulder'}_seed{seed}.png")''')
open('stage2_boulder.py','w').write(s)
EOF
for s in 1 2 3 4; do (HULL=1 nice -n 10 python3 stage2_boulder.py $s > stage2/log$s.txt 2>&1 &) ; done; until [ $(pgrep -f stage2_boulder | wc -l) -eq 0 ]; do sleep 3; done; cat stage2/log*.txt; rm stage2/log*
python3 -c "
from PIL import Image
ims=[Image.open(f'stage2/hull_seed{s}.png') for s in [1,2,3,4]]
rows=[im.crop((0,im.height//2,im.width,im.height)) for im in ims]
out=Image.new('RGB',(rows[0].width,sum(r.height for r in rows)))
y=0
for r in rows: out.paste(r,(0,y)); y+=r.height
out.save('stage2/hull_all.png')"

# output:
ok
ok
ok
ok

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('scene.py').read()
s=s.replace("self.D = D if D is not None else max(3.2, 2.2 * self.bsize)","self.D = D if D is not None else max(3.6, 2.7 * self.bsize)")
open('scene.py','w').write(s)
s=open('stage3_setting.py').read()
s=s.replace("""    gi = np.where(solid, gi * (1 - 0.5 * far) + 6.5 * 0.5 * far, gi)""","""    gi = np.where(solid, gi * (1 - 0.5 * far) + 6.5 * 0.5 * far, gi)
    # distance also closes up the value range of the grass (its shadows lift toward its light)
    fl = np.clip((buf["depth"] - 8) / 40, 0, 0.6)
    gi = np.where(buf["mat"] == LAND, gi + (4.6 - gi) * fl * (gi < 4.6), gi)""")
open('stage3_setting.py','w').write(s)
EOF
J="medium left close 0|medium left 10m 0|medium left 40m 0|large back 10m 1"
echo "$J" | tr '|' '\n' | while read a b c d; do (nice -n 10 python3 stage3_setting.py $a $b $c $d >/dev/null 2>&1 &); done; sleep 2; until [ $(pgrep -f stage3_setting | wc -l) -eq 0 ]; do sleep 3; done
python3 -c "
from PIL import Image
import sys; sys.path.insert(0,'.')
from core import grid, label
n=['medium_left_close_s0','medium_left_10m_s0','medium_left_40m_s0','large_back_10m_s1']
grid([label(Image.open(f'stage3/setting_{x}.png'),x) for x in n],2).save('stage3/v3.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('geometry.py').read()
s=s.replace("- self.curl * np.clip(-a, 0, None) ** 2 / self.D","- self.curl * np.clip(-a, 0, None) ** 2 / (6.0 * self.bsize ** 0.5)")
s=s.replace("self.curl = rng.uniform(2.4, 3.4)","self.curl = rng.uniform(0.9, 1.4)")
open('geometry.py','w').write(s)
s=open('scene.py').read()
s=s.replace("""        face = np.clip(np.sign(slope) * 0.5 + 1.4 * g, -1, 1)  # >0: face turned to the left""","""        # which way a face turns is decided by the crest it hangs from; tracing down-and-out
        # from the crest makes the spur edges lean like real ridges instead of dropping plumb
        ss = np.sign(slope + 1e-9)
        az2 = az + ss * depth_below * np.radians(3.0)
        slope2 = (self.mountains(az2 + dA) - self.mountains(az2 - dA)) / (2 * dA)
        face = np.clip(np.sign(slope2) * 0.5 + 1.2 * g, -1, 1)  # >0: face turned to the left""")
# varied treeline: clumps, gaps, height swell
s=s.replace("""                present = _rnd(i, 2 + r, self.seed) > 0.12
                hgt = np.arctan2((9 + 9 * _rnd(i, 1 + r, self.seed)) * hmul, self.far)""","""                clump = 0.5 + 0.5 * vnoise2(i / dens * 7.0, np.full_like(i, 1.3 + r), self.seed + 40)
                present = _rnd(i, 2 + r, self.seed) < 0.35 + 0.75 * clump
                hgt = np.arctan2((6 + 10 * _rnd(i, 1 + r, self.seed) * (0.5 + clump)) * hmul, self.far)""")
open('scene.py','w').write(s)
EOF
J="medium left 10m 0|medium left 40m 0"
echo "$J" | tr '|' '\n' | while read a b c d; do (nice -n 10 python3 stage3_setting.py $a $b $c $d >/dev/null 2>&1 &); done; sleep 2; until [ $(pgrep -f stage3_setting | wc -l) -eq 0 ]; do sleep 3; done
python3 -c "
from PIL import Image
import sys; sys.path.insert(0,'.')
from core import grid, label
n=['medium_left_10m_s0','medium_left_40m_s0']
grid([label(Image.open(f'stage3/setting_{x}.png'),x) for x in n],1).save('stage3/v4.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
s=open('scene.py').read()
a=s.index("        depth_below = (em - el) / np.radians(4.0)")
b=s.index("        vm = np.where(lit, 0.64, 0.44) + 0.08 * snow")
s=s[:a]+'''        depth_below = (em - el) / np.radians(4.0)                # 0 at the crest, ~1 lower down
        # relief seen face-on: a ridged field stretched vertically, so its crests run DOWN
        # the mountain like spurs; the crest profile's own slope dominates near the top.
        def relief(A, E):
            u = A * 60 + 0.6 * np.sin(E * 90); v2 = E * 18
            r1 = 1 - np.abs(fbm2(u, v2, self.seed + 30, 3))
            r2 = 1 - np.abs(fbm2(u * 2.1 + 5, v2 * 2.1, self.seed + 33, 2))
            return r1 ** 2 + 0.35 * r2 ** 2
        dR = (relief(az + dA * 0.5, el) - relief(az - dA * 0.5, el)) / dA
        top = np.exp(-depth_below * 2.5)
        face = np.clip(np.sign(slope) * top * 1.2 + np.tanh(dR * 0.02) * (1 - 0.5 * top), -1, 1)
        lit = face * np.sign(self.sun_side) > 0.0 if abs(self.sun_side) > 0.2 else (np.abs(face) < 0.35)
        snow = (relief(az, el) > 0.55) & (depth_below < 0.9) & (fbm2(az * 20, el * 12, self.seed + 31, 2) > -0.1)
'''+s[b:]
open('scene.py','w').write(s)
EOF
nice -n 10 python3 stage3_setting.py medium left 10m 0 >/dev/null 2>&1; echo ok

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && sed -i '' 's/u = A \* 60 + 0.6 \* np.sin(E \* 90); v2 = E \* 18/u = A * 16 + 0.5 * np.sin(E * 40); v2 = E * 5.5/' scene.py && grep -n "u = A" scene.py && nice -n 10 python3 stage3_setting.py medium left 10m 0 >/dev/null 2>&1; python3 -c "
from PIL import Image; im=Image.open('stage3/setting_medium_left_10m_s0.png'); im.crop((0,0,960,220)).save('stage3/crop_mtn.png')"

# output:
119:            u = A * 16 + 0.5 * np.sin(E * 40); v2 = E * 5.5

cd ~/work/pixelart-studies/pedagogy/atelier && cat > paint.py <<'EOF'
"""The painter: buffers -> indexed pixel art.

Every material owns a RAMP: an ordered list of palette entries from dark to light.
Painting happens in 'step space' (a float position along the material's ramp); marks
(grass strokes, granite flecks, ripple lines, edge accents) are offsets in step space;
quantising mixes only two ADJACENT entries of ONE ramp, with a checker dither.
Stage 4 uses grey ramps (value only). Stage 5 swaps in colour ramps of the SAME values.
"""
from core import *
from scene import SKY, MTN, TREES, LAND, BOULDER, WATER
from geometry import fbm3, fbm2

# ------------------------------------------------------------------ ramps (values)
# luminance targets per ramp entry, dark -> light. The colour stage keeps these values.
RAMP_LUM = {
    "sky":   [0.60, 0.67, 0.74, 0.81, 0.87, 0.93],
    "mtn":   [0.34, 0.42, 0.52, 0.63, 0.75, 0.86],   # shadow rock .. lit snow
    "trees": [0.05, 0.10, 0.17, 0.26],
    "grass": [0.06, 0.12, 0.20, 0.30, 0.42, 0.54, 0.66],
    "rock":  [0.07, 0.15, 0.25, 0.48, 0.66, 0.84],
    "water": [0.22, 0.30, 0.40, 0.50, 0.60, 0.70],   # the tarn itself: reflected sky, darkened
}
RAMP_ORDER = ["sky", "mtn", "trees", "grass", "rock", "water"]

def palette_layout(ramps=RAMP_LUM):
    off, k = {}, 0
    for r in RAMP_ORDER:
        off[r] = k; k += len(ramps[r])
    return off, k

def grey_palette(ramps=RAMP_LUM):
    cols = []
    for r in RAMP_ORDER:
        for l in ramps[r]:
            g = round(l * 255); cols.append((g, g, g))
    return np.array(cols)

# ------------------------------------------------------------------ helpers
def hash2(x, y, s):
    h = (x * 73856093) ^ (y * 19349663) ^ (s * 83492791)
    h = (h ^ (h >> 13)) * 1274126177
    return ((h ^ (h >> 16)) & 0xFFFF) / 65535.0

def shift(a, dy, dx, fill):
    out = np.full_like(a, fill)
    H, W = a.shape[:2]
    ys, yd = (slice(0, H - dy), slice(dy, H)) if dy >= 0 else (slice(-dy, H), slice(0, H + dy))
    xs, xd = (slice(0, W - dx), slice(dx, W)) if dx >= 0 else (slice(-dx, W), slice(0, W + dx))
    out[yd, xd] = a[ys, xs]
    return out

class Painter:
    def __init__(self, buf, seed=0, ramps=RAMP_LUM, marks=True):
        self.b = buf; self.seed = seed; self.ramps = ramps; self.marks = marks
        self.H, self.W = buf["mat"].shape
        self.off, self.n = palette_layout(ramps)
        yy, xx = np.mgrid[0:self.H, 0:self.W]
        self.yy, self.xx = yy, xx

    # -------------------------------------------------------------- base steps
    def material_steps(self, v, mat):
        """value-model output -> (ramp name per pixel as int code, step float)"""
        H, W = mat.shape
        code = np.full((H, W), -1); step = np.zeros((H, W))
        R = {n: i for i, n in enumerate(RAMP_ORDER)}
        m = mat == SKY;     code[m] = R["sky"];   step[m] = (v[m] - 0.60) / (0.95 - 0.60) * 5
        m = mat == MTN;     code[m] = R["mtn"];   step[m] = (v[m] - 0.36) / (0.80 - 0.36) * 5
        m = mat == TREES;   code[m] = R["trees"]; step[m] = (v[m] - 0.10) / (0.28 - 0.10) * 3
        m = mat == LAND;    code[m] = R["grass"]; step[m] = np.interp(v[m] * 5, np.arange(6), [0.3, 1.2, 2.1, 3.6, 4.7, 5.8])
        m = mat == BOULDER; code[m] = R["rock"];  step[m] = v[m] * 5
        return code, step

    def base(self):
        b = self.b
        code, step = self.material_steps(b["v"], b["mat"])
        # aerial perspective: distant land & rock lose contrast toward a mid step
        d = np.clip(b["depth"], 0, 1e3)
        far = 1 - np.exp(-d / 70.0)
        R = {n: i for i, n in enumerate(RAMP_ORDER)}
        g = code == R["grass"]; r = code == R["rock"]
        step[g] = step[g] + (3.6 - step[g]) * 0.55 * far[g]
        step[r] = step[r] + (3.2 - step[r]) * 0.55 * far[r]
        fl = np.clip((d - 8) / 40, 0, 0.55)
        step[g] = np.where(step[g] < 3.6, step[g] + (3.6 - step[g]) * fl[g], step[g])
        # water: the reflected thing, darker, broken by ripples (painted in marks)
        wcode, wstep = self.water_steps()
        w = b["water"]
        code[w] = wcode[w]; step[w] = wstep[w]
        self.code, self.step = code, step
        return code, step

    def water_steps(self, rv=None, rmat=None):
        b = self.b
        rv = b["rv"] if rv is None else rv
        rmat = b["rmat"] if rmat is None else rmat
        code, step = self.material_steps(rv, rmat)
        R = {n: i for i, n in enumerate(RAMP_ORDER)}
        # Fresnel: grazing water mirrors, water seen from above shows its dark body
        cos_i = np.abs(b["dwy"])
        dark = 1.6 * (1 - np.exp(-cos_i / 0.16))                     # steps of darkening
        # sky & mountains reflect into the WATER ramp (by value), the rest into their own
        sky_like = np.isin(rmat, [SKY, MTN])
        lum = np.zeros_like(step)
        for name in ["sky", "mtn"]:
            m = code == R[name]
            lum[m] = np.interp(step[m], np.arange(len(self.ramps[name])), self.ramps[name])
        wl = np.array(self.ramps["water"])
        wstep = np.interp(lum * 0.92 - 0.02, wl, np.arange(len(wl))) - dark
        code = np.where(sky_like, R["water"], code)
        step = np.where(sky_like, wstep, step - 0.9 - dark * 0.8)
        return code, step

    # -------------------------------------------------------------- marks
    def water_marks(self):
        """broken reflection: rows wobble sideways a little (more toward the viewer), and
        thin horizontal ripple lines lie on the surface -- longer near, shorter far."""
        b = self.b; H, W = self.H, self.W
        w = b["water"]
        if not w.any():
            return
        near = np.clip(np.abs(b["dwy"]) / 0.35, 0, 1)
        # wobble: sample the reflection from a neighbouring column
        amp = np.round(np.sin(self.yy * 0.9 + self.seed) * np.sin(self.yy * 0.23 + 1.7) * (0.6 + 2.0 * near)).astype(int)
        xs = np.clip(self.xx + amp, 0, W - 1)
        src_w = w[self.yy, xs]
        rv2 = np.where(src_w, b["rv"][self.yy, xs], b["rv"])
        rm2 = np.where(src_w, b["rmat"][self.yy, xs], b["rmat"])
        c2, s2 = self.water_steps(rv2, rm2)
        self.code[w] = c2[w]; self.step[w] = s2[w]
        # ripple lines: horizontal dashes of +/- a step
        seglen = np.clip(3 + 22 * near, 3, 26)
        seg = np.floor((self.xx + hash2(self.yy, 7, self.seed) * 50) / seglen).astype(int)
        r = hash2(seg, self.yy, self.seed + 3)
        lite = (r > 0.86) & w
        dark = (r < 0.10) & w
        self.step[lite] += 0.75
        self.step[dark] -= 0.75
        # a dark line where the water meets a bank above it (the turf's shadowed underside)
        up = shift(b["mat"], 1, 0, -1)
        edge = w & np.isin(up, [LAND, BOULDER]) & (shift(b["depth"], 1, 0, 0) < b["depth"] + 0.5)
        self.step[edge] -= 1.2

    def rock_marks(self):
        b = self.b
        rk = b["mat"] == BOULDER
        if not rk.any():
            return
        N = b["N"]; L = b["L"]; p = b["p"]; H, W = self.H, self.W
        depth = b["depth"]
        pix_m = depth * 2 * np.tan(np.radians(26)) / W           # metres per pixel
        size_px = self.b_size / np.maximum(pix_m, 1e-6)
        detail = np.clip((size_px - 12) / 50, 0, 1)              # how much texture the rock can carry
        # 1) edges between planes: convex ridges catch the light, creases go dark
        conv = np.zeros((H, W)); cre = np.zeros((H, W))
        for dy, dx in [(0, 1), (1, 0), (0, -1), (-1, 0)]:
            Nn = np.stack([shift(N[..., i], dy, dx, 0) for i in range(3)], -1)
            pn = np.stack([shift(p[..., i], dy, dx, 0) for i in range(3)], -1)
            same = shift(rk.astype(int), dy, dx, 0) == 1
            c = np.sum((N - Nn) * (p - pn), -1) / np.maximum(pix_m, 1e-6) ** 2 * pix_m
            ang = 1 - np.sum(N * Nn, -1)
            conv += np.where(same & (c < 0), ang, 0)             # neighbour turns away: ridge
            cre += np.where(same & (c > 0), ang, 0)
        ndl = np.sum(N * L, -1)
        ridge = rk & (conv > 0.05) & (ndl > 0.05)
        crease = rk & (cre > 0.12)
        self.step[ridge] += 0.8
        self.step[crease] -= 0.9
        # 2) silhouette: a dark accent where the form turns away from the light; a light rim
        #    where the light comes from behind and grazes the edge
        sil = np.zeros((H, W), bool)
        for dy, dx in [(0, 1), (1, 0), (0, -1), (-1, 0)]:
            nb_m = shift(b["mat"], dy, dx, -1); nb_d = shift(depth, dy, dx, 1e4)
            sil |= rk & ((nb_m != BOULDER) & (nb_d > depth + 0.3))
        back = L[2] * 0 + np.sum(L * self.cam_f) > 0.3             # sun behind the subject
        if back:
            rim = sil & (ndl > -0.45)
            self.step[rim] = np.maximum(self.step[rim], 3.4)
            rim2 = rk & ~sil & (ndl > -0.1) & shift(sil, 0, 1, False)
            self.step[rim2] = np.maximum(self.step[rim2], 2.6)
        else:
            dark_edge = sil & (ndl < 0.05)
            self.step[dark_edge] -= 0.8
        # 3) weathering mottle (sticks to the stone: object space), lichen patches on top planes
        q = self.boulder.local(p[rk])
        mott = fbm3(q / (self.b_size * 0.22), self.seed + 11, 3)
        self.step[rk] += 0.45 * mott * detail[rk]
        # 4) granite: salt-and-pepper flecks, sparse, 1-2 px clusters, quieter in shadow
        hsh = hash2(self.xx, self.yy, self.seed + 5)
        hsh2 = hash2(self.xx // 2, self.yy, self.seed + 6)
        lit = self.step > 2.6
        dk = rk & (np.where(lit, hsh < 0.09, hsh < 0.035)) & (detail > 0.15)
        lt = rk & (np.where(lit, hsh > 0.95, hsh > 0.985)) & (detail > 0.15)
        pair = rk & (hsh2 < 0.03) & lit & (detail > 0.4)
        self.fleck_dark = dk | pair; self.fleck_lite = lt

    def grass_marks(self):
        """vertical blade strokes: lengths from the depth, values from the light."""
        b = self.b; H, W = self.H, self.W
        gr = b["mat"] == LAND
        depth = b["depth"]
        pix_m = depth * 2 * np.tan(np.radians(26)) / W
        blade = np.clip(0.28 / np.maximum(pix_m, 1e-4), 1, 9)       # a 28 cm blade, in px
        rng = np.random.default_rng(self.seed + 77)
        off = np.zeros((H, W))
        for x in range(W):
            y = 0
            while y < H:
                L = blade[y, x]
                n = max(1, int(round(L * rng.uniform(0.5, 1.3))))
                r = rng.random()
                o = -1.0 if r < 0.38 else (0.9 if r > 0.82 else 0.0)
                off[y:y + n, x] = o
                y += n
        lit = self.step > 3.3
        o2 = np.where(lit, off, np.where(off > 0, off * 0.5, off * 0.8))
        self.step[gr] += o2[gr] * np.clip(blade[gr] / 2.5, 0.45, 1.0)
        # tufts: blades poke up over whatever is behind the top edge of a grass mass
        self.tufts = []
        top = gr & (shift(b["mat"], 1, 0, -1) != LAND)
        ys, xs = np.nonzero(top)
        code_g = RAMP_ORDER.index("grass")
        for y, x in zip(ys, xs):
            L = blade[y, x]
            if rng.random() > 0.55 or L < 1.5:
                continue
            n = int(round(L * rng.uniform(0.3, 0.9)))
            s0 = self.step[y, x] + (0.8 if rng.random() < 0.4 else -0.4)
            for k in range(1, n + 1):
                yy = y - k
                if yy < 0 or depth[yy, x] < depth[y, x] + 0.1:
                    break
                self.tufts.append((yy, x, code_g, s0))
        # grass round the boulder's foot overlaps the stone
        foot = (b["mat"] == BOULDER) & (shift(b["mat"], -1, 0, -1) == LAND)
        ys, xs = np.nonzero(foot)
        for y, x in zip(ys, xs):
            L = max(2, blade[min(y + 1, H - 1), x])
            n = int(round(L * rng.uniform(0.4, 1.2)))
            s0 = self.step[min(y + 1, H - 1), x] + (0.7 if rng.random() < 0.35 else -0.2)
            for k in range(0, n):
                yy = y - k
                if yy < 0 or b["mat"][yy, x] != BOULDER:
                    break
                self.tufts.append((yy, x, code_g, s0))

    def tree_marks(self):
        b = self.b
        t = b["mat"] == TREES
        # lit edges of the spires on the sun side: one step up where the left/right
        # neighbour is sky (the sun catches the tips of the branches)
        side = self.sun_side
        if abs(side) > 0.2:
            dx = -1 if side > 0 else 1
            nb = shift(b["mat"], 0, -dx, -1)
            edge = t & np.isin(nb, [SKY, MTN])
            self.step[edge] += 1.2
        hsh = hash2(self.xx, self.yy, self.seed + 9)
        self.step[t & (hsh > 0.93)] += 0.8

    # -------------------------------------------------------------- finish
    def quantize(self):
        idx = np.zeros((self.H, self.W), int)
        for i, name in enumerate(RAMP_ORDER):
            m = self.code == i
            if not m.any():
                continue
            n = len(self.ramps[name])
            v = np.clip(self.step / (n - 1), 0, 1)
            q = quantize(np.where(m, v, 0), n, kind="check", strength=0.55, flat_snap=0.1)
            idx[m] = q[m] + self.off[name]
        return idx

    def paint(self, boulder, cam_f, sun_side, b_size):
        self.boulder, self.cam_f, self.sun_side, self.b_size = boulder, cam_f, sun_side, b_size
        self.base()
        self.fleck_dark = self.fleck_lite = None
        self.tufts = []
        if self.marks:
            self.water_marks(); self.rock_marks(); self.grass_marks(); self.tree_marks()
        idx = self.quantize()
        idx = remove_orphans(idx)
        # marks that are MEANT to be single pixels go on after the clean-up
        if self.fleck_dark is not None:
            o = self.off["rock"]; n = len(self.ramps["rock"])
            rk = (idx >= o) & (idx < o + n)
            idx[self.fleck_dark & rk] = np.maximum(idx[self.fleck_dark & rk] - 1, o)
            idx[self.fleck_lite & rk] = np.minimum(idx[self.fleck_lite & rk] + 1, o + n - 1)
        T = threshold_map(self.H, self.W, "check")
        for (y, x, c, s) in self.tufts:
            name = RAMP_ORDER[c]; n = len(self.ramps[name])
            st = np.clip(s, 0, n - 1)
            lo = int(np.floor(st)); fr = st - lo
            idx[y, x] = self.off[name] + min(lo + (fr > T[y, x]), n - 1)
        return idx
EOF
cat > stage4_material.py <<'EOF'
"""Stage 4 -- material, still in greyscale (value only).
granite: planes kept, ridges catch light, creases & shadow-side silhouette darken, object-space
mottle, salt-and-pepper flecks. grass: vertical blade strokes scaled by distance, tufts over edges.
water: broken reflection + ripple dashes + dark line under the bank. trees: sun-side lit edges.
"""
from core import *
import os, sys
from scene import *
from paint import Painter, grey_palette
from backdrop_util import sun_side_of

OUT = "stage4"; os.makedirs(OUT, exist_ok=True)

def paint_scene(P, ramps=None, marks=True, palette=None):
    from paint import RAMP_LUM
    buf = render_buffers(P)
    pt = Painter(buf, seed=P.seed, ramps=ramps or RAMP_LUM, marks=marks)
    idx = pt.paint(buf["st"].boulder, buf["cam"].f, sun_side_of(buf), P.bsize)
    return idx, buf

if __name__ == "__main__":
    a = sys.argv[1:]
    P = Params(a[0], a[1], a[2], int(a[3]))
    name = f"{a[0]}_{a[1]}_{a[2]}_s{a[3]}"
    idx, buf = paint_scene(P, marks=True)
    to_image(idx, grey_palette(), 3).save(f"{OUT}/mat_{name}.png")
    print(name)
EOF
cat > backdrop_util.py <<'EOF'
from core import dot, norm
import numpy as np
def sun_side_of(buf):
    cam = buf["cam"]; L = buf["L"]
    fh = norm([cam.f[0], 0, cam.f[2]]); rh = np.cross([0, 1, 0], fh)
    return float(-dot(L, rh))     # >0: sun on the left
EOF
time nice -n 10 python3 stage4_material.py medium left close 0

# output:
medium_left_close_s0
nice -n 10 python3 stage4_material.py medium left close 0  29.42s user 14.70s system 74% cpu 59.285 total

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from PIL import Image; import numpy as np
im=Image.open('../shared/ferrari_alpine_tarn_dusk.png').convert('RGB')
a=np.array(im)[::2,::2]   # native ~634x469
nat=Image.fromarray(a); nat.save('ref/native.png'); print(nat.size)
for name,box in [('edge_island',(255,255,420,300)),('edge_mtn',(200,70,330,140)),('edge_tree',(20,20,100,110)),('edge_snow',(290,90,380,160)),('edge_bank',(140,205,260,260))]:
    c=nat.crop(box); c.resize((c.width*6,c.height*6),Image.NEAREST).save('ref/'+name+'.png')
"; python3 -c "import skimage; print('skimage', skimage.__version__)"

# output:
(634, 469)
skimage 0.25.2

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from PIL import Image; import numpy as np
a=np.array(Image.open('ref/native.png')).astype(int)
# sky: pinkish (r>190, g>140, b>140) in the mountain region rows 60-150, cols 170-300
reg=a[60:160,170:300]
sky=(reg[...,0]>195)&(reg[...,1]>135)&(reg[...,2]>135)&(reg[...,0]-reg[...,2]>15)
top=np.array([np.argmax(~sky[:,x]) for x in range(sky.shape[1])])
print('profile', top.tolist())
d=np.diff(top)
# run lengths of constant y
runs=[];c=1
for v in d:
    if v==0: c+=1
    else: runs.append(c); c=1
print('horizontal run lengths', runs)
print('vertical jumps', [int(v) for v in d if v!=0])
"

# output:
profile [64, 64, 64, 63, 63, 62, 62, 61, 61, 61, 60, 60, 59, 59, 60, 60, 61, 62, 62, 61, 61, 54, 54, 48, 54, 57, 60, 60, 60, 60, 60, 60, 59, 59, 58, 58, 57, 56, 55, 55, 54, 53, 52, 51, 50, 49, 49, 48, 47, 46, 45, 45, 44, 45, 46, 46, 47, 48, 49, 49, 50, 50, 49, 49, 48, 48, 47, 47, 47, 46, 46, 47, 47, 48, 48, 49, 49, 48, 47, 45, 44, 44, 44, 43, 43, 43, 43, 43, 43, 42, 42, 42, 42, 42, 41, 41, 41, 41, 40, 40, 40, 40, 39, 39, 39, 39, 40, 40, 39, 36, 34, 34, 35, 35, 35, 34, 34, 33, 33, 32, 32, 31, 31, 30, 29, 28, 27, 26, 26, 25]
horizontal run lengths [3, 2, 2, 3, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 6, 2, 2, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 2, 3, 2, 2, 2, 2, 1, 1, 1, 3, 6, 5, 4, 4, 4, 2, 1, 1, 2, 3, 2, 2, 2, 2, 1, 1, 1, 1, 2]
vertical jumps [-1, -1, -1, -1, -1, 1, 1, 1, -1, -7, -6, 6, 3, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -2, -1, -1, -1, -1, -1, -1, 1, -1, -3, -2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]

cd ~/work/pixelart-studies/pedagogy/atelier && cat > lines.py <<'EOF'
"""Line craft: contours drawn in clean pixel runs.

Measured on Ferrari's mountain skyline (native pixels): vertical steps are almost always
1 px; horizontal run lengths hold steady along a stretch (2,2,2,2 = a 2:1 line; 1,1,1 = 45
deg; 4,4,4) and change gradually round a curve. Ray-marched edges are noisy instead, with
runs like 1,3,1,2 -- 'jaggies'. So every major contour is re-drawn:
  smooth the region's outline -> simplify to a polygon -> snap each edge to a clean ratio
  (flat, 4:1, 3:1, 2:1, 1:1, 1:2, 1:3, 1:4, vertical) -> re-rasterise -> re-fill the pixels
  that changed side from their nearest neighbour on the correct side.
"""
import numpy as np
from PIL import Image, ImageDraw
from scipy import ndimage
from skimage import measure

NICE = [(1, 0), (4, 1), (3, 1), (2, 1), (1, 1), (1, 2), (1, 3), (1, 4), (0, 1)]
DIRS = []
for dx, dy in NICE:
    for sx in (1, -1):
        for sy in (1, -1):
            d = (dx * sx, dy * sy)
            if d not in DIRS:
                DIRS.append(d)
DIRS = np.array(DIRS, float)
DANG = np.arctan2(DIRS[:, 1], DIRS[:, 0])

def snap_polygon(pts, min_len=3.0):
    """walk the polygon; each edge becomes m * (a nice direction), aiming at the original
    next vertex so the drift never accumulates."""
    out = [np.round(pts[0])]
    cur = out[0].copy()
    for i in range(1, len(pts)):
        tgt = pts[i]
        e = tgt - cur
        L = np.hypot(*e)
        if L < min_len:
            cur = np.round(tgt); out.append(cur.copy()); continue
        ang = np.arctan2(e[1], e[0])
        # the two nicest candidate directions around this angle; take the better landing
        diff = np.abs((DANG - ang + np.pi) % (2 * np.pi) - np.pi)
        best = None
        for j in np.argsort(diff)[:2]:
            d = DIRS[j]
            m = max(1, round(np.dot(e, d) / np.dot(d, d)))
            p = cur + m * d
            err = np.hypot(*(p - tgt)) + 0.35 * diff[j] * L
            if best is None or err < best[0]:
                best = (err, p)
        cur = best[1]; out.append(cur.copy())
    return np.array(out)

def clean_mask(mask, sigma=1.0, tol=0.75, min_area=6, snap=True):
    """-> new boolean mask with the same regions, outlines drawn in clean runs."""
    H, W = mask.shape
    pad = 3
    m = np.pad(mask.astype(float), pad, mode="edge")
    sm = ndimage.gaussian_filter(m, sigma)
    out = Image.new("L", (W + 2 * pad, H + 2 * pad), 0)
    dr = ImageDraw.Draw(out)
    lab, n = ndimage.label(mask)
    # outer boundaries fill, holes cut (contours come back oriented; use area sign)
    cs = measure.find_contours(sm, 0.5)
    polys = []
    for c in cs:
        if len(c) < 4:
            continue
        xy = c[:, ::-1]                                     # (x, y)
        poly = measure.approximate_polygon(xy, tolerance=tol)
        if len(poly) < 3:
            continue
        if snap:
            poly = snap_polygon(poly)
        area = 0.5 * np.sum(poly[:-1, 0] * poly[1:, 1] - poly[1:, 0] * poly[:-1, 1])
        if abs(area) < min_area:
            continue
        # decide fill vs hole by sampling the smoothed mask just inside the polygon
        polys.append((abs(area), poly))
    polys.sort(key=lambda t: -t[0])
    for _, poly in polys:
        img = Image.new("L", out.size, 0)
        ImageDraw.Draw(img).polygon([tuple(p) for p in poly], fill=1, outline=1)
        a = np.array(img, bool)
        inside_val = sm[a].mean() if a.any() else 0
        if inside_val >= 0.5:
            out = Image.fromarray(np.maximum(np.array(out), a.astype(np.uint8)))
        else:
            o = np.array(out); o[a & (sm < 0.5)] = 0; out = Image.fromarray(o)
    res = np.array(out, bool)[pad:-pad, pad:-pad]
    return res

def reassign(buf, keys, new_mask, old_mask):
    """pixels that joined the region copy their data from the nearest old-region pixel;
    pixels that left copy from the nearest pixel outside the old region."""
    gained = new_mask & ~old_mask
    lost = old_mask & ~new_mask
    if gained.any():
        _, (iy, ix) = ndimage.distance_transform_edt(~old_mask, return_indices=True)
        for k in keys:
            a = buf[k]
            a[gained] = a[iy[gained], ix[gained]]
    if lost.any():
        _, (iy, ix) = ndimage.distance_transform_edt(old_mask, return_indices=True)
        for k in keys:
            a = buf[k]
            a[lost] = a[iy[lost], ix[lost]]
    return gained.sum(), lost.sum()

# ------------------------------------------------------------------ spruce sprites
def spruce(h, rng, lean=0):
    """A spruce as pixels, built the way Ferrari's are: a 1-px leader with ticks at the tip,
    then tiers of drooping branch strokes (a flat run that steps down once at its end --
    clean 2:1 / 3:1 runs), their length growing down the tree with some raggedness,
    around a solid core. Returns a list of (dy, dx, kind): kind 0 core, 1 branch, 2 tip.
    dy counts upward from the foot (0) to the tip (h-1)."""
    px = []
    tip = max(2, int(h * 0.14))
    half_max = max(1.5, h * rng.uniform(0.17, 0.24))
    for dy in range(h):
        px.append((dy, 0, 2 if dy > h - tip else 0))
    # ticks on the leader
    for dy in range(h - tip, h - 1, 2):
        s = 1 if (dy // 2) % 2 else -1
        px.append((dy, s, 2))
    dy = h - tip
    while dy > 0:
        env = half_max * ((h - dy) / h) ** 0.95
        for s in (-1, 1):
            L = int(round(env * rng.uniform(0.55, 1.15)))
            if L < 1:
                continue
            droop_at = max(1, int(round(L * rng.uniform(0.55, 0.8))))
            for k in range(1, L + 1):
                ddy = 0 if k <= droop_at else -1
                px.append((dy + ddy, s * k, 1))
            if L >= 3 and rng.random() < 0.5:                  # a second, shorter stroke below
                for k in range(1, L - 1):
                    px.append((dy - 1, s * k, 1))
        # core fill between tiers
        core = int(env * 0.45)
        for c in range(-core, core + 1):
            px.append((dy, c, 0))
            px.append((dy - 1, c, 0))
        dy -= rng.choice([2, 2, 3])
    return px
EOF
echo ok

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && cat > drawing.py <<'EOF'
"""Drawing pass over the rendered buffers: clean contours and hand-built trees.
render_buffers() -> drawing.finish(buf) -> Painter."""
import numpy as np
from core import norm, dot
from scene import SKY, MTN, TREES, LAND, BOULDER, WATER, _rnd
from lines import clean_mask, reassign, spruce

KEYS = ["v", "mat", "depth", "p", "N", "shadow", "ao", "lit", "lamt", "water", "rv", "rmat",
        "rdist", "dwy", "tw"]
RKEYS = ["rv", "rmat", "rdist"]

def clean_region(buf, mask, keys, within=None, sigma=1.0, tol=0.75):
    new = clean_mask(mask, sigma=sigma, tol=tol)
    if within is not None:
        new &= within
    return reassign(buf, keys, new, mask)

def horizon_row(cam):
    far = cam.pos + norm([cam.f[0], 0, cam.f[2]]) * 1e5
    x, y, z = cam.project(far)
    return float(y)

def elev_to_row(cam, e, yh):
    k = (cam.H / 2) / (cam.tanh * cam.H / cam.W)
    # small-angle is fine at these elevations; exact for the centre column
    pitch = np.arcsin(-cam.f[1])
    return cam.H / 2 - np.tan(e + pitch) * k

def finish(buf, trees=True, lines=True):
    cam = buf["cam"]; st = buf["st"]; P = buf["P"]
    mat = buf["mat"]
    if lines:
        # the bank against the water, the boulder against everything, its reflection
        clean_region(buf, (buf["mat"] == LAND), KEYS, sigma=1.1, tol=0.8)
        clean_region(buf, (buf["mat"] == BOULDER), KEYS, sigma=0.9, tol=0.7)
        w = buf["water"]
        clean_region(buf, w & (buf["rmat"] == BOULDER), RKEYS, within=w, sigma=0.9, tol=0.7)
        clean_region(buf, w & (buf["rmat"] == LAND), RKEYS, within=w, sigma=1.1, tol=0.8)
        clean_region(buf, (buf["mat"] == MTN), KEYS, sigma=0.8, tol=0.6)
        clean_region(buf, w & (buf["rmat"] == MTN), RKEYS, within=w, sigma=0.8, tol=0.6)
    if trees:
        draw_trees(buf)
    return buf

def draw_trees(buf):
    cam = buf["cam"]; st = buf["st"]; P = buf["P"]; L = buf["L"]
    H, W = buf["mat"].shape
    yh = horizon_row(cam)
    eye_h = cam.pos[1]
    e_bank = np.arctan2(2.5 - eye_h, st.far)
    rad_per_px = 2 * cam.tanh / W
    fh = norm([cam.f[0], 0, cam.f[2]]); rh = np.cross([0, 1, 0], fh)
    sun_side = -dot(L, rh)
    back = dot(L, fh) > 0.3
    rng = np.random.default_rng(P.seed + 555)
    mat, v = buf["mat"], buf["v"]
    water, rmat, rv = buf["water"], buf["rmat"], buf["rv"]
    # two rows of spruce: the back row paler (air), the front row darker and taller
    for row, (spacing, hmul, val) in enumerate([(9.0, 0.8, 0.20), (13.0, 1.0, 0.13)]):
        x = rng.uniform(-8, 0)
        while x < W + 8:
            clump = 0.5 + 0.5 * np.sin(x * 0.045 + P.seed * 1.7 + row) * np.sin(x * 0.017 + 2.1)
            if rng.random() < 0.25 + 0.7 * clump:
                h_m = (7 + 11 * rng.random()) * hmul * (0.7 + 0.5 * clump)
                h = int(round(np.arctan2(h_m, st.far) / rad_per_px))
                base_e = e_bank + np.radians(0.25) * np.sin(x * 0.13 + P.seed)
                yb = int(round(elev_to_row(cam, base_e, yh))) + rng.integers(0, 2)
                xi = int(round(x))
                for dy, dx, kind in spruce(max(h, 4), rng):
                    yy, xx = yb - dy, xi + dx
                    if not (0 <= xx < W):
                        continue
                    lit = (abs(sun_side) > 0.2 and kind == 1 and np.sign(dx) == -np.sign(sun_side)
                           and not back)
                    tv = val + (0.07 if lit else 0.0) + (0.03 if kind == 2 else 0.0)
                    if 0 <= yy < H and mat[yy, xx] in (SKY, MTN, TREES) and not water[yy, xx]:
                        mat[yy, xx] = TREES; v[yy, xx] = tv
                    # reflection: mirror about the horizon row
                    ry = int(round(2 * yh - yy))
                    if 0 <= ry < H and water[ry, xx] and rmat[ry, xx] in (SKY, MTN, TREES):
                        rmat[ry, xx] = TREES; rv[ry, xx] = tv - 0.02
            x += spacing * rng.uniform(0.35, 1.0) * (1.3 - 0.6 * clump)
EOF
python3 - <<'EOF'
import re
p='~/work/pixelart-studies/pedagogy/atelier/scene.py'
s=open(p).read()
# the backdrop no longer draws trees (hand-built sprites do); keep the far bank band
s=s.replace("""        tm, row, side = self.trees(az, el)
        tv = np.where(row > 0.5, 0.12, 0.19)
        tv = tv + 0.06 * ((side * np.sign(self.sun_side)) < 0) * (abs(self.sun_side) > 0.2)
        v = np.where(tm, tv, v); mat = np.where(tm, TREES, mat)
""","")
open(p,'w').write(s)
p='~/work/pixelart-studies/pedagogy/atelier/stage4_material.py'
s=open(p).read()
s=s.replace("    buf = render_buffers(P)\n","    import drawing\n    buf = drawing.finish(render_buffers(P))\n")
open(p,'w').write(s)
EOF
time nice -n 10 python3 stage4_material.py medium left close 0

# output:
medium_left_close_s0
nice -n 10 python3 stage4_material.py medium left close 0  12.01s user 3.90s system 92% cpu 17.183 total

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from core import *
from scene import *
import drawing
P=Params('medium','left','close',0)
buf=render_buffers(P)
mm=np.array([[210,210,230],[150,120,150],[40,60,40],[110,140,60],[230,120,90],[120,150,220]])
def show(m,name):
    im=mm[np.clip(m,0,5)]; im[m<0]=0
    Image.fromarray(im.astype(np.uint8)).resize((960,540),Image.NEAREST).save('stage4/'+name)
show(buf['mat'],'dbg_mat_before.png'); show(np.where(buf['water'],buf['rmat'],buf['mat']),'dbg_r_before.png')
drawing.finish(buf)
show(buf['mat'],'dbg_mat_after.png'); show(np.where(buf['water'],buf['rmat'],buf['mat']),'dbg_r_after.png')
print(drawing.horizon_row(buf['cam']))
"

# output:
38.04220376899461

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='lines.py'; s=open(p).read()
s=s.replace("""    pad = 3
    m = np.pad(mask.astype(float), pad, mode="edge")""","""    pad = 7                                                 # edge-replicate, then a ring of 0
    m = np.pad(np.pad(mask.astype(float), pad - 1, mode="edge"), 1, mode="constant")""")
s=s.replace('''def reassign(buf, keys, new_mask, old_mask):
    """pixels that joined the region copy their data from the nearest old-region pixel;
    pixels that left copy from the nearest pixel outside the old region."""''','''def reassign(buf, keys, new_mask, old_mask, within=None):
    """pixels that joined the region copy their data from the nearest old-region pixel;
    pixels that left copy from the nearest pixel outside the old region (inside `within`)."""''')
s=s.replace("""        _, (iy, ix) = ndimage.distance_transform_edt(old_mask, return_indices=True)""","""        src_out = ~old_mask if within is None else (within & ~old_mask)
        _, (iy, ix) = ndimage.distance_transform_edt(~src_out, return_indices=True)""")
open(p,'w').write(s)
p='drawing.py'; s=open(p).read()
s=s.replace("    return reassign(buf, keys, new, mask)","    return reassign(buf, keys, new, mask, within)")
open(p,'w').write(s)
EOF
nice -n 10 python3 stage4_material.py medium left close 0

# output:
medium_left_close_s0

cd ~/work/pixelart-studies/pedagogy/atelier && python3 -c "
from PIL import Image
im=Image.open('stage4/mat_medium_left_close_s0.png')
im.crop((0,0,480,240)).resize((960,480),Image.NEAREST).save('stage4/zoom_trees.png')
im.crop((0,240,480,540)).resize((960,600),Image.NEAREST).save('stage4/zoom_boulder.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='lines.py'; s=open(p).read()
s=s.replace("""    half_max = max(1.5, h * rng.uniform(0.17, 0.24))""","""    half_max = max(1.5, h * rng.uniform(0.22, 0.30))""")
s=s.replace("""        core = int(env * 0.45)""","""        core = int(round(env * rng.uniform(0.55, 0.75)))""")
s=s.replace("""        dy -= rng.choice([2, 2, 3])""","""        dy -= rng.choice([1, 2, 2])""")
s+='''
def blade(y, x, n, lean, ratio):
    """one grass blade from its foot (y, x) upward, n px long, leaning `lean` (-1, 0, 1)
    one pixel sideways every `ratio` px -- a clean 1:ratio line."""
    return [(y - k, x + lean * (k // ratio)) for k in range(n)]

def tuft(y, x, n, rng, spread=None):
    """a clump of 3-6 blades fanning out from one foot: the outer ones lean, the middle
    ones stand; lengths vary, the tallest near the middle."""
    k = rng.integers(3, 7) if spread is None else spread
    out = []
    for i in range(k):
        t = (i - (k - 1) / 2) / max(1, (k - 1) / 2)            # -1 .. 1
        lean = int(np.sign(t)) if abs(t) > 0.3 else 0
        ratio = int(rng.choice([2, 3, 3, 4])) if abs(t) < 0.9 else int(rng.choice([1, 2, 2]))
        L = max(1, int(round(n * (1 - 0.45 * abs(t)) * rng.uniform(0.6, 1.1))))
        out.append(blade(y, x + int(round(t * 1.2)), L, lean, ratio))
    return out
'''
open(p,'w').write(s)

p='paint.py'; s=open(p).read()
# calmer water
s=s.replace("""        r = hash2(seg, self.yy, self.seed + 3)
        lite = (r > 0.86) & w
        dark = (r < 0.10) & w""","""        r = hash2(seg, self.yy, self.seed + 3)
        calm = 0.5 + 0.5 * np.sin(self.yy * 0.21 + self.seed) * np.sin(self.yy * 0.057 + 1.3)
        lite = (r > 0.955 - 0.05 * calm) & w
        dark = (r < 0.025 + 0.04 * calm) & w""")
# granite flecks: gated by a patch noise, sparser
s=s.replace("""        lit = self.step > 2.6
        dk = rk & (np.where(lit, hsh < 0.09, hsh < 0.035)) & (detail > 0.15)
        lt = rk & (np.where(lit, hsh > 0.95, hsh > 0.985)) & (detail > 0.15)
        pair = rk & (hsh2 < 0.03) & lit & (detail > 0.4)""","""        lit = self.step > 2.6
        patch = np.zeros((self.H, self.W))
        patch[rk] = 0.5 + 0.5 * fbm3(q / (self.b_size * 0.12), self.seed + 12, 2)   # where grains gather
        dens = np.clip(patch * 1.6 - 0.3, 0, 1)
        dk = rk & (np.where(lit, hsh < 0.075 * dens, hsh < 0.02 * dens)) & (detail > 0.15)
        lt = rk & (np.where(lit, hsh > 1 - 0.04 * dens, hsh > 0.99)) & (detail > 0.15)
        pair = rk & (hsh2 < 0.025 * dens) & lit & (detail > 0.4)""")
s=s.replace("""        self.step[rk] += 0.45 * mott * detail[rk]""","""        self.step[rk] += 0.6 * mott * detail[rk]""")
# tufts rebuilt with lines.tuft
a=s.index("        # tufts: blades poke up over whatever is behind the top edge of a grass mass")
b=s.index("    def tree_marks(self):")
s=s[:a]+'''        # tufts: clumps of blades poke up over whatever is behind the top edge of a grass mass
        from lines import tuft
        self.tufts = []
        N = b["N"]; L = b["L"]
        top = gr & (shift(b["mat"], 1, 0, -1) != LAND)
        ys, xs = np.nonzero(top)
        code_g = RAMP_ORDER.index("grass")
        for y, x in zip(ys, xs):
            Lb = blade[y, x]
            if Lb < 1.5 or rng.random() > 0.22:
                continue
            base = self.step[y, x]
            for bl in tuft(y, x, int(round(Lb * rng.uniform(0.6, 1.2))), rng):
                for j, (yy, xx) in enumerate(bl):
                    if not (0 <= yy < H and 0 <= xx < W):
                        break
                    if depth[yy, xx] < depth[y, x] - 0.05 and b["mat"][yy, xx] == LAND:
                        continue
                    tipv = 0.9 if (j >= len(bl) - 2 and base > 3.3) else 0.0
                    self.tufts.append((yy, xx, code_g, base + tipv - 0.3 * (j == 0)))
        # grass round the boulder's foot overlaps the stone; lit where the stone is lit
        foot = (b["mat"] == BOULDER) & (shift(b["mat"], -1, 0, -1) == LAND)
        ys, xs = np.nonzero(foot)
        for y, x in zip(ys, xs):
            if rng.random() > 0.3:
                continue
            Lb = max(2, blade[min(y + 1, H - 1), x])
            ndl = float(np.dot(N[y, x], L))
            base = rng.uniform(4.3, 5.4) if ndl > 0.1 else rng.uniform(1.4, 2.4)
            for bl in tuft(y + 1, x, int(round(Lb * rng.uniform(0.7, 1.4))), rng):
                for j, (yy, xx) in enumerate(bl):
                    if not (0 <= yy < H and 0 <= xx < W):
                        break
                    tipv = 0.8 if j >= len(bl) - 2 and base > 3.3 else 0.0
                    self.tufts.append((yy, xx, code_g, base + tipv))

'''+s[b:]
# grass strokes: shorter, leaning sometimes
s=s.replace("""                off[y:y + n, x] = o
                y += n""","""                off[y:y + n, x] = o
                y += n + (1 if rng.random() < 0.25 else 0)""")
open(p,'w').write(s)
EOF
nice -n 10 python3 stage4_material.py medium left close 0 && python3 -c "
from PIL import Image
im=Image.open('stage4/mat_medium_left_close_s0.png')
im.crop((0,0,480,240)).resize((960,480),Image.NEAREST).save('stage4/zoom_trees.png')
im.crop((0,240,480,540)).resize((960,600),Image.NEAREST).save('stage4/zoom_boulder.png')"

# output:
medium_left_close_s0

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='paint.py'; s=open(p).read()
s=s.replace("""            ndl = float(np.dot(N[y, x], L))
            base = rng.uniform(4.3, 5.4) if ndl > 0.1 else rng.uniform(1.4, 2.4)""","""            base = rng.uniform(4.3, 5.4) if self.step[y, x] > 2.6 else rng.uniform(1.4, 2.4)""")
s=s.replace("blade = np.clip(0.28 / np.maximum(pix_m, 1e-4), 1, 9)       # a 28 cm blade, in px","blade = np.clip(0.2 / np.maximum(pix_m, 1e-4), 1, 7)        # a 20 cm stroke, in px")
open(p,'w').write(s)
EOF
J="medium left 10m 0|medium left 40m 0|large back 10m 1|small right close 2"
echo "$J" | tr '|' '\n' | while read a b c d; do (nice -n 10 python3 stage4_material.py $a $b $c $d >/dev/null 2>&1 &); done; sleep 2; until [ $(pgrep -f stage4_material | wc -l) -eq 0 ]; do sleep 3; done
python3 -c "
from PIL import Image
import sys; sys.path.insert(0,'.')
from core import grid, label
n=['medium_left_10m_s0','medium_left_40m_s0','large_back_10m_s1','small_right_close_s2']
grid([label(Image.open(f'stage4/mat_{x}.png'),x) for x in n],2).save('stage4/v1_sheet.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && grep -n "back = \|rim\b\|rim =" paint.py | head; nice -n 10 python3 -c "
from core import *
from scene import *
import drawing
from paint import Painter
from backdrop_util import sun_side_of
P=Params('large','back','10m',1)
buf=drawing.finish(render_buffers(P))
print('L',buf['L'],'f',buf['cam'].f, np.dot(buf['L'],buf['cam'].f))
N=buf['N']; rk=buf['mat']==4
ndl=(N*buf['L']).sum(-1)
print('ndl boulder pct', np.percentile(ndl[rk],[5,50,95,99]))
"

# output:
169:        # 2) silhouette: a dark accent where the form turns away from the light; a light rim
175:        back = L[2] * 0 + np.sum(L * self.cam_f) > 0.3             # sun behind the subject
177:            rim = sil & (ndl > -0.45)
178:            self.step[rim] = np.maximum(self.step[rim], 3.4)
L [0.22034887 0.15643447 0.96279522] f [ 0.22034887 -0.15643447  0.96279522] 0.9510565162951538
ndl boulder pct [-0.94455397 -0.71448584  0.03963216  0.20416874]

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='paint.py'; s=open(p).read()
s=s.replace("""            rim = sil & (ndl > -0.45)
            self.step[rim] = np.maximum(self.step[rim], 3.4)
            rim2 = rk & ~sil & (ndl > -0.1) & shift(sil, 0, 1, False)
            self.step[rim2] = np.maximum(self.step[rim2], 2.6)""","""            # contre-jour: the silhouette facing up/sideways into the light burns bright,
            # the pixel inside it half as much (a 2-px rim that thins as it turns away)
            up_edge = N[..., 1] > -0.2
            rim = sil & (ndl > -0.5) & up_edge
            self.step[rim] = np.maximum(self.step[rim], np.where(ndl[rim] > -0.2, 4.8, 3.6))
            inner = rk & ~sil & up_edge & (ndl > -0.35) & (shift(rim, 1, 0, False) | shift(rim, 0, 1, False) | shift(rim, 0, -1, False))
            self.step[inner] = np.maximum(self.step[inner], 3.2)""")
# back-lit grass: blade tips glow
s=s.replace("""        lit = self.step > 3.3
        o2 = np.where(lit, off, np.where(off > 0, off * 0.5, off * 0.8))""","""        lit = self.step > 3.3
        o2 = np.where(lit, off, np.where(off > 0, off * 0.5, off * 0.8))
        if np.sum(b["L"] * self.cam_f) > 0.3:
            # contre-jour grass: translucent blades light up; the light strokes carry it
            o2 = np.where(off > 0, 2.2, o2)""")
open(p,'w').write(s)
p='scene.py'; s=open(p).read()
s=s.replace("""        v = 0.88 - 0.26 * e ** 0.8 + 0.07 * np.clip(toward_sun, 0, 1) ** 3 * (1 - e)""","""        v = 0.88 - 0.26 * e ** 0.8 + 0.07 * np.clip(toward_sun, 0, 1) ** 3 * (1 - e)
        # the glow round the low sun (in view only when the light is behind the subject)
        ang = np.arccos(np.clip(dot(d, self.L), -1, 1))
        v = v + 0.10 * np.exp(-ang / 0.12) + 0.05 * np.exp(-ang / 0.03)""")
open(p,'w').write(s)
EOF
nice -n 10 python3 stage4_material.py large back 10m 1 >/dev/null && echo ok

# output:
ok

cd ~/work/pixelart-studies/pedagogy/atelier && cat > palette.py <<'EOF'
"""The palette: 32 colours chosen by hand, each placed to hit a VALUE the grey studies
already proved. Ramps are lists of palette indices (dark -> light); ramps share colours
where a value and a hue coincide (the water borrows the rock's and mountains' darks and
the sky's lilacs -- it is, after all, mostly reflecting them).

Dusk logic: the light is warm (peach, orange, ochre, cream), the shadows are cool and
violet (the sky fills them); greens are pushed to olive and ochre because the low sun is
orange and there is a lot of pink sky in the bounce.
"""
import numpy as np

COLORS = [
    # sky: lilac zenith -> peach horizon -> cream glow                         (0-5)
    (170, 140, 180), (190, 155, 190), (212, 168, 190), (236, 190, 186), (248, 210, 186), (255, 228, 196),
    # mountain rock: violet shadow -> rose -> orange-lit -> peach              (6-11)
    (98, 78, 110), (122, 96, 126), (170, 114, 120), (222, 136, 118), (238, 170, 150), (250, 208, 190),
    # spruce: near-black -> brown-olive -> lit needles                         (12-15)
    (16, 14, 8), (34, 26, 18), (54, 44, 24), (88, 64, 34),
    # grass: black-olive -> olive -> ochre                                     (16-22)
    (14, 18, 6), (34, 36, 14), (56, 54, 24), (88, 76, 34), (132, 100, 48), (178, 126, 64), (214, 158, 86),
    # granite: violet-black -> violet-grey -> mauve halftone -> warm pink-grey  (23-28)
    (22, 18, 28), (44, 36, 52), (70, 58, 76), (140, 112, 116), (196, 152, 140), (246, 204, 178),
    # water's own: mid violet; deep peaty violet; a glint                      (29-31)
    (140, 116, 148), (38, 32, 52), (255, 242, 214),
]
RAMPS = {
    "sky":   [0, 1, 2, 3, 4, 5],
    "mtn":   [6, 7, 8, 9, 10, 11],
    "trees": [12, 13, 14, 15],
    "grass": [16, 17, 18, 19, 20, 21, 22],
    "rock":  [23, 24, 25, 26, 27, 28],
    "water": [30, 25, 6, 7, 29, 0, 1, 31],
}

def lum(c):
    c = np.asarray(c, float)
    return (0.299 * c[..., 0] + 0.587 * c[..., 1] + 0.114 * c[..., 2]) / 255.0

def ramp_lums(colors=COLORS, ramps=RAMPS):
    return {k: [float(lum(colors[i])) for i in v] for k, v in ramps.items()}

def grey_colors(colors=COLORS):
    out = []
    for c in colors:
        g = round(float(lum(c)) * 255); out.append((g, g, g))
    return out

if __name__ == "__main__":
    assert len(COLORS) <= 32, len(COLORS)
    for k, v in ramp_lums().items():
        print(k, [round(x, 3) for x in v])
    print(len(COLORS), "colours")
EOF
python3 palette.py

# output:
sky [0.602, 0.665, 0.72, 0.797, 0.857, 0.911]
mtn [0.344, 0.42, 0.515, 0.626, 0.737, 0.857]
trees [0.055, 0.108, 0.175, 0.266]
grass [0.061, 0.129, 0.201, 0.293, 0.406, 0.527, 0.653]
rock [0.08, 0.158, 0.25, 0.474, 0.642, 0.838]
water [0.141, 0.25, 0.344, 0.42, 0.497, 0.602, 0.665, 0.952]
32 colours

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='paint.py'; s=open(p).read()
s=s.replace('from scene import SKY, MTN, TREES, LAND, BOULDER, WATER','from scene import SKY, MTN, TREES, LAND, BOULDER, WATER, SNOW\nimport palette as PAL')
a=s.index("# ------------------------------------------------------------------ ramps (values)")
b=s.index("# ------------------------------------------------------------------ helpers")
s=s[:a]+'''# ------------------------------------------------------------------ ramps
# each ramp is a list of palette indices, dark -> light (palette.py). Their luminances are
# the values the grey stages were built on.
RAMP_ORDER = ["sky", "mtn", "trees", "grass", "rock", "water"]
RAMP_LUM = PAL.ramp_lums()

def grey_palette():
    return np.array(PAL.grey_colors())

def colour_palette():
    return np.array(PAL.COLORS)

'''+s[b:]
s=s.replace("""        self.off, self.n = palette_layout(ramps)
""","""        self.ridx = {k: np.array(v) for k, v in PAL.RAMPS.items()}
""")
s=s.replace("""        m = mat == SKY;     code[m] = R["sky"];   step[m] = (v[m] - 0.60) / (0.95 - 0.60) * 5""","""        m = mat == SKY;     code[m] = R["sky"];   step[m] = (v[m] - 0.60) / (0.95 - 0.60) * 5
        m = mat == SNOW;    code[m] = R["sky"];   step[m] = (v[m] - 0.60) / (0.95 - 0.60) * 5""")
s=s.replace("sky_like = np.isin(rmat, [SKY, MTN])","sky_like = np.isin(rmat, [SKY, MTN, SNOW])")
# quantize -> levels
a=s.index("    def quantize(self):")
b=s.index("    def paint(self, boulder, cam_f, sun_side, b_size):")
s=s[:a]+'''    def quantize(self):
        lev = np.zeros((self.H, self.W), int)
        for i, name in enumerate(RAMP_ORDER):
            m = self.code == i
            if not m.any():
                continue
            n = len(self.ramps[name])
            v = np.clip(self.step / (n - 1), 0, 1)
            q = quantize(np.where(m, v, 0), n, kind="check", strength=0.55, flat_snap=0.1)
            lev[m] = q[m]
        return lev

    def to_index(self, code, lev):
        idx = np.zeros_like(lev)
        for i, name in enumerate(RAMP_ORDER):
            m = code == i
            idx[m] = self.ridx[name][np.clip(lev[m], 0, len(self.ridx[name]) - 1)]
        return idx

'''+s[b:]
a=s.index("        idx = self.quantize()")
b=s.index("        return idx\n", a)
s=s[:a]+'''        lev = self.quantize()
        code = self.code.copy()
        # clean strays on the combined (code, level) map
        key = remove_orphans(code * 16 + lev)
        code, lev = key // 16, key % 16
        # marks that are MEANT to be single pixels go on after the clean-up
        if self.fleck_dark is not None:
            rk = code == RAMP_ORDER.index("rock")
            n = len(self.ramps["rock"])
            lev[self.fleck_dark & rk] = np.maximum(lev[self.fleck_dark & rk] - 1, 0)
            lev[self.fleck_lite & rk] = np.minimum(lev[self.fleck_lite & rk] + 1, n - 1)
        T = threshold_map(self.H, self.W, "check")
        for (y, x, c, s_) in self.tufts:
            name = RAMP_ORDER[c]; n = len(self.ramps[name])
            st = np.clip(s_, 0, n - 1)
            lo = int(np.floor(st)); fr = st - lo
            code[y, x] = c; lev[y, x] = min(lo + (fr > T[y, x]), n - 1)
        self.final_code, self.final_lev = code, lev
        idx = self.to_index(code, lev)
'''+s[b:]
open(p,'w').write(s)

p='scene.py'; s=open(p).read()
s=s.replace("SKY, MTN, TREES, LAND, BOULDER, WATER = 0, 1, 2, 3, 4, 5","SKY, MTN, TREES, LAND, BOULDER, WATER, SNOW = 0, 1, 2, 3, 4, 5, 6")
s=s.replace("""        vm = np.where(lit, 0.64, 0.44) + 0.08 * snow
        vm = vm * 0.75 + v * 0.25                               # aerial perspective: toward sky
        v = np.where(m, vm, v); mat = np.where(m, MTN, mat)""","""        vm = np.where(lit, 0.64, 0.44)
        vm = vm * 0.8 + v * 0.2                                 # aerial perspective: toward sky
        vs = np.where(lit, 0.90, 0.64)                          # snow: lilac in shade, cream lit
        v = np.where(m, np.where(snow, vs, vm), v)
        mat = np.where(m, np.where(snow, SNOW, MTN), mat)""")
open(p,'w').write(s)
p='drawing.py'; s=open(p).read()
s=s.replace("from scene import SKY, MTN, TREES, LAND, BOULDER, WATER, _rnd","from scene import SKY, MTN, TREES, LAND, BOULDER, WATER, SNOW, _rnd")
s=s.replace("""        clean_region(buf, (buf["mat"] == MTN), KEYS, sigma=0.8, tol=0.6)
        clean_region(buf, w & (buf["rmat"] == MTN), RKEYS, within=w, sigma=0.8, tol=0.6)""","""        clean_region(buf, np.isin(buf["mat"], [MTN, SNOW]), KEYS, sigma=0.8, tol=0.6)
        clean_region(buf, (buf["mat"] == SNOW), KEYS, sigma=0.8, tol=0.6)
        clean_region(buf, w & np.isin(buf["rmat"], [MTN, SNOW]), RKEYS, within=w, sigma=0.8, tol=0.6)""")
s=s.replace("mat[yy, xx] in (SKY, MTN, TREES)","mat[yy, xx] in (SKY, MTN, SNOW, TREES)")
s=s.replace("rmat[ry, xx] in (SKY, MTN, TREES)","rmat[ry, xx] in (SKY, MTN, SNOW, TREES)")
open(p,'w').write(s)
p='stage4_material.py'; s=open(p).read()
s=s.replace("from paint import Painter, grey_palette","from paint import Painter, grey_palette, colour_palette")
s=s.replace('''    to_image(idx, grey_palette(), 3).save(f"{OUT}/mat_{name}.png")''','''    to_image(idx, grey_palette(), 3).save(f"{OUT}/mat_{name}.png")
    if len(a) > 4 and a[4] == "colour":
        os.makedirs("stage5", exist_ok=True)
        to_image(idx, colour_palette(), 3).save(f"stage5/col_{name}.png")''')
open(p,'w').write(s)
EOF
grep -n "SKY, MTN" paint.py | head -3; nice -n 10 python3 stage4_material.py medium left close 0 colour

# output:
10:from scene import SKY, MTN, TREES, LAND, BOULDER, WATER, SNOW
91:        sky_like = np.isin(rmat, [SKY, MTN, SNOW])
260:            edge = t & np.isin(nb, [SKY, MTN])
medium_left_close_s0

cd ~/work/pixelart-studies/pedagogy/atelier && cat > method.py <<'EOF'
"""THE METHOD: paint 'a granite boulder in a grass bank at the edge of a still mountain tarn,
at dusk' from parameters, as 320x180 indexed pixel art (32 colours), shown 3x nearest.

  python3 method.py one --size medium --light left --dist 10m --seed 0 [--grey]
  python3 method.py sheets            # every parameter sheet (4 worker processes, 1 thread each)

Pipeline (each stage is one of the atelier lessons):

  1. geometry  (geometry.py)  boulder = convex hull of a lopsided mass, weathered edges,
                              spalls, a joint; bank = cove shoreline + hummocks + turf lip
  2. value     (valuemodel.py, scene.py)  six-step value plan per material: light family /
                              shadow family kept apart, core, reflected light, occlusion;
                              grass keyed to the sun; water = reflection x Fresnel
  3. drawing   (drawing.py, lines.py)  contours re-drawn in clean pixel runs; spruce built
                              from pixel strokes; reflections mirrored about the horizon
  4. material  (paint.py)     granite ridges/creases/flecks/mottle; grass strokes & tufts;
                              ripple dashes; lit tree edges; contre-jour rims
  5. colour    (palette.py)   32 hand-placed colours on the proven values
"""
from core import *
import os, sys, argparse
from multiprocessing import Pool
from scene import Params, render_buffers, LIGHTS, SIZES, DISTS
import drawing
from paint import Painter, grey_palette, colour_palette
from backdrop_util import sun_side_of
import palette as PAL

OUT = "sheets"

def paint(size="medium", light="left", dist="10m", seed=0, grey=False):
    P = Params(size, light, dist, seed)
    buf = drawing.finish(render_buffers(P))
    pt = Painter(buf, seed=seed)
    idx = pt.paint(buf["st"].boulder, buf["cam"].f, sun_side_of(buf), P.bsize)
    return idx

def cell_path(size, light, dist, seed, grey=False):
    os.makedirs(f"{OUT}/cells", exist_ok=True)
    return f"{OUT}/cells/{size}_{light}_{dist}_s{seed}{'_grey' if grey else ''}.png"

def render_cell(args):
    size, light, dist, seed = args
    path = cell_path(size, light, dist, seed)
    if not os.path.exists(path):
        idx = paint(size, light, dist, seed)
        np.save(path.replace(".png", ".npy"), idx.astype(np.uint8))
        to_image(idx, colour_palette(), 3).save(path)
    return path

def sheet(cells, cols, title, name, row_labels=None, col_labels=None):
    ims = []
    for c in cells:
        im = Image.open(cell_path(*c))
        ims.append(label(im, f"{c[0]} / {c[1]} / {c[2]} / seed {c[3]}"))
    g = grid(ims, cols, pad=8, bg=(24, 20, 28))
    from PIL import ImageDraw
    out = Image.new("RGB", (g.width, g.height + 30), (24, 20, 28))
    out.paste(g, (0, 30))
    ImageDraw.Draw(out).text((10, 8), title, fill=(235, 220, 230))
    out.save(f"{OUT}/{name}.png")
    # a half-size preview for quick viewing
    out.resize((out.width // 3, out.height // 3), Image.NEAREST).save(f"{OUT}/{name}_preview.png")

SHEETS = {
    "sheet_light":   ("light direction (medium boulder, 10 m, seeds 0 and 1)",
                      [("medium", l, "10m", s) for s in (0, 1) for l in LIGHTS], 4),
    "sheet_size_distance": ("boulder size (rows) x distance (columns), low sun from the left, seed 2",
                      [(z, "left", d, 2) for z in SIZES for d in DISTS], 3),
    "sheet_light_distance": ("light direction (rows) x distance (columns), medium boulder, seed 3",
                      [("medium", l, d, 3) for l in LIGHTS for d in DISTS], 3),
    "sheet_seed":    ("random seed (medium boulder, 10 m, light from the left)",
                      [("medium", "left", "10m", s) for s in range(8)], 4),
    "sheet_seed_close": ("random seed, close, back-lit",
                      [("large", "back", "close", s) for s in range(4)], 2),
}

def palette_card():
    os.makedirs(OUT, exist_ok=True)
    from PIL import ImageDraw
    rows = list(PAL.RAMPS.items())
    sw = 40
    im = Image.new("RGB", (120 + 8 * sw, len(rows) * (sw * 2 + 10) + 10), (24, 20, 28))
    d = ImageDraw.Draw(im)
    for r, (name, ids) in enumerate(rows):
        y = 10 + r * (sw * 2 + 10)
        d.text((8, y + sw // 2), name, fill=(230, 220, 230))
        for k, i in enumerate(ids):
            c = tuple(PAL.COLORS[i]); g = round(float(PAL.lum(c)) * 255)
            d.rectangle([110 + k * sw, y, 110 + (k + 1) * sw - 2, y + sw - 2], fill=c)
            d.rectangle([110 + k * sw, y + sw, 110 + (k + 1) * sw - 2, y + 2 * sw - 2], fill=(g, g, g))
    im.save(f"{OUT}/palette_card.png")

if __name__ == "__main__":
    ap = argparse.ArgumentParser()
    ap.add_argument("cmd", choices=["one", "sheets"])
    ap.add_argument("--size", default="medium"); ap.add_argument("--light", default="left")
    ap.add_argument("--dist", default="10m"); ap.add_argument("--seed", type=int, default=0)
    ap.add_argument("--grey", action="store_true"); ap.add_argument("--only", default=None)
    a = ap.parse_args()
    if a.cmd == "one":
        idx = paint(a.size, a.light, a.dist, a.seed)
        pal = grey_palette() if a.grey else colour_palette()
        p = cell_path(a.size, a.light, a.dist, a.seed, a.grey)
        to_image(idx, pal, 3).save(p); print(p)
    else:
        palette_card()
        names = [a.only] if a.only else list(SHEETS)
        jobs = sorted({c for n in names for c in SHEETS[n][1]})
        with Pool(4) as pool:                       # studio rule: <= 4 threads in total
            for p in pool.imap_unordered(render_cell, jobs):
                print(p, flush=True)
        for n in names:
            title, cells, cols = SHEETS[n]
            sheet(cells, cols, title, n)
        print("sheets done")
EOF
rm -rf sheets; time nice -n 10 python3 method.py sheets 2>&1 | tail -3

# output:
sheets/cells/small_left_40m_s2.png
sheets/cells/small_left_close_s2.png
sheets done
nice -n 10 python3 method.py sheets 2>&1  519.91s user 271.89s system 361% cpu 3:38.93 total
tail -3  0.00s user 0.00s system 0% cpu 3:38.93 total

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='scene.py'; s=open(p).read()
a=s.index("    def mountains(self, az):")
b=s.index("    def trees(self, az, el):")
s=s[:a]+'''    def build_massifs(self):
        """a few big peaks, each a pyramid of PLANES: flanks are straight lines, and the faces
        between them fan out from the summit along clean pixel ratios (1:1, 2:1, 1:2, 3:1)."""
        rng = np.random.default_rng(self.seed + 4242)
        self.peaks = []
        az = -0.75
        while az < 0.75:
            h = np.radians(rng.uniform(2.6, 6.2))
            sl = rng.choice([0.5, 0.5, 1.0, 0.33])            # flank slopes (el per az): 2:1, 1:1, 3:1
            sr = rng.choice([0.5, 0.5, 1.0, 0.33])
            # spur directions (as slopes from the summit), sorted outward, each side
            nice = np.array([0.25, 0.33, 0.5, 1.0, 2.0])
            spl = np.sort(rng.choice(nice, rng.integers(1, 3), replace=False))
            spr = np.sort(rng.choice(nice, rng.integers(1, 3), replace=False))
            snow = rng.uniform(0.35, 0.7)
            self.peaks.append((az, h, sl, sr, spl, spr, snow, rng.integers(0, 2)))
            az += rng.uniform(0.12, 0.3)

    def mountains(self, az, detail=True):
        if not hasattr(self, "peaks"):
            self.build_massifs()
        em = np.full(az.shape, -1.0); own = np.zeros(az.shape, int)
        for i, (c, h, sl, sr, *_ ) in enumerate(self.peaks):
            e = np.where(az < c, h - (c - az) * sl, h - (az - c) * sr)
            own = np.where(e > em, i, own); em = np.maximum(em, e)
        self._own = own
        return em

    def massif_faces(self, az, el):
        """-> (face: +1 turned left / -1 turned right / 0 facing us, snow bool)"""
        em = self.mountains(az); own = self._own
        face = np.zeros(az.shape); snow = np.zeros(az.shape, bool)
        for i, (c, h, sl, sr, spl, spr, sn, flip) in enumerate(self.peaks):
            m = own == i
            if not m.any():
                continue
            dx = az[m] - c; dy = h - el[m]                        # below-summit offset
            slope = dy / np.maximum(np.abs(dx), 1e-6)            # steepness of the ray from the summit
            left = dx < 0
            f = np.where(left, 1.0, -1.0)                        # outer flank faces its own side
            # inner sectors (steeper than a spur) alternate, like buttresses and gullies
            for k, s_ in enumerate(spl):
                f = np.where(left & (slope > s_ * 1.0 / max(sl, 1e-3) * sl * 1.6 + 0.0) , (-1.0) ** (k + 1 + flip), f)
            for k, s_ in enumerate(spr):
                f = np.where(~left & (slope > s_ * 1.6), (-1.0) ** (k + flip), f)
            face[m] = f
            # snow: the upper part, deeper in the gullies, in straight-edged tongues
            depth = sn * h * (1.0 + 0.6 * (f < 0))
            tongue = (np.floor(np.abs(dx) / np.radians(0.5)) % 3) == 0
            snow[m] = (dy < depth * np.where(tongue, 1.6, 1.0)) & (dy < h * 0.9)
        return em, face, snow

'''+s[b:]
a=s.index("        em = self.mountains(az)\n        m = el < em")
b=s.index("        vm = np.where(lit, 0.64, 0.44)")
s=s[:a]+'''        em, face, snow = self.massif_faces(az, el)
        m = el < em
        if abs(self.sun_side) > 0.2:
            lit = face * np.sign(self.sun_side) > 0
        else:
            lit = face == face                                    # front light: all faces lit
            if dot(self.L, self.fh) > 0.3:
                lit = np.zeros_like(lit)                          # back light: all in shade
'''+s[b:]
open(p,'w').write(s)
EOF
nice -n 10 python3 method.py one --size medium --light left --dist 10m --seed 0 && python3 -c "
from PIL import Image; im=Image.open('sheets/cells/medium_left_10m_s0.png'); im.crop((0,0,960,240)).save('stage5/crop_mtn.png')"

# output:
sheets/cells/medium_left_10m_s0.png

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='scene.py'; s=open(p).read()
s=s.replace("""                f = np.where(left & (slope > s_ * 1.0 / max(sl, 1e-3) * sl * 1.6 + 0.0) , (-1.0) ** (k + 1 + flip), f)""","""                f = np.where(left & (slope > s_ * 1.6), (-1.0) ** (k + 1 + flip), f)""")
s=s.replace("""            depth = sn * h * (1.0 + 0.6 * (f < 0))
            tongue = (np.floor(np.abs(dx) / np.radians(0.5)) % 3) == 0
            snow[m] = (dy < depth * np.where(tongue, 1.6, 1.0)) & (dy < h * 0.9)""","""            # snow line: a zigzag of straight couloirs (a triangle wave = clean 1:1 / 2:1 runs)
            period = np.radians(0.9)
            tri = np.abs(((np.abs(dx) / period) % 1.0) - 0.5) * 2          # 0..1..0
            depth = sn * h * (0.55 + 0.45 * (f < 0)) + tri * period * (0.5 if flip else 1.0)
            snow[m] = dy < depth""")
s=s.replace("snow = rng.uniform(0.35, 0.7)","snow = rng.uniform(0.18, 0.4)")
open(p,'w').write(s)
EOF
nice -n 10 python3 method.py one --size medium --light left --dist 10m --seed 0 >/dev/null && python3 -c "
from PIL import Image; im=Image.open('sheets/cells/medium_left_10m_s0.png'); im.crop((0,0,960,240)).save('stage5/crop_mtn.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='scene.py'; s=open(p).read()
s=s.replace("""        vm = np.where(lit, 0.64, 0.44)
        vm = vm * 0.8 + v * 0.2 """,""" 
        vm = np.where(lit, 0.67, 0.39)
        vm = vm * 0.88 + v * 0.12 """)
open(p,'w').write(s)
p='paint.py'; s=open(p).read()
s=s.replace("""        self.step[lite] += 0.75
        self.step[dark] -= 0.75""","""        # ripples only where the surface is disturbed: fewer toward the far shore
        far_calm = np.clip(1.2 - near * 2.5, 0, 1)
        keep = hash2(self.xx // 7, self.yy, self.seed + 4) > far_calm * 0.8
        self.step[lite & keep] += 0.75
        self.step[dark & keep] -= 0.75
        # a catspaw: a band of breeze-roughened water that shows a paler, flatter sky
        rng = np.random.default_rng(self.seed + 31)
        rows = np.nonzero(w.any(1))[0]
        if len(rows) > 20:
            for _ in range(rng.integers(1, 3)):
                cy = rng.uniform(rows.min() + 0.25 * len(rows), rows.max() - 0.1 * len(rows))
                hw = rng.uniform(2, 5) * (0.5 + near[int(cy)].mean())
                x0 = rng.uniform(-0.3, 0.6) * W; x1 = x0 + rng.uniform(0.4, 0.9) * W
                edge = np.clip(1 - np.abs(self.yy - cy) / hw, 0, 1) * np.clip(np.minimum(self.xx - x0, x1 - self.xx) / 30, 0, 1)
                band = w & (edge > 0)
                self.step[band] += 0.9 * edge[band]""")
open(p,'w').write(s)
EOF
nice -n 10 python3 method.py one --size medium --light left --dist 10m --seed 0 >/dev/null; echo

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='paint.py'; s=open(p).read()
s=s.replace("""    def rock_marks(self):""","""    def cap_water(self):
        # the glint colour is reserved for the sun's own reflection
        b = self.b; w = b["water"] & (self.code == RAMP_ORDER.index("water"))
        n = len(self.ramps["water"])
        sun = b["rv"] > 0.95
        self.step[w & ~sun] = np.minimum(self.step[w & ~sun], n - 2 + 0.2)

    def rock_marks(self):""")
s=s.replace("self.water_marks(); self.rock_marks();","self.water_marks(); self.cap_water(); self.rock_marks();")
open(p,'w').write(s)
EOF
rm -rf sheets/cells; nice -n 10 python3 method.py sheets 2>&1 | tail -1

# output:
sheets done

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from core import *
from scene import *
import drawing
P=Params('medium','back','close',3)
buf=render_buffers(P)
m0=buf['mat'].copy()
drawing.finish(buf)
m1=buf['mat']
for y in range(12,22):
    print(y, ''.join(str(v) for v in m1[y,100:200]))
"

# output:
12 2220002222200022222226622226622226666666662662222200000000000000200002000000000000000000000000000000
13 2220002222206622222226662226622222266666622262222266600000000000200002000000000000000000000000000000
14 2220002222266222222226622222222222266666622262222266666600000002220022200000000000000000000000000000
15 2222222222226622222226622222222222226666622262222226666666600000200002000000000000000000000000000000
16 2222222222226622222222222222222222226666222222222226666666666602220022200000000000000000000000000000
17 2222222222222222222226622222222222226666622222222226666666266602220122211111111111111116660066666666
18 2222222222226222222222222222222222221666222222222266666661266612222122220000000000000006666666666666
19 2222222222222222222222222222222222222266222222222222112612226112221122200000000000000166661666666666
20 2222222222222222222222222222222222222226222222222221112212226122221222211100000200011116616666666666
21 2222222222222222222222222222222222222212222222222222222212221112221122221111102221111111111161116666

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='lines.py'; s=open(p).read()
a=s.index("    polys.sort(key=lambda t: -t[0])")
b=s.index("    res = np.array(out, bool)[pad:-pad, pad:-pad]")
s=s[:a]+'''    # even-odd composition: outer outlines and holes are nested closed contours, so XOR
    acc = np.zeros((H + 2 * pad, W + 2 * pad), bool)
    for _, poly in polys:
        img = Image.new("L", (W + 2 * pad, H + 2 * pad), 0)
        ImageDraw.Draw(img).polygon([tuple(p) for p in poly], fill=1)
        acc ^= np.array(img, bool)
    out = Image.fromarray(acc.astype(np.uint8))
'''+s[b:]
open(p,'w').write(s)
EOF
nice -n 10 python3 -c "
from core import *
from scene import *
import drawing
P=Params('medium','back','close',3)
buf=drawing.finish(render_buffers(P))
m1=buf['mat']
for y in range(12,22):
    print(y, ''.join(str(v) for v in m1[y,100:200]))
"

# output:
12 2220002222200022222226622226622226666666662662222200000000000000200002000000000000000000000000000000
13 2220002222206622222226662226622222266666622262222266600000000000200002000000000000000000000000000000
14 2220002222266222222226622222222222266666622262222266666600000002220022200000000000000000000000000000
15 2222222222226622222226622222222222226666622262222226666666600000200002000000000000000000000000000000
16 2222222222226622222222222222222222226666222222222226666666666602220022200000000000000000000000000000
17 2222222222222222222226622222222222226666622222222226666666266602220122211111111111111116660066666666
18 2222222222226222222222222222222222221666222222222266666661266612222122220000000000000006666666666666
19 2222222222222222222222222222222222222266222222222222112612226112221122200000000000000166661666666666
20 2222222222222222222222222222222222222226222222222221112212226122221222211100000200011116616666666666
21 2222222222222222222222222222222222222212222222222222222212221112221122221111102221111111111161116666

cd ~/work/pixelart-studies/pedagogy/atelier && nice -n 10 python3 -c "
from core import *
from scene import *
P=Params('medium','back','close',3)
buf=render_buffers(P)
m1=buf['mat']
for y in range(14,22):
    print(y, ''.join(str(v) for v in m1[y,160:230]))
"

# output:
14 0000000000000000000000000000000000000000000000000000000000000000000000
15 0000000000000000000000000000000000000000000000000000000000000000000000
16 0000000000000000000000000000000666600000000000000000000000000000000000
17 6610000000000000000000000000066666666000000000000000000000000000000000
18 6611660000000000000000000006661666666660000000000000000000000000000000
19 6111161110000000000000000166611666616666600000000000000000000000000000
20 6111161111610000000000011116116666116666666000000000000000000000000000
21 1111111111111110000001111111161661111666116660000000000000000000000000

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='lines.py'; s=open(p).read()
s=s.replace("""    res = np.array(out, bool)[pad:-pad, pad:-pad]
    return res""","""    res = np.array(out, bool)[pad:-pad, pad:-pad]
    # snapping can fold an edge back on itself and leave a 1-px sliver: drop NEW pixels
    # that are one pixel thick (nothing above and below, or nothing left and right)
    for _ in range(2):
        p_ = np.pad(res, 1)
        up, dn, lf, rt = p_[:-2, 1:-1], p_[2:, 1:-1], p_[1:-1, :-2], p_[1:-1, 2:]
        thin = res & ~mask & ((~up & ~dn) | (~lf & ~rt))
        res = res & ~thin
    return res""")
open(p,'w').write(s)
p='paint.py'; s=open(p).read()
s=s.replace("""    def cap_water(self):""","""    def glitter(self):
        \"\"\"contre-jour only: under the sun the ripples that face it flash -- a column of
        glints widening toward the viewer.\"\"\"
        b = self.b; w = b["water"]
        if not (np.sum(b["L"] * self.cam_f) > 0.3) or not w.any():
            return
        cam = self.cam
        x_s, y_s, z = cam.project(cam.pos + b["L"] * 1e4)
        near = np.clip(np.abs(b["dwy"]) / 0.35, 0, 1)
        half = 6 + 40 * near
        col = np.abs(self.xx - x_s) < half * (0.6 + 0.4 * hash2(self.yy, 3, self.seed))
        h = hash2(self.xx // 3, self.yy, self.seed + 21)
        g = w & col & (h > 0.93 - 0.05 * (1 - np.abs(self.xx - x_s) / half))
        self.glint = g

    def cap_water(self):""")
s=s.replace("self.water_marks(); self.cap_water();","self.water_marks(); self.cap_water(); self.glitter();")
s=s.replace("""        self.fleck_dark = self.fleck_lite = None
        self.tufts = []""","""        self.fleck_dark = self.fleck_lite = None
        self.tufts = []; self.glint = None""")
s=s.replace("""        T = threshold_map(self.H, self.W, "check")
        for (y, x, c, s_) in self.tufts:""","""        if self.glint is not None:
            wc = RAMP_ORDER.index("water")
            g = self.glint & (code == wc)
            lev[g] = len(self.ramps["water"]) - 1
        T = threshold_map(self.H, self.W, "check")
        for (y, x, c, s_) in self.tufts:""")
s=s.replace("    def paint(self, boulder, cam_f, sun_side, b_size):\n        self.boulder, self.cam_f, self.sun_side, self.b_size = boulder, cam_f, sun_side, b_size",
            "    def paint(self, boulder, cam_f, sun_side, b_size, cam=None):\n        self.boulder, self.cam_f, self.sun_side, self.b_size = boulder, cam_f, sun_side, b_size\n        self.cam = cam")
open(p,'w').write(s)
p='method.py'; s=open(p).read()
s=s.replace("idx = pt.paint(buf[\"st\"].boulder, buf[\"cam\"].f, sun_side_of(buf), P.bsize)","idx = pt.paint(buf[\"st\"].boulder, buf[\"cam\"].f, sun_side_of(buf), P.bsize, cam=buf[\"cam\"])")
open(p,'w').write(s)
EOF
grep -n "cam=buf" method.py; nice -n 10 python3 method.py one --size medium --light back --dist close --seed 3 >/dev/null; echo

# output:
34:    idx = pt.paint(buf["st"].boulder, buf["cam"].f, sun_side_of(buf), P.bsize, cam=buf["cam"])

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='paint.py'; s=open(p).read()
s=s.replace("""        x_s, y_s, z = cam.project(cam.pos + b["L"] * 1e4)""","""        x_s, y_s, z = cam.project(cam.pos + b["L"] * 1e4)
        _, y_r, _ = cam.project(cam.pos + b["L"] * np.array([1, -1, 1]) * 1e4)   # the sun's image""")
s=s.replace("""        g = w & col & (h > 0.93 - 0.05 * (1 - np.abs(self.xx - x_s) / half))""","""        g = w & col & (h > 0.93 - 0.05 * (1 - np.abs(self.xx - x_s) / half)) & (self.yy > y_r - 14)""")
open(p,'w').write(s)
EOF
rm -rf sheets/cells; nice -n 10 python3 method.py sheets 2>&1 | tail -1

# output:
sheets done

cd ~/work/pixelart-studies/pedagogy/atelier && cat > lessons_figs.py <<'EOF'
"""figures for LESSONS.md, all from the final pipeline"""
from core import *
import os
from scene import *
import drawing
from paint import Painter, grey_palette, colour_palette
from backdrop_util import sun_side_of
import stage3_setting as S3

os.makedirs("figs", exist_ok=True)

def painted(P, lines=True, trees=True, marks=True):
    buf = drawing.finish(render_buffers(P), trees=trees, lines=lines)
    pt = Painter(buf, seed=P.seed, marks=marks)
    idx = pt.paint(buf["st"].boulder, buf["cam"].f, sun_side_of(buf), P.bsize, cam=buf["cam"])
    return idx, buf

# stage 3: the setting as value only (no drawing pass, no marks), + three-value notan
for args in [("medium", "left", "10m", 0), ("large", "back", "close", 1), ("medium", "right", "40m", 3)]:
    P = Params(*args); buf = render_buffers(P)
    gi = S3.value_image(buf); q = S3.quantize_global(gi)
    n = "_".join(map(str, args))
    to_image(q, S3.GREY10, 2).save(f"figs/s3_value_{n}.png")
    to_image(np.digitize(gi, [3.2, 5.6]), np.array([[45]*3, [128]*3, [215]*3]), 2).save(f"figs/s3_notan_{n}.png")

# stage 4 vs 5: the same picture, grey then colour; and marks off vs on
P = Params("medium", "left", "close", 0)
idx0, _ = painted(P, marks=False)
idx1, buf = painted(P)
to_image(idx0, grey_palette(), 3).save("figs/s4_nomarks_grey.png")
to_image(idx1, grey_palette(), 3).save("figs/s4_marks_grey.png")
to_image(idx1, colour_palette(), 3).save("figs/s5_colour.png")
# zoomed material crops (native 6x)
def zoom(idx, pal, box, k=6):
    im = to_image(idx, pal, 1).crop(box); return im.resize((im.width * k, im.height * k), Image.NEAREST)
zoom(idx1, colour_palette(), (40, 90, 160, 160)).save("figs/s4_zoom_granite_grass.png")
zoom(idx1, colour_palette(), (160, 30, 300, 80)).save("figs/s4_zoom_water_trees.png")

# line craft: contours as rendered vs re-drawn in clean runs
P = Params("medium", "left", "10m", 0)
a, _ = painted(P, lines=False, trees=False)
b, _ = painted(P, lines=True, trees=True)
zoom(a, colour_palette(), (60, 80, 150, 130)).save("figs/lines_before.png")
zoom(b, colour_palette(), (60, 80, 150, 130)).save("figs/lines_after.png")
zoom(a, colour_palette(), (0, 0, 160, 45)).save("figs/lines_mtn_before.png")
zoom(b, colour_palette(), (0, 0, 160, 45)).save("figs/lines_mtn_after.png")

# honest comparison: Ferrari native crops beside mine at the same pixel scale (3x)
ref = Image.open("ref/native.png").convert("RGB")
mine = to_image(b, colour_palette(), 1)
pairs = [((150, 190, 470, 370), (0, 0, 320, 180))]
r = ref.crop((150, 190, 470, 370)).resize((960, 540), Image.NEAREST)
m = mine.resize((960, 540), Image.NEAREST)
out = Image.new("RGB", (960 * 2 + 10, 540), (20, 18, 24)); out.paste(r, (0, 0)); out.paste(m, (970, 0))
out.save("figs/compare_full.png")
# detail at 6x: his grass & hummock vs my grass & boulder
r2 = ref.crop((0, 300, 110, 360)).resize((660, 360), Image.NEAREST)
m2 = mine.crop((0, 90, 110, 150)).resize((660, 360), Image.NEAREST)
out = Image.new("RGB", (660 * 2 + 10, 360), (20, 18, 24)); out.paste(r2, (0, 0)); out.paste(m2, (670, 0))
out.save("figs/compare_grass.png")
r3 = ref.crop((200, 70, 330, 140)).resize((780, 420), Image.NEAREST)
m3 = mine.crop((20, 0, 150, 70)).resize((780, 420), Image.NEAREST)
out = Image.new("RGB", (780 * 2 + 10, 420), (20, 18, 24)); out.paste(r3, (0, 0)); out.paste(m3, (790, 0))
out.save("figs/compare_mountains.png")
print("figs done")
EOF
nice -n 10 python3 lessons_figs.py 2>&1 | tail -2; cp stage1/forms_6v.png figs/s1_forms.png; cp stage1/forms_lowsun_dirs.png figs/s1_lowsun.png; ls figs

# output:
figs done
compare_full.png
compare_grass.png
compare_mountains.png
lines_after.png
lines_before.png
lines_mtn_after.png
lines_mtn_before.png
s1_forms.png
s1_lowsun.png
s3_notan_large_back_close_1.png
s3_notan_medium_left_10m_0.png
s3_notan_medium_right_40m_3.png
s3_value_large_back_close_1.png
s3_value_medium_left_10m_0.png
s3_value_medium_right_40m_3.png
s4_marks_grey.png
s4_nomarks_grey.png
s4_zoom_granite_grass.png
s4_zoom_water_trees.png
s5_colour.png

cd ~/work/pixelart-studies/pedagogy/atelier && python3 - <<'EOF'
p='paint.py'; s=open(p).read()
s=s.replace("dark = 1.6 * (1 - np.exp(-cos_i / 0.16))                     # steps of darkening","dark = 2.6 * (1 - np.exp(-cos_i / 0.22))                     # steps of darkening")
# denser, higher-contrast grass strokes (Ferrari: ~40% dark strokes in lit grass)
s=s.replace("""                r = rng.random()
                o = -1.0 if r < 0.38 else (0.9 if r > 0.82 else 0.0)""","""                r = rng.random()
                o = -1.6 if r < 0.30 else (-0.9 if r < 0.46 else (0.9 if r > 0.84 else 0.0))""")
s=s.replace("""        self.step[gr] += o2[gr] * np.clip(blade[gr] / 2.5, 0.45, 1.0)""","""        self.step[gr] += o2[gr] * np.clip(blade[gr] / 2.0, 0.55, 1.0)""")
open(p,'w').write(s)
EOF
nice -n 10 python3 method.py one --size medium --light left --dist 10m --seed 0 >/dev/null; nice -n 10 python3 method.py one --size medium --light left --dist close --seed 0 >/dev/null; echo

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/pedagogy/atelier && rm -rf sheets/cells && nice -n 10 python3 method.py sheets 2>&1 | tail -1 && nice -n 10 python3 lessons_figs.py | tail -1; ls sheets; du -sh sheets figs

# output:
sheets done
figs done
cells
palette_card.png
sheet_light.png
sheet_light_distance.png
sheet_light_distance_preview.png
sheet_light_preview.png
sheet_seed.png
sheet_seed_close.png
sheet_seed_close_preview.png
sheet_seed_preview.png
sheet_size_distance.png
sheet_size_distance_preview.png
4.6M    sheets
336K    figs

# Lessons from the atelier: a boulder by a tarn at dusk, after Mark Ferrari

I worked value first and colour last. Each stage below says what I learned, points to the images, and says what went wrong. After the stages come the method (the program) and an honest comparison of my sheets with Ferrari's tarn.

Everything is in `~/work/pixelart-studies/pedagogy/atelier/`. To regenerate every sheet, run `nice -n 10 python3 method.py sheets`: four single-thread worker processes, about 4 minutes.

**Research, and its limits.** I studied Ferrari's `ferrari_alpine_tarn_dusk.png` closely:

- palette statistics;
- crops at his native pixel scale (the file is his 634×469 frame doubled);
- measured run lengths of his contours.

I read his Q&A (EffectGames) on:

- palette economy;
- patterned dither for smooth ramps and diffusion dither for chaos;
- colour cycling.

On form and light I used Proko's value families (the terminator, core, reflected light that never reaches the halftone, cast shadow, occlusion). I also read Pedro Medeiros and Slynyrd on clustering and on rocks, and the notes on Kazuo Oga's poster-colour method (big soft masses first, detail last, about 21 colours).

I did **not** look closely at Ferrari's other colour-cycling scenes or at Ghibli plates. Saving images from the web counts as a download, which needs permission I could not ask for in this session, so I only had text descriptions of them. That gap is real and shows in the results (see the comparison).

---

## Stage 1: light on simple forms (greyscale)

![forms](figs/s1_forms.png)
![low sun](figs/s1_lowsun.png)

The top row of the first sheet is the two-value **notan** (light family against shadow family). The bottom row is the six-value painting.

What I learned:

- **The two families must never touch.** Early versions let ambient occlusion push lit ground into the shadow family. The notan showed a black pool round every base, and the forms stopped reading. The fix is a rule, not a tweak: occlusion may dim the light family, but it only crosses into shadow at a thin contact line. I found the same bug a second time, in its SDF normalisation, when the terrain distance field was scaled. The notan caught it both times.
- **A value plan beats a lighting formula.** Mapping N·L linearly onto greys gave the lit face of the cube the same grey as the ground. What worked was a six-step plan:
  - light: steps 5 and 4;
  - halftone: 3;
  - *(the gap at the terminator)*;
  - reflected light: 2;
  - shadow: 1;
  - core, occlusion and the root of the cast shadow: 0.

  The greys are spaced by eye (0.07, 0.17, 0.27 | 0.52, 0.70, 0.90), with a big jump at the terminator. A *keyed* material like the ground gets its own local value (an albedo) so it turns fewer steps.
- **The core falls out of geometry.** Point the bounce light from below and away from the sun, and the band just past the terminator faces neither light nor bounce. The cylinder then shows light, halftone, terminator, a dark core, a reflected-light strip and the rim, as in any Bargue-style cast drawing.
- **Dither is for turning surfaces, not flat ones.** A flat plane that lands between two greys got a checkerboard, which was ugly. `quantize(flat_snap=…)` rounds wherever the value is not changing, ignoring jumps at plane edges. Curved forms get one checker band per step, which is how Ferrari's hummocks and snowfields turn.
- **A low sun from each direction** (second sheet), the lesson I carried into the task:
  - from the side, the form is split in two;
  - from the front, the form flattens and survives only in the halftone at its edges;
  - from behind, everything sits in the shadow family and only rims can speak (the rims needed a later stage).

## Stage 2: the boulder as form (greyscale)

![first boulders](stage2/boulder_v2_all.png)
![hull boulders](stage2/hull_all.png)

First attempt: a lopsided ellipsoid with planes shaved off, rounded edges and a joint crack. At a distance the notan read. In close-up it read as a **pillow or loaf**: the flat faces and even rounding looked manufactured. One seed read as a skull (two lobes plus a crack).

Second attempt, `HullBoulder`: I stopped drawing an ellipsoid and started drawing *a potato with facets*.

- It takes a dozen points on a lopsided mass, plus a ring of points at the buried waist so the foot is broad (it *sits*, it does not balance).
- The convex hull of those points gives 20–30 facets of **very unequal size**: a few big planes and many small ones.
- A log-sum-exp blend rounds the edges (weathering, since granite exfoliates in shells).
- A gentle low-frequency swell, one or two shallow spalls (concave scoops) and one partial joint finish it.

Sunk about 40% into the ground, it reads as heavy in five or six values: a lit plane, a halftone top, a clear terminator ridge and a quiet shadow.

What I learned:

- **Keep the shadows simple and put the detail in the light.** My first shadow side had reflected-light patches and a core on every facet that faced slightly away. The result was busy and face-like. Weakening the bounce and narrowing the core made the shadow read as one plane.
- **Where the boulder meets the ground decides whether it is heavy.** An ellipsoid's bulge above the ground makes a black occlusion ring, like a rubber base. Sinking the boulder past its widest point, and later growing grass over the foot, turns that ring into a *seat*.

## Stage 3: the setting as form (greyscale)

![value](figs/s3_value_medium_left_10m_0.png) ![notan](figs/s3_notan_medium_left_10m_0.png)
![value](figs/s3_value_large_back_close_1.png) ![notan](figs/s3_notan_large_back_close_1.png)

The first layout put land behind the boulder and water in front of the camera. I had to work out (on a top-down debug map) where the water must be for the boulder's **reflection** to be visible: *between* the camera and the boulder's foot. The final layout is a cove:

- the boulder sits at the tip of a promontory, on the turf lip;
- its front shore curls round to join the bank the camera stands on;
- its far shore recedes back-left;
- the far bank, spruce and mountains lie 150–200 m off.

Being that far away, they live in a *direction-only* backdrop, so their reflection is automatically the mirror image about the horizon, as it is in nature.

What I learned:

- **Design the notan of the whole picture; don't let the renderer choose it.** In the first honest notan, land and water were the same middle value and the picture could not be read. Ferrari's pattern is:
  - light sky and water;
  - a dark band of trees and far bank;
  - middle-value grass;
  - the focal boulder carrying the full range.

  I got there with value decisions, not physics:
  - weaker Fresnel darkening, so the water stays light;
  - grass *keyed*: flat turf sits on the halftone at any sun height, and slopes climb or fall from there (vertical blades catch a low sun that a flat surface would barely get);
  - grass shadows lifted to a mid-dark, because sky fill is strong at dusk.
- **Contrast closes up with distance.** At 40 m the hummock shadows made dark blotches. Aerial perspective, and shadows that lift toward the light with distance, fixed it.
- **Choices I made that the brief did not state:**
  - pitch 9° rather than exactly 10°, which leaves a strip of sky above the far trees;
  - sun elevation 9°;
  - horizontal field of view 52°;
  - "left" means azimuth 22° (side light from a little in front), because at 8° the camera saw almost only the shadow side;
  - "close" means about 2.7 boulder widths, minimum 3.6 m.

## Stage 4: material, still in grey, and line craft

![no marks](figs/s4_nomarks_grey.png)
![marks](figs/s4_marks_grey.png)

Painting happens in *step space*. Every material owns a ramp, and marks are offsets along it. Quantising mixes only two **adjacent** entries of **one** ramp, with a checker dither. Marks that are meant to be single pixels (granite flecks, glints, tufts) go on *after* the stray-pixel cleanup.

- **Granite.** First I tried salt-and-pepper flecks everywhere, at even density. It read as terrazzo or concrete. What read as granite:
  - flecks sparse and gathered in patches (a noise field decides where grains gather);
  - flecks mostly dark, in the light and halftone, and rare in shadow;
  - a mottle that sticks to the stone (object space, not screen space);
  - convex ridges between planes catching light (+1 step);
  - creases and the shadow-side silhouette going dark;
  - under back light, a 2-px rim that thins as the form turns away.
- **Grass.** Vertical strokes whose length is a fixed 20 cm in the world, converted to pixels by depth: long in front, dithered specks at 40 m. The first stroke value mix was too soft and stripy next to Ferrari (his lit grass is about 40% dark olive strokes on ochre), so I made it denser and higher in contrast. Tufts were first a picket fence of identical sticks. They became **clumps of 3–6 blades fanning from one foot**, each blade a clean 1:2, 1:3 or 1:4 line, tallest in the middle, with lit tips. Grass round the boulder's foot takes its light from the stone behind it.
- **Water.** A broken reflection (rows wobble sideways, more toward the viewer), sparse ripple dashes that thin toward the far shore, one or two catspaw bands, and a dark line under every bank. Under back light, a glitter path of glint pixels widens toward the viewer. The glint colour is reserved for the sun.

![granite and grass](figs/s4_zoom_granite_grass.png)
![water and trees](figs/s4_zoom_water_trees.png)

### Line craft (the critic's note)

The critic pointed out that the pixel medium favours short straight runs and clean 45° and 2:1 steps. I measured Ferrari's mountain skyline at native scale:

- vertical steps are almost always **1 px**;
- horizontal runs **hold steady along a stretch** (2,2,2,2 is a 2:1 line; 1,1,1 is 45°; 4,4,4 is 4:1);
- runs change **gradually** round a curve (6,5,4,4,4).

My ray-marched edges gave noisy runs like 1,3,1,2. `lines.py` redraws every major contour: the bank against the water, the boulder, their reflections, the skyline and the snowline. It smooths the outline, simplifies it to a polygon and snaps each edge to a clean ratio (flat, 4:1, 3:1, 2:1, 1:1, and the steep equivalents), always aiming at the original next vertex so drift cannot build up. It then re-rasterises and refills the changed pixels from their nearest neighbour on the correct side. Snapping can fold an edge back on itself, so one-pixel slivers are removed.

![before](figs/lines_mtn_before.png)
![after](figs/lines_mtn_after.png)
![boulder before](figs/lines_before.png) ![boulder after](figs/lines_after.png)

The note also changed how I *build* things, not only how I clean them:

- **Mountains** are now a few massifs with straight flanks (2:1, 1:1, 3:1 slopes). Their facets fan from each summit along clean ratios, and the snowline is a triangle wave, so the couloirs are clean zigzags. My noise-based mountains before this looked like curtains, scratches and, at one point, a city skyline.
- **Spruce** are built from pixel strokes, the way his are: a 1-px leader with ticks at the tip, then tiers of drooping branch strokes (a flat run that steps down once at its end) round a solid core, with lit branch ends on the sun side.

## Stage 5: colour last

![palette](sheets/palette_card.png)
![grey](figs/s4_marks_grey.png)
![colour](figs/s5_colour.png)

The palette has 32 hand-placed colours (`palette.py`). Each colour was chosen to hit the luminance its ramp position had already proved in grey; the palette card shows each colour above its own grey. Ramps are lists of palette indices, so they can **share**: the water ramp borrows the rock's and mountains' violets and the sky's lilacs, which is honest, since the water is mostly reflecting them. That sharing freed room for a deep peaty violet and one reserved glint.

The dusk logic, taken from Ferrari's own colours (sampled from the tarn painting):

- **Light is warm:** peach, orange, ochre, cream.
- **Shadow is cool and violet**, because the pink and violet sky fills it.
- **Greens are pushed to olive and ochre**, with near-black olive in the darkest turf.
- **Snow in shade uses the sky's lilacs; snow in sun, the sky's creams.** Snow is a mirror of the sky; giving it the mountain ramp gave rose-coloured snow in shade, which was wrong.
- **Trees go to near-black brown-olive**, with bronze lit needles.

Colour changed no value decision. The grey and colour versions have the same structure, which was the point of the order.

---

## The method

`method.py` paints the task from parameters:


python3 method.py one --size {small,medium,large} --light {left,front,right,back} \
                      --dist {close,10m,40m} --seed N [--grey]
python3 method.py sheets



| stage | module | what it does |
|---|---|---|
| geometry | `geometry.py` | `HullBoulder` (lopsided hull, weathered edges, spalls, joint, sunk 38%); `Setting` (cove shoreline in the boulder's view frame, hummocks, turf lip, near bank) |
| value | `valuemodel.py`, `scene.py` | numpy SDF ray-march; six-step value plan (families kept apart, core, reflected light, contact occlusion, keyed grass); water = reflection × Fresnel; backdrop of sky gradient + sun glow, faceted massifs, far bank, all direction-only |
| drawing | `drawing.py`, `lines.py` | contours snapped to clean pixel runs; spruce built from strokes; tree reflections mirrored about the horizon |
| material | `paint.py` | ramps in step space; granite ridges/creases/mottle/flecks, contre-jour rim; grass strokes, fanning tufts, backlit glow; water wobble, ripples, catspaw, bank line, glitter; lit tree edges; stray-pixel cleanup before deliberate single pixels |
| colour | `palette.py` | 32 colours on the proven values; shared ramps |

The parameter sheets (in `sheets/`; each has a `_preview` at 1/3 size, and every cell is also saved alone in `sheets/cells/` as a 3× PNG and a `.npy` index map):

- `sheet_light.png`: light left, front, right, back × 2 seeds (medium, 10 m);
- `sheet_size_distance.png`: size (rows) × distance (columns);
- `sheet_light_distance.png`: light (rows) × distance (columns);
- `sheet_seed.png`: 8 seeds; `sheet_seed_close.png`: 4 seeds, large, close, back-lit.

Every cell is 320×180 at native size, shown at 3× nearest neighbour, with at most 32 colours by construction (every pixel is a palette index). It takes about 15–25 s per cell on one thread.

![light](sheets/sheet_light_preview.png)
![size x distance](sheets/sheet_size_distance_preview.png)
![light x distance](sheets/sheet_light_distance_preview.png)
![seeds](sheets/sheet_seed_preview.png)
![seeds, back-lit close](sheets/sheet_seed_close_preview.png)

---

## Honest comparison with the reference

![whole](figs/compare_full.png)
![grass](figs/compare_grass.png)
![mountains](figs/compare_mountains.png)

Ferrari is on the left, my sheet cell on the right, both at the same pixel scale.

**What holds up:**

- The value architecture reads in every cell and under every light: a light sky and water, a dark tree band, middle-value grass and a focal boulder. Back-lit cells work: glowing sky, dark silhouettes, a cream rim on the stone, glowing grass and a glitter path.
- The boulder reads as heavy, faceted and seated, and the granite reads as granite at close and 10 m range. The reflection hangs beneath it.
- The palette is warm light over violet shadow, with olive and ochre grass, and is recognisably from the same family as his.
- The contours are drawn in clean runs, the spruce are built like his, and the mountain flanks step in 1:1 and 2:1.

**Where it falls short, plainly:**

1. **Texture density.** Ferrari's surfaces are textured *everywhere*: dither is a texture in its own right, not only a transition, and the grass, trees and snow all vibrate. My large areas (the middle water, the lit mountain faces, much of the sky) are flat fields with a transition band. His picture is richer by a wide margin; mine reads cleaner and more graphic, closer to a modern pixel-art poster than to his painterly density.
2. **Water.** His water is the richest part of the tarn painting: a full, detailed, slightly darkened mirror of mountains and trees. Mine reflects correctly, but at 9° down most of my water mirrors the plain upper sky, so a lot of the frame is lilac with ripple dashes. The Fresnel gradient and catspaws help but do not solve it.
3. **Hummocks.** His banks are big rounded mounds that turn from lit ochre to olive shadow through long checker gradients (the island in the tarn). My bank's hummocks are low and read only as texture. I never got a mound to read the way his do.
4. **Trees.** Mine are clean sprites; his have dense, irregular needle texture with bronze light inside the mass, not only at the edges.
5. **Composition.** Every cell uses the same cove recipe, so the sheets are consistent but samey. At 40 m the boulder is honestly small (about 12 px) and the land spit dominates.
6. **Research gap.** I studied one Ferrari painting in depth and knew his other scenes and Ghibli backgrounds only from text. The lessons above are learned against that one painting.

If I had another week, I would put it into items 1 to 3, in that order: dither as texture across the flat masses, richer reflections (lower the water's share of the frame, or give the sky more structure to mirror), and mounds built on purpose as forms in the bank.

cd ~/work/pixelart-studies/pedagogy/atelier && grep -o '([a-z0-9_/]*\.png)' LESSONS.md | tr -d '()' | while read f; do [ -f "$f" ] || echo MISSING $f; done; echo checked; nice -n 10 python3 -c "
import numpy as np, glob
for f in glob.glob('sheets/cells/*.npy')[:40]:
    a=np.load(f); assert a.max()<32
print('all cells <=32 colours')"

# output:
checked
all cells <=32 colours