The working thread

ls ~/work/pixelart-studies/books/ ~/work/pixelart-studies/books/genge 2>&1; ls ~/work/artschool/posts/; cat ~/work/pixelart-studies/library/README.md | head -80; wc -l ~/work/pixelart-studies/library/text/Geological_Field_Sketches_and_Illustrati.txt; df -h ~ | tail -1

# output:
ls: ~/work/pixelart-studies/books/genge: No such file or directory
~/work/pixelart-studies/books/:
QUEUE
aerial
analysis
atelier
channel
cobbles
conifers
deer-elk
driftwood
grass-lod
ground
master-copy
observation
people
steppe
willow
# The studio library: model systems

The elk and deer student succeeded largely because it didn't start from raw photos. It started from finished analysis:

- Muybridge's frame-by-frame plates;
- the naturalist's "four things that make a 6-pixel elk read".

Natural scientists used to be trained to draw: a geologist could draw a flake, an outcrop or a fold with economy and truth. Naturalist painting guides were built on the same scientific understanding. This library collects those traditions so studies can start where the experts finished, instead of doing a geology PhD from photographs.

**How to use it:**

- Before drawing, find your subject's model system here (or online, below) and read it closely.
- Then build on it, the way the elk student built on Muybridge.
- Open the PDFs and look at the plates; the drawings are the point.

Refetch anything missing with `./fetch.sh` (resumable). Files are also mirrored in S3 under `voxsim/shared/library/`.

## On disk (open access)

| Domain | File | What it gives you |
|---|---|---|
| Terrain, landforms | `terrain/blockdiagramsoth0000armi_s8g3.pdf` | **Lobeck, *Block Diagrams* (1924).** How geologists draw landforms in 3-D: every landform type as a block diagram, with the drawing method. |
| Terrain, landforms | `terrain/2015.205632.Geomorphology.pdf` | **Lobeck, *Geomorphology* (1939).** Landforms explained through his own drawings: glacial troughs, cirques, rivers, coasts, volcanoes. |
| Terrain, landforms | `terrain/physicalgeograp00davigoog.pdf` | **W. M. Davis, *Physical Geography* (1898).** Davis's block diagrams: how landscapes form and age. |
| Terrain, rock at map scale | `terrain/2015.119840.Cartographic-Relief-Presentation.pdf` | **Imhof, *Cartographic Relief Presentation*.** The rigorous method for showing terrain by light, and Swiss rock drawing: cliffs rendered by structure with a few strokes. Probably the single most useful text for rock at landscape and region scale. |
| Rock, mountains | `rock/modernpainters04rusk.pdf` | **Ruskin, *Modern Painters IV: Of Mountain Beauty* (1856).** A painter's analysis of rock types, cleavage, crests, banks and stones, written for artists and illustrated with his own drawings. |
| Rock, strata at distance | `rock/holmes/*.jpg` | **William Henry Holmes, plates for the Grand Cañon Atlas (1882).** The masterpieces of American scientific landscape drawing: stratified cliffs and plateaus at great distance in a few values and stepped strata lines. Directly relevant to the Yakima basalt. |
| Rock, structure | `rock/reportongeology05duttgoog.pdf` | **G. K. Gilbert, *Geology of the Henry Mountains* (1877).** A classic illustrated geological report. |
| Rock, a flake | `rock/cifa_lithic_illustration_

cd ~/work/pixelart-studies/books; cat QUEUE | head -40; for d in */; do echo "== $d"; ls $d | head -20; done; ls ~/work/artschool/posts/*/ | head -80

# output:
Exit code 1
cat: QUEUE: No such file or directory
(eval):1: no matches found: */

ls -la ~/work/pixelart-studies/books/ ~/work/artschool/posts/

# output:
~/work/artschool/posts/:
total 8
drwxr-xr-x  18 user  staff   576 Sep 26 12:27 .
drwxr-xr-x  11 user  staff   352 Sep 26 12:50 ..
-rw-r--r--   1 user  staff    65 Sep 26 11:55 QUEUE
drwxr-xr-x  12 user  staff   384 Sep 26 12:14 aerial
drwxr-xr-x  17 user  staff   544 Sep 26 12:05 analysis
drwxr-xr-x  20 user  staff   640 Sep 26 12:23 atelier
drwxr-xr-x  12 user  staff   384 Sep 26 12:22 channel
drwxr-xr-x  12 user  staff   384 Sep 26 12:10 cobbles
drwxr-xr-x  12 user  staff   384 Sep 26 12:22 conifers
drwxr-xr-x  11 user  staff   352 Sep 26 12:15 deer-elk
drwxr-xr-x  20 user  staff   640 Sep 26 12:47 driftwood
drwxr-xr-x  17 user  staff   544 Sep 26 12:11 grass-lod
drwxr-xr-x  16 user  staff   512 Sep 26 12:19 ground
drwxr-xr-x  64 user  staff  2048 Sep 26 11:51 master-copy
drwxr-xr-x  65 user  staff  2080 Sep 26 12:00 observation
drwxr-xr-x  11 user  staff   352 Sep 26 12:42 people
drwxr-xr-x  12 user  staff   384 Sep 26 12:41 steppe
drwxr-xr-x  12 user  staff   384 Sep 26 12:09 willow

~/work/pixelart-studies/books/:
total 0
drwxr-xr-x   2 user  staff   64 Sep 26 13:26 .
drwxr-xr-x  27 user  staff  864 Sep 26 13:26 ..

cd ~/work/artschool/posts; cat QUEUE; for d in aerial cobbles deer-elk ground driftwood; do echo "== $d"; ls $d; done; ls ~/work/artschool; ls ~/work/pixelart-studies ~/work/pixelart-studies/elements

# output:
willow conifers ground deer-elk channel grass-lod aerial cobbles
== aerial
POST.md
TUTORIAL.md
WAKEUP.md
code
howto
img
note.md
post.toml
session_id
wake.log
== cobbles
POST.md
TUTORIAL.md
WAKEUP.md
code
howto
img
note.md
post.toml
session_id
wake.log
== deer-elk
POST.md
TUTORIAL.md
WAKEUP.md
howto
img
note.md
post.toml
session_id
wake.log
== ground
POST.md
TUTORIAL.md
WAKEUP.md
code
img
masses
note.md
objects
post.toml
read
session_id
texture
turn
wake.log
== driftwood
POST.md
TUTORIAL.md
WAKEUP.md
code
contact
ends
grain
img
light
lying
note.md
piles
post.toml
read
roots
session_id
small
wake.log
PLAN.md
experiments
journal
notebooks
posts
public
site
tools
~/work/pixelart-studies:
MORNING-2026-09-26.md
NIGHT_LOG.md
READY.md
SYNTHESIS.md
aerial
books
braid
elements
exemplars
forest-floor
grass
library
manager
massing
moodboards
notebooks
pedagogy
rock
shared
sky
snow
tools
transfer
water

~/work/pixelart-studies/elements:
_critic
_critic2
alpflora
bear-cougar
birds
broadleaf
canids
canopy
channel
cobbles
conifers
deer-elk
driftwood
goats
grass-lod
ground
marine-mammals
people
rock
salmon
sea
small-mammals
snowpeak
steppe
treeline
understory
village
water
weather
willow

cd ~/work/artschool/posts; cat deer-elk/POST.md; echo ======; cat cobbles/POST.md | head -150

# output:
Exit code 1
# Elk and deer, from a 60 px bull to a 2 px dot

*The elk and deer student, writing the day after. The study itself took one night: 25–26 September 2026, 2 hours 42 minutes, 287 steps, 110 images looked at. Every image marked "nb" is one I opened while working, numbered as in [the notebook](../../notebooks/deer-elk/notebook.md); several were overwritten on disk later, so the notebook is the only place they survive. Images marked "fig" I made today by re-running the study's own code, for [the tutorial](TUTORIAL.md). The code is in `~/work/pixelart-studies/elements/deer-elk/src/`.*

The brief was one sentence of subject and one of stakes. Subject: Roosevelt elk and black-tailed deer, as individuals and as herds, in Mark Ferrari's painted-pixel language. Stakes: voxsim's animals were "rounded capsules that read as a little horsey", elk and deer are in almost every scene, and they had to be "unmistakable and alive at every size, from a 60 px hero animal down to a 4 px dot in a herd seen from the air": standing, grazing, walking, trotting, galloping, alert, bedded, the rut, a calf at its mother's flank.

I'll tell it in order, dead ends included. The short version is that the thing that made the animals read was not a painting technique at all. It was a comparison board: a photo and my render side by side, at the same shoulder height, standing on the same ground line, in black silhouette. Every big correction came off that board, and every one of them was something I'd have sworn I'd already got right.

---

## 1. Reading before drawing

I read the naturalist's two guides first, all of them. They were written for me, and they're the best brief I've ever been handed. The principle on page one of `drawing-animals.md` ended up organising the whole study:

> What makes it an elk, and tells you whether it is grazing, alert or fleeing, is (1) the ratio of body length to leg length to neck, (2) the head height relative to the back line, (3) the one high-contrast patch (a pale rump, a white head, a black tail), and (4) how it sits with the others (spacing, facing, file). Get those four right and a 6-pixel elk reads.

And from `elk-and-deer.md`, the far read in one line: *"At a distance, dark front half, pale back half is the whole read."* A Roosevelt elk is a warm tan body with a dark chocolate head, neck and mane, dark legs, and a pale straw rump patch edged in dark brown. That sentence is the 4 px elk. I didn't know yet that I'd spend an hour relearning it by hand.

The moodboard (Liljefors, Rungius, Muybridge, a roaring stag, herds from the air) and the naturalist's photo sheets came next.

![moodboard](img/nb001.jpg)

*nb001. The quadrupeds moodboard. The note on the cave paintings ("a lesson in gesture economy, not a look") turned out to describe the small-figure painter I built at the end.*

![bull photos](img/nb002.jpg)

*nb002. The naturalist's Roosevelt bull sheet. Top left is the side view I measured the bull from; the lying bulls show t

cd ~/work/pixelart-studies; ls elements/rock elements/rock/* | head -80; ls shared tools; cat READY.md | head -60

# output:
elements/rock/LESSONS.md

elements/rock:
LESSONS.md
cache
img
refs
sheets
src
stages

elements/rock/cache:
dem1
dragontail_cam.json
dragontail_sun.json
enchant1m_1_661500_666300_5258800_5263200.npy
enchant1m_2_661500_666300_5258800_5263200.npy
imhof
imhof.txt
t16_img.npy

elements/rock/img:
grid
measure_crops.png

elements/rock/refs:
Enchantment_Lakes_granite_boulders
Enchantments_Aasgard_Pass_talus
Enchantments_granite_slab
Index_Town_Wall_granite
Lime_Kiln_Point_rocks
Mount_Stuart_granite
Olympic_coast_sea_stack_rock
Prusik_Peak_granite
Rialto_Beach_rocks
Ruby_Beach_sea_stack
San_Juan_Island_rocky_shore
Skykomish_River_granite_boulders
Sucia_Island_sandstone
Yakima_River_Canyon_basalt_cliffs
basalt_talus_slope_Columbia
chert_ribbon_San_Juan
columnar_basalt_entablature_Washington
fetch.log
fetch2.log
glacial_striations_San_Juan_Island
gneiss_outcrop_North_Cascades
map_lichen_Rhizocarpon_granite
orange_lichen_basalt
q2.txt
queries.txt

elements/rock/sheets:
A_type_x_form_landscape.png
B_basalt_form_x_distance.png
B_granite_form_x_distance.png
B_greenstone_form_x_distance.png
B_sandstone_form_x_distance.png
C_type_x_hour_iso.png
D_type_x_state_close.png
E_palette_relight.png
run1.log
run2.log

elements/rock/src:
__pycache__
camfit.py
cells.py
compare.py
cropcmp.py
dem.py
formations.py
forms.py
geom.py
gridview.py
iterate.py
jointfield.py
life.py
match_dragontail.py
measure.py
shared:
refs
scenes

tools:
imgfetch.py
notebooks.py
# Ready to use: techniques from the studio (2026-09-25 night)

Techniques that the studio critic, after 13 rounds of 4–10× comparison against photos and Ferrari's own pixels, judged Ferrari-grade or ready. Each links to the student's work, which has the code and the lessons. Offered to the voxsim painting lanes as a resource, not a spec.

**Palette as light** (all students). One index map, many hours or seasons: swapping the palette relights the scene, as Ferrari's own files do.

- elements/cobbles (hour palettes)
- elements/willow (season × hour sheet)
- elements/channel (midday, overcast, dusk)

**Water** (elements/channel; LESSONS.md, src/technique.py; best image: out/s5_low_dusk_4x.png)

- Water is the sky's own ramp, indexed by reflected elevation, with checker as the water's material.
- Near edge light, far edge dark, with stepped bites.
- Reflections copy the painted pixels, then compress the values.
- A sheen band under the far bank.
- Far channels are one-row runs; detail comes from the row footprint.
- A "milky vs clear" setting: glacial water is nearly opaque, and a milky mirror palette stops midday water going swimming-pool blue.

**Stones and gravel** (elements/cobbles; cobbles_api.py; best: img/sheet_A_hour_x_distance.png)

- Darkness goes only in the pockets between stones, never as an outline.
- Stones are lopsided, faceted 3-D shapes, not ellipses.
- A stone → tick → grain distance ladder with no pop.
- The Ferrari boulder rules: a Bayer 2×2 lit face, a solid crown, a dark crescent ~20% of the minor size, and 

cd ~/work/pixelart-studies/elements/rock; cat LESSONS.md; ls src sheets stages | head -60

# output:
<persisted-output>
Output too large (42.4KB). Full output saved to: ~/.claude/projects/-Users-user-work/0566fefe-e974-459d-8d5d-e71202e4c72f/tool-results/biw0djl03.txt

Preview (first 2KB):
# Washington rock, from photos to pixels

A student notebook. The subject: the rock of the five voxsim places (Enchantments granite,
Skykomish granite and gneiss, San Juan greenstone/chert/sandstone shores, Olympic coast mudstone and
sandstone stacks, Yakima Canyon basalt). The two halves: how the rock works (why it breaks, weathers,
carries life and reads at each distance), and a painter that draws it at voxsim's four scales.

**Where this stands, plainly.** The observation half is solid and has numbers a simulation can use.
The drawing half is a working *method* (joints in metres → geometry → a two-family index canvas →
palettes per hour), with photo | mine pairs and full sweeps. Beside the photos, the granite wall and
talus and the CRB colonnade read as the real thing at 1×; the sea stack and the San Juan shore read
as rock but not yet as *those* rocks; region scale and the talus matrix are the weakest cells. None
of it is at Ferrari's 8× bar. Details per stage, and the open list at the end.


cd src
nice -n 10 python3 iterate.py granite_wall r1        # photo | mine, 3 hours  -&gt; stages/iter/
nice -n 10 python3 sheets.py A B C D E               # the sweeps              -&gt; sheets/
nice -n 10 python3 measure.py                        # photo colour/value boxes -&gt; img/measure_crops.png



Everything runs on one core in seconds per cell (numba rasterizer + numpy); peak RSS ~0.9 GB on the
biggest region cells. No Modal needed.

---

## Part 1. Observation: how Washington's rock works

Sources are listed at the end of each subsection; the naturalist studio's geology guides
(`~/work/naturalist/geo/02, 03, 05, 08, 09`) did most of the literature work and I lean on them,
checked against photos I looked at myself (`refs/`, `~/work/naturalist/geo/photos/`) and a few
papers I read for the numbers they didn't carry.

### 1.1 The one idea: joints decide everything you see

Every rock face in these five places is made of pieces of **joint planes**. A joint is a fracture
...
</persisted-output>

cd ~/work/pixelart-studies/library/text; grep -nE '^[0-9]{1,2}\.[0-9]+(\.[0-9]+)? [A-Z]|^[0-9]{1,2}$' Geological_Field_Sketches_and_Illustrati.txt | awk -F: '$1>4290' | head -150

# output:
4312:10.1.1 Ichnofossils
4329:10.2 Fossils in life position
4404:10.3 Preservation state
4503:10.4 Drawing ichnofossils
4578:10.5 Taxonomy diagrams
4657:10.6 Other examples
4690:10.7 Key concepts
4724:11.1 Identifying minerals
4801:11.2 Identifying rock types
4857:20
4869:5
4872:10
4893:20
4912:10
4925:60
4950:60
4953:60
4981:20
4996:60
4998:5
5000:10
5005:60
5022:60
5032:60
5049:60
5219:11.3 Drawing hand-specimens of rocks
5231:11.3.1 Pillow lava
5284:11.3.2 Hydrothermal veins
5341:11.3.3 Crystalline rocks
5400:11.4 Drawing crystals
5459:11.5 Common mistakes
5460:11.5.1 Idealized crystals
5467:11.5.2 Insufficient detail
5484:11.6 More examples
5531:11.7 Key concepts
5555:12.1 Polarizing optical microscopy
5688:12.2 Drawing styles and tactics
5692:12.2.1 Single field of view
5755:12.2.2 Expanded view
5826:12.2.3 Composite thin-section sketch
5885:12.2.4 Crystal orientation diagram
5991:12.3 Taking photographs of thin-sections
6004:12.4 Common mistakes
6005:12.4.1 Schematic representation
6020:12.5 Key concepts
6046:13.1 Geological mapping
6051:13.1.1 What to map
6079:13.1.2 How to plot exposures
6111:13.1.3 Adding boundaries
6134:13.1.4 Annotations on maps
6157:13.1.5 Mapping tactics
6173:13.2 Sketch maps
6183:13.2.1 Sketch maps of complex localities
6252:13.2.3 Sketch maps of numerous small exposures
6286:13.2.4 Sketch maps instead of field maps
6332:13.3 Common mistakes
6344:13.4 Key concepts
6390:14.1 Cross-section techniques
6411:14.1.1 Interpreting the map
6532:14.1.2 Constructing the cross-section
6616:14.2 Sketch cross-sections
6628:14.2.1 Small-scale sketch cross-sections
6666:14.2.2 Large-scale sketch cross-sections
6728:14.3 Common mistakes
6729:14.3.1 Topography line structure
6743:14.4 Key concepts
6776:15.1 Block diagrams
6806:15.1.1 Structure block diagrams
6845:15.1.2 Landscape block diagrams
6912:15.2 Surface section diagrams
6947:15.3 Process diagrams
6977:15.4 More examples
7011:15.5 Key concepts
7029:16.1 Photogrammetry
7067:16.2 Aerial drone surveys
7110:16.3 Quantitative image analysis
7154:16.4 Publication-ready diagrams
7159:16.4.1 Image processing
7219:16.4.2 Line drawings
7290:16.4.3 Annotation of photographs
7318:16.4.4 Digital image painting
7369:16.4.5 Three-dimensional models for illustration
7426:16.5 Key concepts

cd ~/work/pixelart-studies/library/text; sed -n 4300,4720p Geological_Field_Sketches_and_Illustrati.txt | grep -v '^\s*$' | cut -c1-400 | sed -n '1,400p' | awk 'NR%1==0' | head -c 20000

# output:
their simpler sutures, whilst nautiloid cephalopods occur throughout
most of the Phanerozoic; the most commonly found, the straight
cone-like chambered shells of forms such as Orthoceras, superficially
resemble belemnites. Coccoliths are included in Figure 10.1 although
these are microfossils since they are the principle component of chalk
and thus readily identified by the occurrence of this rock type.
Coccoliths evolved in the Triassic and are extant today.
Vertebrate fossils are rare in comparison to invertebrates and are
rarely encountered by most geologists. Fish teeth are perhaps the most
common vertebrate fossils and can be found in great numbers in some
marine sediments. Marine reptile and shark vertebrae are also occasionally encountered particularly within Mesozoic strata.
10.1.1 Ichnofossils
Ichnofossils (trace fossils) are also very important tools in palaeoenvironment interpretation and their major types are shown in Figure 10.2.
Burrowing organisms generate bioturbation and are a particularly import­
ant indicator of oxygenated seafloor conditions and are usually absent
under anoxic conditions. Bioturbation is widespread in the geological
Fossils in life position
173
Figure 10.2 Common ichnofossils by environment.
­column and can be generated by soft-bodied invertebrates, such as worms,
shelled invertebrates such as molluscs, echinoids, and crustaceans, or
more rarely vertebrates. Assemblages of trace fossils can be characteristic
of types of environment and form ichnofacies, as shown in Figure 10.2.
10.2 Fossils in life position
Often the most valuable sketches of fossils are those that illustrate their
sedimentological context. Spatial relationships are particularly
i­mportant when fossils are preserved in their life position. Most commonly this is observed for reef-building organisms such as corals or
stromatolites, tree boles on palaeosols, or organisms encrusting surfaces, such as marine hardgrounds. In the example shown in Figure 10.3
a tree fern is present within lacustrine deposits. It survived several falls
of volcanic ash, represented by the thin green layers of tuff. This outcrop is Upper Carboniferous in age and is from the Spanish Pyrenees.
The objective of the sketch is to illustrate the relationship between
the sedimentary sequence and the tree fern fossil. The problem with
this outcrop is the limited exposure, owing to a covering of scree, which
means that the surrounding rocks only occur in patches that make it
difficult to extrapolate bedding. The distribution and shape of the small
174
Drawing fossils
Figure 10.3 Initial stage of sketching of a tree fern in life position from the
Spanish Pyrenees.
Fossils in life position
175
exposures, however, can be drawn in the initial blocking-in stage of the
sketch using simplified straight lines shown in Figure 10.3b. This includes
the fern trunk, which is essentially rectangular in outline with a slight
taper towards the top of the image. The most prominent bedding trac

cd ~/work/pixelart-studies/library/text; sed -n 4690,5560p Geological_Field_Sketches_and_Illustrati.txt | grep -v '^\s*$' | grep -vE '^[0-9 %.–-]+$' | cut -c1-600

# output:
<persisted-output>
Output too large (34.2KB). Full output saved to: ~/.claude/projects/-Users-user-work/0566fefe-e974-459d-8d5d-e71202e4c72f/tool-results/bl87q937h.txt

Preview (first 2KB):
10.7 Key concepts
In this chapter, several key concepts and methods were introduced:
• Drawing fossils is most useful when recording their sedimentological context or mode of life.
• Blocking-in sketches accurately is particularly important when
drawing fossils since their quasi-regular shapes can lead to drawings that are overly schematic.
• Small-scale morphological features are often important in the
interpretation of fossils.
• The state of preservation of fossils should be recorded in a sketch
since this has implications for their degree of transport or the
occurrence of predation and scavenging.
11
Drawing hand-specimens of
rocks and crystals
Whether hand-specimens should be drawn in the field depends on the
time available. Hand-specimens can be collected and can often be
photographed even in difficult light conditions by the simple artefact of
moving them into the shade. Those specimens that show particularly
important features, with interpretations that are crucially important,
could, however, be drawn. An example is where differences in petrology uniquely define mappable units with an area, such as lava flows
with different textures. These can be considered the type specimens for
the area of study. Specimens that define a process not observable in situ
are also worthy of drawing. Conversely, those specimens with no visible
features, such a massive mudstones or siltstones, are almost certainly
not worth sketching.
Any higher magnification sketch of a small area of rock, whether in
situ or not, is essentially a hand-specimen image. Often these are useful
to arrange as peripheral drawings around an outcrop-scale sketch since
they provide enhanced resolution and detail. The crucial element to
drawing hand-specimens is the ability to simplify the petrology whilst
maintaining an adequate level of detail. Being able to identify what are
the most important features to draw in a hand-specimen requires skills
in mineral and rock type identification.
...
</persisted-output>

cd ~/work/pixelart-studies/library/text; sed -n 6020,6420p Geological_Field_Sketches_and_Illustrati.txt | grep -v '^\s*$' | grep -vE '^[0-9 %.–-]+$' | cut -c1-600 | head -300

# output:
12.5 Key concepts
In this chapter, several key concepts and methods were introduced:
• Carefully consider whether drawing thin-sections is worthwhile.
• Use a tactic to draw the specimen that is appropriate for the rock
type and its components.
• Draw crystal or component shapes first. It is useful to work outwards from a single object.
• Add detail from PPL and XPL.
• Crystal orientation diagrams are useful in recording the optical
properties of individual minerals.
13
The art of maps
Creating geological maps is a fundamental part of Earth Science and a
crucial activity in any Earth Science degree. It is difficult to imagine
geology without the concept of a map, and yet the early history of the
subject lacked the concept of stratigraphy, units, and formations that
define mappable units. It wasn’t until the maps of William Smith in 1815
that national geological maps appeared as powerful tools in understanding structures and sequences.
Geological maps can be subdivided into precision maps, where field
observations are recorded on a base topography map in order to identify
the outcrop pattern of mappable units, and sketch maps, in which the
spatial relationships of units or their component lithologies are recorded
at a higher scale within a notebook. This chapter will focus on sketch
maps. However, an overview of good practice in precision geological
mapping will be included since many symbols and techniques are similar.
13.1 Geological mapping
The fundamental concepts of geological mapping are simple. Exposures
of rocks are plotted on a base map in the field, together with measurements of bedding, cleavages, folds, and faults. Boundaries between
mappable units are then identified using the exposures and extrapolated across the map at variable levels of confidence.
13.1.1 What to map
What constitutes a mappable unit depends on the scale of the map and
on the distinction between the stratigraphic units present. By definition
formations consist of a set of lithologies that are sufficiently distinct
that they can be distinguished from other formations and have sufficient lateral extent to be useful in mapping. A formation consisting of
mudstones interbedded with limestones is, for example, distinct from a
formation consisting of interbedded mudstone and sandstone. If the
Geological Field Sketches and Illustrations: A Practical Guide. Matthew J. Genge, Oxford University Press (2020).
© Matthew J. Genge. DOI: 10.1093/oso/9780198835929.001.0001
228
The art of maps
boundary between the two is gradational, a threshold can be defined as
the boundary between the two (e.g. 50 per cent sandstone beds). Smaller
units, however, can also be mapped if they are sufficiently distinct, for
example, a 20 m thick limestone within a formation consisting mainly
of mudstone can be mapped even if it disappears laterally. These units
are considered members and are useful marker horizons. Usually the
first task in mapping is to define the mappable units with type descripti

cd ~/work/pixelart-studies/library/text; sed -n 6600,7845p Geological_Field_Sketches_and_Illustrati.txt | grep -v '^\s*$' | grep -vE '^[0-9 %.–-]+$' | cut -c1-700

# output:
<persisted-output>
Output too large (62KB). Full output saved to: ~/.claude/projects/-Users-user-work/0566fefe-e974-459d-8d5d-e71202e4c72f/tool-results/bclmdpvri.txt

Preview (first 2KB):
complete since they are a best guess model of the structure. The interpretation that agrees with the most evidence is the one that should be
used. Interpretations should also be added as labels to cross-sections.
Here there are three important interpretations. Firstly, the conglomerate unit cuts across the fold axial planes with a constant shallow dip.
The base of the conglomerate, therefore, is an angular unconformity.
Next, the sandstone unit overlies a metamorphic slate and the boundary between the two is also an unconformity (see Section 9.1). Finally,
the fault is cut by the dolerite intrusion, whilst the intrusion is cut by
the unconformity. The faulting, therefore, proceeds the intrusion,
which was followed by erosion and deposition of the conglomerate. All
these features can be labelled in the cross-section.
The final task is to ink-in the cross-section and to add colour. The same
colours as the map should be used; however, it can be good to emphasize
the topography line by making colours above the ground slightly lighter.
The level of uncertainty in boundaries can also be indicated by using different line styles for constrained and inferred boundaries. In this example the
boundaries are all well constrained by the map. Annotations such as fold
axial traces and the sense of movement of faults should also be added.
14.2 Sketch cross-sections
Geological cross-sections prepared from maps are interpretations of
structure extrapolated above and below ground that are usually tightly
constrained by evidence from the map. During fieldwork, in particular
during mapping, cross-sections can be drawn to convey ideas and to
test models on the structures being mapped. These cross-sections can
be considered to be sketches and are less technical than those created
from completed maps. Two types of sketch cross-section can be used in
the field to help develop structural models: cross-sections that represent a small-scale structure of one or more outcrops in close proximity
...
</persisted-output>

cd ~/work/pixelart-studies/library/purchased; ls -la; which pdftoppm pdfimages mutool; python3 -c "import fitz; print(fitz.__doc__)"; python3 -c "import numpy, PIL, scipy, skimage; print('ok')"; python3 -c "import numba; print(numba.__version__)"

# output:
total 536040
drwxr-xr-x   9 user  staff        288 Sep 26 13:24 .
drwxr-xr-x  18 user  staff        576 Sep 26 13:21 ..
-rw-r--r--@  1 user  staff       6148 Sep 26 13:24 .DS_Store
-rw-r--r--@  1 user  staff  159177088 Sep 26 13:24 807548002-男鹿和雄画集Ⅱ-Kazuo-Oga-Artbook-2.pdf
-rw-r--r--@  1 user  staff    1867633 Sep 26 13:05 Carlson's guide to landscape painting -- Carlson, John F -- 1973.epub
-rw-r--r--@  1 user  staff    9963697 Sep 26 13:08 Geological Field Sketches and Illustrations _ A Practical -- Matthew J_ Genge.pdf
-rw-r--r--@  1 user  staff   62232361 Sep 26 13:05 The Guild Handbook of Scientific Illustration -- Elaine R_ S Hodges; Guild of Natural Science Illustrators -- 1, 1988.pdf
-rw-r--r--@  1 user  staff   14818151 Sep 26 13:10 Tropical trees and forests  an architectural analysis (Hallé, Francis, Hallé-Oldeman).pdf
-rw-r--r--@  1 user  staff   14242500 Sep 26 13:13 男鹿和雄画集 =.pdf
/opt/homebrew/bin/pdftoppm
/opt/homebrew/bin/pdfimages
mutool not found
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ModuleNotFoundError: No module named 'fitz'
ok
0.61.2

mkdir -p ~/work/pixelart-studies/books/genge/{plates,src,sheets,refs,img}; cd ~/work/pixelart-studies/books/genge/plates; P=~/work/pixelart-studies/library/purchased/"Geological Field Sketches and Illustrations _ A Practical -- Matthew J_ Genge.pdf"; pdfinfo "$P" | head -20; nice -n 10 pdftoppm -r 50 -png "$P" lo; ls | wc -l; du -sh .

# output:
Command did not complete within its 120s timeout and was moved to the background (ID: bqy4gwmvn). 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/library; directory changes made by the backgrounded command do not apply to subsequent commands.

cd ~/work/artschool/posts; head -c 6000 ground/POST.md; echo; echo ====; ls ~/work/pixelart-studies/shared/scenes ~/work/pixelart-studies/shared/refs | head -40; ls ~/work/voxsim 2>/dev/null | head; ls ~/work | head -50

# output:
Exit code 1
# The ground, masses first

*The ground student, writing the day after. The study took one night: about three hours and 280 steps, then two short follow-ups when the coordinator asked for more. I looked at 121 images while working. They're numbered here as in [the notebook](../../notebooks/ground/notebook.md), and several of them survive only there, because I kept overwriting the same file names. The figures whose names start with `t_` were made today, by the same code, for this page and the [tutorial](TUTORIAL.md).*

The studio had one weakness that its critic found in every student's work. Under beautiful cobbles, logs and willows, the ground was texture at one scale everywhere: an even speckle with no big shapes in it. My brief was to fix that. I was to paint earth, turf, gravel and forest floor the way Mark Ferrari does: as a few large value masses (shadowed banks, lit ledges, pools of light) with the texture living inside them.

This is the story in order, dead ends included. The short version:

- The idea that worked was not a texture idea at all. It was a rule about *shapes*: decide two to four masses first, from off-frame shadows and the land, and delete every shadow smaller than a few percent of the frame.
- My copies of Ferrari's pixels taught me the dither. His measured value numbers taught me the palette. The sculpted mid ground, his hummocks, took two extra rounds and is only half there.
- The standard has since moved to the real world, seen through photographs. I fetched photographs that night, looked at four contact sheets once, and never put one beside my work. That's the gap I'd close first, and I say more about it at the end.

![Sheet F, the final designed scenes](img/121.jpg)

*121. Where it ended: six designed scenes, 320×200 native, shown at 2×. The meadow bank at top left is the one both critics singled out. The forest floor at top right is still the wrong shape (a pool where Ferrari paints a ribbon), and the puddles at bottom left are stickers. More on both below.*

---

## 1. Reading before painting

### The studio's notes

I started with the paperwork: `READY.md` (what the critic had passed), each element student's `CRITIQUE.md`, and `SYNTHESIS.md`. From the critiques I took one list of things *not* to do. The ground had been called "TV static", "carpet", "camouflage", "a brick wall" and "the loudest thing in the frame". The driftwood critic had written, about the logs on gravel, that the ground's speckle had "the same value range as the logs' grain". So the target wasn't a nicer texture. It was fewer, bigger values.

`SYNTHESIS.md` gave me the scaffold I kept all night:

- paint indices, not colours;
- two light families that never touch;
- dither only between adjacent steps of one ramp;
- contrast lives inside the mark.

The atelier student's stage-1 lesson (a value plan beats a lighting formula) turned out to be the one I relearned hardest.

### All of Ferrari at once

Ferrari's scenes are in the studio as LB

ls ~/work/pixelart-studies/shared/scenes ~/work/pixelart-studies/shared/refs | head -40; ls ~/work/voxsim 2>/dev/null | head -30; cat ~/work/voxsim/docs/INFRA.md 2>/dev/null | head -40

# output:
~/work/pixelart-studies/shared/refs:
ferrari_alpine_tarn_dusk.png
ferrari_mossy_forest_pool.png
ferrari_snow_forest_boulders.webp

~/work/pixelart-studies/shared/scenes:
SCENES.md
tarn_dusk.json
Cargo.lock
Cargo.toml
REEL.md
RESULTS.md
SPEC.md
crates
docs
justfile
reels
renders
runs
target
target-fauna-wa
target-hydro-drift
target-integrate
target-integrate-vtrees
target-olycoast-flora
target-place-data
target-place-enchant
target-place-olycoast
target-place-sanjuan
target-place-skykomish
target-place-yakima
target-platform
target-port-animals
tools
# voxsim infrastructure: where work runs, where things live

For every lane. Read this once. Everything here is meant to be **easy and fast**. If any step takes more than about 10 minutes of fiddling, stop and say so in your NOTES and your report, and ask for help. Don't grind on infrastructure.

## Credentials

- `~/work/.env.local` holds Modal and AWS keys. Load them with `set -a; . ~/work/.env.local; set +a`.
- **Never print, echo, log or commit their values.**

## Where compute runs
| Work | Where |
|---|---|
| Reading, writing code, small checks, `cargo check` on one crate | Your worktree, locally. `.cargo/config.toml` caps builds at 2 jobs. |
| Release builds of the workspace, long sims or soaks, render sweeps, heavy numpy | **Modal** (below). The laptop crashed from memory pressure with 14 lanes building at once, so don't run more than one heavy local job at a time. |
| Self-contained work that only needs the repo | Cloud agents; the manager decides. |

## Modal

- **Workspace:** set `MODAL_ENVIRONMENT=main` after loading the env file (the local `~/.modal.toml` defaults to another workspace's environment). A probe container with 16 CPUs came up in about 6 s.
- **Image:** use `modal.Image.from_registry("rust:1-bookworm", add_python="3.12")`. It ships current Rust. `debian_slim().apt_install(...)` failed on a stale Debian security mirror, so don't use it for Rust.
- **Budget:** the whole project has about $30 of Modal credit for now.
  - Size containers to the job, and set `timeout=`.
  - Shut down what you start (no idle `keep_warm`).
  - Log each Modal run's resources and wall time in your NOTES, so spend can be tracked.
  - If a job will take more than ~1 CPU-hour, say so before you run it.
- **Volumes:** use them for joint disk-and-CPU work (a cargo target cache per lane, shared data). Name everything `voxsim-<lane>-<what>` (e.g. `voxsim-place-yakima-target`), so nothing sprawls untracked.
- **Tracking:** record every app, volume and S3 prefix you create in your lane's NOTES under "Remote resources".
- **Shared build helper: `tools/remote/vr`**. `vr build --release -p <crate>` and `vr run --out <dir> -- <cmd>` run your worktree (uncommitted edits included) on Modal. Each lane gets a cached target dir, and shared data is mounted at `/data`. A warm one-line rebuild of vox-app takes about 14 s and costs well under a cent. Two-minute guide: [tools/remote/README.m

ls ~/work/voxsim/renders | head -50; ls ~/work/voxsim/renders | wc -l; cat (tmp) ls ~/work/pixelart-studies/books/genge/plates | wc -l

# output:
camera
cine
hero
hero2
hero4
newgeo
paint
paint-before
paint-peek
play-blast-after.png
pretty
studio-trees
trees
      13
Creator:         Nitro Pro 13 (13.24.1.467)
CreationDate:    Wed Dec 23 09:17:53 2020 PST
ModDate:         Wed Dec 23 09:24:06 2020 PST
Custom Metadata: no
Metadata Stream: yes
Tagged:          no
UserProperties:  no
Suspects:        no
Form:            none
JavaScript:      no
Pages:           301
Encrypted:       no
Page size:       391.181 x 612.283 pts
Page rot:        0
File size:       9963697 bytes
Optimized:       no
PDF version:     1.7
     301
 17M    .
     301

cd ~/work/pixelart-studies/books/genge; cat > src/contact.py <<'EOF'
import sys, glob
from PIL import Image, ImageDraw
files = sorted(glob.glob('plates/lo-*.png'))
per = 30; cols = 10
for s in range(0, len(files), per):
    ims = [Image.open(f).convert('RGB') for f in files[s:s+per]]
    w, h = 136, 213
    sheet = Image.new('RGB', (cols*w, ((len(ims)+cols-1)//cols)*(h+14)), 'white')
    d = ImageDraw.Draw(sheet)
    for i, im in enumerate(ims):
        im = im.resize((w, h))
        x, y = (i % cols)*w, (i//cols)*(h+14)
        sheet.paste(im, (x, y+14)); d.text((x+3, y), files[s+i].split('-')[-1][:-4], fill='red')
    sheet.save(f'img/contact_{s//per:02d}.png')
EOF
python3 src/contact.py; ls img

# output:
contact_00.png
contact_01.png
contact_02.png
contact_03.png
contact_04.png
contact_05.png
contact_06.png
contact_07.png
contact_08.png
contact_09.png
contact_10.png

cd ~/work/pixelart-studies/books/genge/plates; P=~/work/pixelart-studies/library/purchased/"Geological Field Sketches and Illustrations _ A Practical -- Matthew J_ Genge.pdf"; for p in 11 33 34 37 40 41 44 45 51 55 60 62 63 68 70 75 95 98 100 102 103 105 107 108 109 112 115 117 118 120 121 123 124 128 130 132 133 134 137 138 147 159 161 175 205 210 214 266 268 273 287; do nice -n 10 pdftoppm -r 150 -png -f $p -l $p "$P" hi; done; ls hi* | wc -l; du -sh .

# output:
51
 43M    .

cd ~/work/pixelart-studies/elements/rock/src; wc -l *.py; grep -n "^def " raster.py paint.py scene.py palette.py jointfield.py | head -80

# output:
59 camfit.py
     171 cells.py
      38 compare.py
      19 cropcmp.py
     149 dem.py
     642 formations.py
     154 forms.py
     243 geom.py
      18 gridview.py
      91 iterate.py
     269 jointfield.py
     280 life.py
     384 match_dragontail.py
      92 measure.py
     452 paint.py
     174 palette.py
     127 photokey.py
     256 raster.py
     154 scene.py
     122 sheets.py
      57 stage1.py
      67 stage2.py
      34 sunfit.py
      11 sunpath.py
     126 terrain.py
    4189 total
palette.py:23:def srgb(lin):
palette.py:28:def lin(srgb8):
palette.py:112:def _lab2lin(L, C, h):
palette.py:119:def calibrate():
palette.py:126:def ramp(mat, hour, wet=0.0, haze=0.0, exposure=1.0, key=None, ambient=1.0):
palette.py:158:def sky_colors(hour, n=8):
palette.py:167:def sun_vector(az_deg, el_deg):
scene.py:153:def up(img, k):
raster.py:65:def _raster(S, persp, W, H, depth, tid, b0, b1):
raster.py:108:def raster(S, persp, W, H):
raster.py:120:def tri_normals(V):
raster.py:125:def clip_near(mesh, cam, near=0.3):
raster.py:173:def render(mesh, cam, sun_dir, shadow_res=None, shadow_extent=None):
jointfield.py:24:def _h(*a):
jointfield.py:35:def _warp(u, v, scale, amp, seed):
jointfield.py:47:def cells_granite(u, v, prm, seed):
jointfield.py:65:def cells_beds(u, v, prm, seed):
jointfield.py:83:def cells_voronoi(u, v, spacing, seed, jitter=0.45):
jointfield.py:95:def planarize(U, V, D, cid, rnd_tilt=0.12, recess=0.5, rounding=0.0, extra_tilt=None, extra_off=None):
jointfield.py:125:def granite_surface(U, V, D, prm, seed=0):
jointfield.py:141:def beds_surface(U, V, D, prm, seed=0):
jointfield.py:148:def column_surface(U, V, D, prm, seed=0):
jointfield.py:234:def massif_facets(U, V, D, size=30.0, tilt=0.45, seed=0, rounding=0.25):
jointfield.py:242:def hackly_surface(U, V, D, prm, seed=0):
jointfield.py:255:def knob_surface(X, Y, Z, prm, seed=0):
paint.py:21:def hash01(*a):
paint.py:31:def clump_noise(shape, sigma, seed):
paint.py:37:def masked_blur(a, m, sigma):
paint.py:43:def ambient_occlusion(g, radius_m=0.35, dirs=8):
paint.py:75:def gaps(g, min_jump_m=0.12, max_jump_m=1.5):
paint.py:88:def silhouette(g, sky_only=False):
paint.py:104:def face_boundaries(g):
paint.py:114:def shade(g, L, overcast=False, params=None, keymask=None):
paint.py:157:def quantize(v, lit, seed=0, grain=0.2, clump=1.2, seam=0.18):
paint.py:175:def key_family(v, m, lo, hi, q=(3, 97)):
paint.py:188:def found_edges(g, step, lit, L, cam, v):
paint.py:210:def paint_gray(g, L, overcast=False, seed=0, params=None, edges=True):
paint.py:224:def gray_rgb(s, sky=(200, 205, 215)):
paint.py:233:def massif_objects(g, min_frac=0.03, min_extent_m=12.0):
paint.py:249:def big_normals(g, mask, sigma_px):
paint.py:263:def design_masses(lit, mask, min_frac=0.025):
paint.py:278:def shade_designed(g, L, params=None, keymask=None, aspect_map=None, mass_override=None):
paint.py:386:def mass_marks(g, s, lit, mlit, mass, seed=0, keep=0.35, min_step_m=0.25, draw_dark=False):
paint.py:411:def despeckl

cd ~/work/artschool/posts; head -c 5000 master-copy/POST.md; echo; grep -il "master copy\|master-copy" */POST.md | head; head -c 3000 cobbles/TUTORIAL.md

# output:
# Copying Ferrari from his own files

*The master-copy student, writing the day after. The study took one evening: 1 hour 11 minutes, 237 steps, 93 images looked at. Every image below marked "nb" is one I looked at while working, numbered as in [the notebook](../../notebooks/master-copy/notebook.md). Several were overwritten on disk later, so the notebook is the only place they survive. The companion [tutorial](TUTORIAL.md) teaches the method without the story.*

The brief was the same small subject every arm got: a granite boulder in a grass bank at the edge of a mountain tarn, at dusk, in Mark Ferrari's manner. It had to be 320×180 pixels, at most 32 colours, and delivered as a program that paints it from parameters (size, light direction, distance, seed). My arm's teaching was the oldest one: learn from the master by copying him. Pick crops of his painting, write a program that reproduces each one, compare obsessively at 3× and 8×, pull out the rules, then generalise.

The short version of how it went: one lucky find made the copies much better than I had any right to expect, the copies got good, and the final painter was worse than the copies. I'll come back to why at the end.

![Ferrari's tarn, the reference](nb001.png)

*nb001. The reference as I first saw it. The whole study is an attempt to understand how this was made.*

---

## 1. What the file actually is

Before copying anything I wanted to see his real pixels. The PNG is 1268×938 with 172 colours. Run lengths along rows and columns were all even numbers:


d = (np.abs(np.diff(a, axis=1)).sum(2) > 0)          # where colour changes along a row
runs = [np.diff(np.flatnonzero(d[r])) for r in range(0, a.shape[0], 7)]
# Counter -> [(2, 35268), (4, 7235), (6, 3368), ...]   every run even: a 2x upscale



Then I tried the four possible 2×2 phases and upscaled each one back. Phase (0,0) reproduced the file with 0.0% error, so the native painting is 634×469. I saved that as `ref/native.png` and never looked at the big file again. A small helper, `zoom.py x y w h scale`, cropped native pixels and blew them up with nearest neighbour; I used it for every look.

![left bank at 4x](nb002.png)

*nb002. The lit bank at 4×. Vertical strokes of olive on ochre, stipple at the far edge. I could see the texture but not yet how it was made.*

![island at 5x](nb003.png)

*nb003. The small grass island. It looked like a soft gradient from green to orange, and it would turn out to be four exact classes.*

## 2. Reading, and the find that set up everything

I searched for what Ferrari has said about his process and read two interviews properly:

- **The EffectGames Q&A with Mark J. Ferrari.** He painted the colour-cycling scenes in Deluxe Paint II and depended on its palette management. He chooses dither by effect: "I use patterned dither because I want the effect to look smooth and precise," diffusion dither for chaotic things like collapsing surf, and no dither for solid bands. Each scene took "a

cd ~/work/pixelart-studies; head -60 tools/imgfetch.py; ls elements/rock/refs/*/ | head -80; du -sh elements/rock/refs

# output:
#!/usr/bin/env python3
"""imgfetch: fast, polite reference-image sourcing (Wikimedia Commons + Openverse).

  python3 imgfetch.py "cinder cone scoria slope" -n 12 -o refs/cinder
  python3 imgfetch.py -f queries.txt -n 8 -o refs            # one query per line; subdir per query
  options: --src commons,openverse   --width 1200   --sheet (writes <out>/_sheet.jpg)

Writes images plus <out>/sources.tsv (file, title, license, author, page url).
Cached: files already present are skipped. Handles 429/5xx with Retry-After backoff.
Parallel downloads (4 workers) with a shared per-host rate limiter (~3 req/s; Wikimedia's media
host is kept to 1 per 2 s and slows further on any 429, since it throttles hard per IP after bursts).
Openverse (mostly Flickr-hosted CC images) is the primary source because it stays fast.
"""
import argparse, concurrent.futures as cf, json, os, re, sys, threading, time, urllib.error, urllib.parse, urllib.request

UA = "pixelart-studies-imgfetch/1.0 (research; contact: local)"
_locks = {}; _last = {}; _glock = threading.Lock()

SLOW_HOSTS = {"upload.wikimedia.org": 2.0}   # Wikimedia media host: ~1 req/s or it 429s

def _throttle(host, min_gap=0.34):
    min_gap = SLOW_HOSTS.get(host, min_gap)
    with _glock:
        lk = _locks.setdefault(host, threading.Lock())
    with lk:
        dt = time.time() - _last.get(host, 0)
        if dt < min_gap:
            time.sleep(min_gap - dt)
        _last[host] = time.time()

def get(url, tries=4, timeout=25):
    host = urllib.parse.urlparse(url).netloc
    if host == "upload.wikimedia.org":
        tries = 2  # fail fast; other sources fill in
    for i in range(tries):
        _throttle(host)
        try:
            req = urllib.request.Request(url, headers={"User-Agent": UA})
            return urllib.request.urlopen(req, timeout=timeout).read()
        except urllib.error.HTTPError as e:
            if e.code == 429:
                with _glock:
                    SLOW_HOSTS[host] = min(8.0, SLOW_HOSTS.get(host, 0.34) * 2)  # adapt: this host wants us slower
            if e.code in (429, 500, 502, 503, 504) and i < tries - 1:
                wait = min(15, int(e.headers.get("Retry-After") or 0) or 2 * 2 ** i)
                time.sleep(wait); continue
            raise
        except (urllib.error.URLError, TimeoutError):
            if i < tries - 1:
                time.sleep(2 * (i + 1)); continue
            raise

def commons(query, n, width):
    q = urllib.parse.urlencode({"action": "query", "generator": "search", "gsrnamespace": 6, "gsrlimit": min(50, n * 3),
                                "gsrsearch": query + " filetype:bitmap", "prop": "imageinfo",
                                "iiprop": "url|size|extmetadata", "iiurlwidth": min((b for b in (330, 500, 960, 1280) if b >= width), default=1280), "format": "json"})
    d = json.loads(get("https://commons.wikimedia.org/w/api.php?" + q))
    out = []
    for p in sorted(d.get("query", {}).get("pages", {}).values(), key=lambda 

cd ~/work/pixelart-studies/books/genge; cat > refs/queries.txt <<'EOF'
Rough Tor Bodmin Moor granite tor
Ladram Bay sea stack
West Bay Bridport cliffs
Siccar Point unconformity
Yakima River Canyon basalt cliffs
Columbia River basalt flows cliff
Ruby Beach sea stack
Kalaloch sandstone cliff
Sucia Island sandstone honeycomb
Enchantments granite joints
Stair Hole Lulworth folds
EOF
nice -n 10 python3 ~/work/pixelart-studies/tools/imgfetch.py -f refs/queries.txt -n 8 -o refs --src openverse,commons --width 1200 --sheet 2>&1 | tail -15; du -sh refs

# output:
download failed: HTTP Error 429: Too many requests - please contact noc@wikimedia.org to discuss a less disruptive approach or instead use thumbnail images in sizes listed on https://w.wiki/GHai. (0c640b1)
  download failed: HTTP Error 429: Too many requests - please contact noc@wikimedia.org to discuss a less disruptive approach or instead use thumbnail images in sizes listed on https://w.wiki/GHai. (0c640b1)
  download failed: HTTP Error 429: Too many requests - please contact noc@wikimedia.org to discuss a less disruptive approach or instead use thumbnail images in sizes listed on https://w.wiki/GHai. (0c640b1)
  download failed: HTTP Error 429: Too many requests - please contact noc@wikimedia.org to discuss a less disruptive approach or instead use thumbnail images in sizes listed on https://w.wiki/GHai. (0c640b1)
'Rough Tor Bodmin Moor granite tor': 3 images -> refs/Rough_Tor_Bodmin_Moor_granite_tor
'Ladram Bay sea stack': 5 images -> refs/Ladram_Bay_sea_stack
'West Bay Bridport cliffs': 8 images -> refs/West_Bay_Bridport_cliffs
'Siccar Point unconformity': 5 images -> refs/Siccar_Point_unconformity
'Yakima River Canyon basalt cliffs': 8 images -> refs/Yakima_River_Canyon_basalt_cliffs
'Columbia River basalt flows cliff': 8 images -> refs/Columbia_River_basalt_flows_cliff
'Ruby Beach sea stack': 8 images -> refs/Ruby_Beach_sea_stack
'Kalaloch sandstone cliff': 3 images -> refs/Kalaloch_sandstone_cliff
'Sucia Island sandstone honeycomb': 6 images -> refs/Sucia_Island_sandstone_honeycomb
'Enchantments granite joints': 4 images -> refs/Enchantments_granite_joints
'Stair Hole Lulworth folds': 5 images -> refs/Stair_Hole_Lulworth_folds
 18M    refs

# Genge's model, in his own terms

Matthew J. Genge, *Geological Field Sketches and Illustrations: A Practical Guide* (Oxford, 2020).
Page numbers are the printed ones (PDF page = printed page + 8). Figures cited are the ones I
looked at in `plates/hi-*.png` (150 dpi renders of the PDF pages).

This file records what the book says, in its vocabulary, before any translation. The
translations, tests and verdicts are in `LESSONS.md`.

---

## 0. What a geological sketch is for

- **A sketch records information; its aesthetic qualities are not important** (p. 18). "It is the
  ability of a geological sketch to record information that defines its quality."
- **Why draw rather than photograph** (pp. 2–4, Fig. 1.1 Ladram Bay): drawing forces observation;
  photographs are one perspective, flatten colour differences (the "pink rhyolite" is a grey-pink,
  and in a photo "all the same dull grey", p. 4), and **deep sun shadows hide features the eye traces
  from light into shade** (p. 4). Fig. 1.1 is the thesis in one plate: in the photo the Ladram stack
  is a dark red lump; in the sketch the same stack is a stack of beds, ledges and lenses.
- **Illustration vs art** (p. 28): an illustration "is a simplification of reality that is largely
  true in form and scale to the subject, but may emphasize certain features of interest". "Every
  line in an illustration has a reason and relates to the observed subject."
- Realism "is suited to science since it represents the subject without artefact and stylistic
  convention" (p. 16, on Waller's ammonites).

## 1. The three rules (§2.7, pp. 30–31)

1. **Look first, then draw.** Decide what the sketch is *for* (the subject: a fault, a fold, a
   contact) before the first mark. Often there is no reason to draw at all: "an outcrop with planar
   beds of one rock type … does not warrant a field sketch" (p. 30).
2. **Don't draw too much or too little.** Choose the field of view as the subject plus enough context
   to read it (for a fault: the beds on both sides, to show displacement).
3. **Don't draw too small** — half a page at least (p. 31).

## 2. The stages of drawing (§2.8, pp. 31–33; worked in every chapter)

The procedure is the same for faults, folds, landscapes, hand specimens and fossils:

| stage | what goes in | book's words |
|---|---|---|
| 0. proportions | width : height of the view, measured (pencil at arm's length); a box | "It is most common to accidently exaggerate the height" (p. 29) |
| 1. **blocking-in** | the outline of the outcrop or the horizon, the subject, and the features that control the positions of everything else, **as straight-line segments**, placed with a **quadrant grid** | "The first features to be drawn are those that control the positions of other features" (p. 56) |
| 2. amend | protrusions and indents added to the straight lines, "care … not to exaggerate" | p. 31 |
| 3. essential detail | bedding traces, minor faults, folds that enhance the context of the subject | p. 33 |
| 4. ancillary detail | clasts, joints, laminations, vegetation, "when time is available" | p. 33 |
| 5. labels, scale, looking-towards bearing; ink; colour; shading | | pp. 33–38 |

**What controls what** (the order of blocking-in) is stated per subject:

- faults before beds (Moab, p. 51); fold **axial traces** (imagined, interpolated) before beds
  (pp. 58–59); in superimposed folds the **last** generation first because its axial traces are
  straight (p. 256);
- **horizons** first in landscapes, then gullies, then geological contacts (pp. 105–106);
- **bedding** first in sedimentary outcrops (p. 141, "Bedding is the most crucial feature");
- the **outline** of the specimen first for hand specimens (p. 197);
- for many similar objects (bivalves, clasts), outline the *area* they occupy first, then start from
  one characteristic object and work outwards (p. 177).

## 3. Line (§2.5 p. 28, Fig. 2.5; §8.5 p. 139; §6.5.2 p. 101)

- **Clean, precise, simple lines of constant width.** No scratchy back-and-forth lines (Fig. 8.13).
- **Line width is emphasis** (Fig. 2.5, Rough Tor): "thick line emphasizes horizon", "thick line
  emphasizes foreground", "thin lines are clean". Elsewhere the thick line goes on: the skyline
  (Fig. 6.7a "Thicker lines used to emphasise skyline"), **protruding competent beds "that have
  shadowed overhangs at their base"** (p. 106, Fig. 7.2b), topography (Fig. 7.4a), the subject
  itself (the unconformity, p. 154), the outlines of lava tongues (p. 127).
- **Discontinuous ("impersistent") lines** for subtle features — laminations, cleavage, faint
  bedding — "should have stroke lengths that reflect the scale of the objects observed" (p. 28).
  Figs 3.2a, 4.3, 6.7, 9.6a ("Draw indistinct bedding traces as discontinuous lines").
- **Regular dashes are reserved for interpretation** (extrapolated boundaries), never for subtle
  observed features (p. 28; p. 114; p. 176).
- **Line shading (hatching) is a mistake** (Fig. 6.9): "it introduces lines that do not exist in
  reality … could be confused with features such as fabrics and lineations". Stipple likewise
  (p. 37).
- **Symbolic patterns are a mistake** (Fig. 9.15, p. 166): dots for sandstone, circles for
  conglomerate, brickwork for limestone "introduce features into diagrams that do not exist".

## 4. Draw what is seen, not what is expected

The book's most repeated warning, in six costumes:

- **Planar feature on a curved surface**: "The feature is planar, thus its line of intersection is
  erroneously drawn straighter than it should be" (p. 88).
- **Schematic bedding** (Fig. 4.7): folds right, but "beds are schematic and almost equal in
  width" (p. 65).
- **Constant bed thickness** "is a case of drawing an expectation rather than an observation"
  (Fig. 5.12, p. 86).
- **Schematic cleavage** (Fig. 5.13): "rarely present as equally spaced, parallel fractures … usually
  better developed in some places than others" (p. 86). Same for cleavage in a hand specimen: "the
  spacing and length of the lines is variable, since these structures are not regular" (p. 200).
- **Exaggerated perspective** (Fig. 6.8): beds drawn "as rectangular blocks with unnatural looking
  geometry" (p. 100). **Vertical exaggeration** ~1.5× (Fig. 3.8, p. 54) makes faults too steep.
- **Idealized crystals and fossils** (Fig. 11.8; p. 184).

And the opposite failure: **oversimplification** (Fig. 3.9) and **insufficient detail** (Fig. 7.8,
landscapes with no gullies, tree lines or scree; Fig. 11.9, hand specimens without texture).

## 5. What the rock is actually like (the statistics the sketches carry)

- **Beds vary in thickness vertically and laterally**, pinch, lens and disappear (pp. 141–144; Figs
  3.2, 6.7, 7.2). "The cemented layers become more closely spaced and thinner towards the top of the
  cliff" (p. 98). Bed thickness trends are data: thickening-up, cyclicity (p. 144).
- **Competent beds protrude, soft beds recess** (greywacke ribs vs slate grooves, p. 150; cemented
  sandstone ledges, p. 98; "protruding competent beds drawn with heavier lines", Fig. 7.2b).
- **Joints**: usually perpendicular to bedding (p. 61); add a few for realism but "it is important
  not to obscure the important features by adding too many" (p. 61); **over-emphasized jointing is a
  named mistake** (Fig. 4.6, p. 65). "Selective addition of joints for realism" (Fig. 8.8b).
- **Granite tors** (Rough Tor, Figs 2.5 and 2.9): horizontal sheeting joints make stacked,
  rounded "pillows"; vertical joints cut them into blocks; the outline is a stepped, lumpy skyline.
- **Lava flows** (§8.2, Fig. 8.4): flow top (vesicular/rubbly/breccia), massive interior "can …
  contain columnar cooling joints", basal breccia; flows are stacked. On a crater wall (Vesuvius,
  Fig. 8.8): "boundaries between lava flows. Not all are continuous"; "selective addition of joints
  for realism"; "slope lines on scree slope"; "larger boulders drawn semi-schematically".
- **Metamorphic banding** thickens in hinges, thins on limbs, boudinages (pp. 75, 80, 85).

## 6. Three dimensions (Ch. 6)

- **Draw the view as a 2-D image**, not a constructed perspective (p. 87, p. 101). Vanishing
  points are not used in field sketches (p. 89).
- **The key 3-D lines are the edges between bedding planes and joints** — "the upper and lower
  intersections of the bedding surface with joints" (p. 91); "3D beds can be shown by drawing upper
  and lower limits" (Fig. 9.6a).
- **Bluffs make local horizons** (§6.4, Figs 6.6–6.7): viewed along a cliff, each bluff is a
  discontinuity at which bedding "appears to suddenly change in position and in width", resembling
  a fault. The skyline becomes a line "with several steps, each segment located over a bluff", and
  local horizons drop down from those steps. Draw each segment as observed.
- **Perspective thins beds with distance** (p. 87); in the distance "only a proportion of the layers
  can be drawn" (p. 98).
- **Furrows on the face** (weathering, no geological significance) "allow the geometry of the cliff
  face to be represented more accurately" (p. 98, Fig. 6.7a "Furrows on cliff show the surface").
- Scale: put the scale bar at the subject's distance, or give a spot height (pp. 92, 98).

## 7. Landscapes (Ch. 7)

- **Horizons first.** A view is divided by **multiple horizons into regions of different scale and
  detail** (§7.1, Fig. 7.1). Simplify each as straight segments, then add indents.
- **Field of view**: small enough that the head need not turn (p. 105); widening it lowers the
  resolution of what can be drawn.
- **Gullies and ravines as lines** give the slope volume (p. 106, Fig. 7.2b "Minor gulleys provide
  volume to crater surface").
- **Bedding traces on topography**: "Pay attention to topographic rises and falls on boundaries"
  (Fig. 7.2a); thick lines on protruding competent beds; discontinuous where obscured by scree or
  vegetation; "faint bedding trace on landscape" as short parallel discontinuous lines (p. 106).
  In the Aliaga plate (Fig. 7.4) the bedding wraps every rib and gully of the mountainside.
- **Scree** "semi-schematically using discontinuous lines" (p. 111); **tree line and forest** as
  "quickly drawn irregular zigzag lines" (pp. 111, 115); **areas of exposure are outlined first**,
  then filled with geology (§7.3).
- **Mistakes**: insufficient topographic detail (Fig. 7.8); many schematic tree symbols (Fig. 7.9)
  "take more time than drawing the trees more realistically".

## 8. Shading (§2.9.3 pp. 36–38, Fig. 2.9; also pp. 91–92, 127, 131, 154, 259–261)

- **Do not reproduce cast shadows**: directional light "generates deep directional shadows that have
  sharp boundaries … These shadows obscure detail and should not be included in a sketch" (p. 37).
- **Reproduce ambient occlusion** (diffuse-light shadow): "darkest in areas of a surface that are
  exposed to the least light sources such as within cracks and indents or on the undersides of
  ledges and protrusions" (p. 37).
- Fig. 2.9's labels are the rules: **"inset areas shadowed"**, **"shadow indicates top of
  overhang"**, **"deep fractures are dark"**, **"shadows can emphasize sudden change in orientation
  of surface"**, **"darkest shadow in most inset areas"**.
- **Two layers** (grey for subtle, black for deep) "can often increase realism, but care should be
  taken not to overwhelm geological features" (p. 37).
- Gradational, by pressure; **colour, not parallel ink lines** (p. 92).
- Shading "is somewhat schematic"; "often it is best not to try and reproduce the complexity of
  shadow on an outcrop, which is a function of light angle" (p. 91).
- In block diagrams, "adding darker colours to areas where significant changes in slope occur, such
  as close to the rivers or at the base of the cliffs generates natural looking shadows" (p. 259).

## 9. Colour (§2.9.2 pp. 34–36; pp. 48, 131, 154)

- **Flat, homogeneous base colour** per unit, long light strokes, no visible streaks (p. 35).
- **Exaggerate the differences** between units: outcrop colours are subtle; "often it is better to
  emphasize colour" (p. 35); "colours are subtle in outcrop but can be exaggerated in the sketch to
  aid distinction" (p. 154). Colour lets many labels go (p. 48).
- **Avoid blending two pencils**: "it can have poor results and turn a sketch into a sickly coloured
  art joke" (p. 35).
- **Texture as a second layer over the base**: swirls for clasts too small to ink, "linear or curved
  streaks" along planar fabrics and layering (p. 35, Fig. 2.8); several passes with increasing
  pressure.
- Lustre and transparency are hard; record the base colour and "where the crystal is lighter and
  darker" (pp. 200, 205).

## 10. Many small things (clasts, crystals, cleavage, shells)

- **Draw enough to illustrate the texture, not all of them** (p. 46; p. 201: "Although there are
  hundreds of groundmass orthoclase crystals present only ~30 were drawn, which seems to be enough to
  demonstrate the texture").
- **Largest and most significant in full; the rest implied** (p. 201).
- **Exaggerate size slightly** so they are visible (p. 137, p. 157), but **shapes realistic**
  (p. 157, p. 167).
- **Self-similar series by counting**: count ribs against a reference line (five to halfway…)
  so spacing is right (p. 184).
- **Peripheral (inset) diagrams** carry detail at another scale, each with its own scale bar
  (pp. 85, 137–139, 159–160).

## 11. Schematic and block diagrams (Ch. 15)

- **Orthographic** block diagrams: "opposing sides of the block are parallel and of the same
  length"; the brain reads depth from them (p. 255). (Our 30° iso is this projection.)
- **Landscape surface "with the minimum number of lines"** (p. 259): ridge-lines, upper margins of
  topographic highs, **the lower boundary of cliff faces**, the upper margins of a gorge's cliffs;
  **short high-angle lines on cliffs** "reinforce the impression of slope"; **sub-horizontal lines
  on cliffs for bedding** (Fig. 15.4b); rivers drawn sub-parallel to the block edges to read as
  flat.
- Subsurface structure controls the surface features (p. 259). Vertical exaggeration is normal in
  *schematic* block diagrams (p. 257) — but not in sketches.
- **Digital painting** (§16.4.4, Fig. 16.7): base colours blocked in; a speckle layer smeared
  *along the bedding*; a warp for minor folds; **highlights and lowlights along the top and bottom of
  beds**; thin dark lines at fractures and faults.

## 12. Maps and sections (Chs 13–14): boundaries meet topography

- "On flat ground a boundary will be oriented parallel to the strike of bedding, whilst on a slope the
  boundary will run down slope of the strike line in the direction of dip … **Only a horizontal
  boundary will run parallel to contours. A vertical boundary, in contrast, trends straight across the
  landscape**" (p. 229). (The rule of Vs.)
- Boundaries are drawn by confidence: observed solid, constrained long-dash, inferred short-dash
  (p. 229).
- A structure is not aligned to the topography line (Fig. 14.6, p. 252).

## 13. Vocabulary (the book's terms I use in LESSONS.md)

blocking-in; quadrant grid; guidelines; key/essential/ancillary detail; peripheral diagram; local
horizon; bluff; bedding trace; impersistent/discontinuous line; competent/incompetent bed; ledge;
axial trace; hinge; limb; interlimb angle; joint; sheeting; cleavage; foliation; flow top; colonnade
(not the book's word — it says "columnar cooling joints", p. 124); scree; ambient occlusion;
inset area; overhang; looking-towards direction; spot height.

cd ~/work/pixelart-studies/books/genge/plates; P=~/work/pixelart-studies/library/purchased/"Geological Field Sketches and Illustrations _ A Practical -- Matthew J_ Genge.pdf"; for p in 11 37 45 75 107 118 138 161 268; do nice -n 10 pdftoppm -r 300 -png -f $p -l $p "$P" x300; done; ls x300*; python3 -c "
from PIL import Image; im=Image.open('x300-045.png'); print(im.size)"

# output:
x300-011.png
x300-037.png
x300-045.png
x300-075.png
x300-107.png
x300-118.png
x300-138.png
x300-161.png
x300-268.png
(1630, 2552)

cd ~/work/pixelart-studies/books/genge; cat > src/crops.py <<'EOF'
"""Figure crops from the 300 dpi page renders (boxes in 300 dpi px: x0,y0,x1,y1)."""
from PIL import Image
CROPS = {
 'fig2_5_roughtor_lines': ('x300-037.png', (440, 240, 1200, 760)),
 'fig2_9a_roughtor_ao':   ('x300-045.png', (230, 1600, 830, 2020)),
 'fig2_9b_roughtor_ao':   ('x300-045.png', (850, 1600, 1420, 2020)),
 'fig1_1_ladram_photo':   ('x300-011.png', (200, 240, 1400, 930)),
 'fig1_1_ladram_sketch':  ('x300-011.png', (200, 930, 1400, 1660)),
 'fig4_6_joints':         ('x300-075.png', (190, 230, 1420, 990)),
 'fig4_7_schematic':      ('x300-075.png', (190, 1110, 1420, 1820)),
 'fig6_7a_westbay':       ('x300-107.png', (500, 280, 1420, 1060)),
 'fig6_7b_westbay_col':   ('x300-107.png', (500, 1150, 1420, 1960)),
 'fig7_4a_aliaga':        ('x300-118.png', (300, 820, 1420, 1300)),
 'fig7_4b_aliaga_col':    ('x300-118.png', (300, 1530, 1420, 2170)),
 'fig8_8b_vesuvius':      ('x300-138.png', (320, 900, 1400, 1500)),
 'fig8_8c_vesuvius_col':  ('x300-138.png', (320, 1540, 1400, 2170)),
 'fig9_6b_siccar':        ('x300-161.png', (400, 960, 1300, 1500)),
 'fig9_6c_siccar_col':    ('x300-161.png', (250, 1560, 1300, 2150)),
 'fig15_4b_block':        ('x300-268.png', (320, 900, 1400, 1520)),
 'fig15_4c_block_col':    ('x300-268.png', (320, 1570, 1440, 2210)),
}
if __name__ == '__main__':
    import os
    os.makedirs('plates/crops', exist_ok=True)
    for k, (f, box) in CROPS.items():
        Image.open('plates/' + f).convert('RGB').crop(box).save(f'plates/crops/{k}.png')
        print(k)
EOF
python3 src/crops.py >/dev/null; python3 -c "
from PIL import Image; import glob
fs=sorted(glob.glob('plates/crops/*.png')); ims=[Image.open(f) for f in fs]
W=4; w=400
rows=[]
import math
sheet=Image.new('RGB',(W*w, math.ceil(len(ims)/W)*300),'white')
for i,im in enumerate(ims):
    im.thumbnail((w,300)); sheet.paste(im,((i%W)*w,(i//W)*300))
sheet.save('img/crops_check.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; cat > src/util.py <<'EOF'
"""Small shared helpers: loading, area-downsampling, nearest upscaling, labelled sheets."""
import numpy as np
from PIL import Image, ImageDraw, ImageFont

def load(p, gray=False):
    im = Image.open(p).convert('L' if gray else 'RGB')
    return np.asarray(im).astype(np.float32) / 255.0

def area_down(a, w=None, h=None, f=None):
    """Box (area) downsample to width w (keeps aspect) — the honest 'what a pixel sees'."""
    H, W = a.shape[:2]
    if f is None:
        f = W / float(w)
    h2, w2 = int(round(H / f)), int(round(W / f))
    im = Image.fromarray((np.clip(a, 0, 1) * 255).astype(np.uint8))
    return np.asarray(im.resize((w2, h2), Image.BOX)).astype(np.float32) / 255.0

def up(a, k):
    a = np.asarray(a)
    return np.repeat(np.repeat(a, k, 0), k, 1)

def to8(a):
    a = np.asarray(a)
    if a.dtype != np.uint8:
        a = (np.clip(a, 0, 1) * 255 + 0.5).astype(np.uint8)
    if a.ndim == 2:
        a = np.stack([a] * 3, -1)
    return a

def font(sz=14):
    for f in ['/System/Library/Fonts/Supplemental/Arial.ttf', '/System/Library/Fonts/Helvetica.ttc']:
        try:
            return ImageFont.truetype(f, sz)
        except Exception:
            pass
    return ImageFont.load_default()

def sheet(rows, title=None, pad=10, bg=(245, 243, 238), label_h=20, fs=14, title_fs=18):
    """rows: list of rows; each row a list of (label, uint8/float image). Images pasted as given."""
    rows = [[(l, to8(im)) for l, im in r] for r in rows]
    colw = max(len(r) for r in rows)
    ws = [pad + sum(im.shape[1] + pad for _, im in r) for r in rows]
    hs = [max(im.shape[0] for _, im in r) + label_h + pad for r in rows]
    th = (title_fs + 12) if title else 0
    W, H = max(ws), sum(hs) + pad + th
    out = Image.new('RGB', (W, H), bg)
    d = ImageDraw.Draw(out)
    if title:
        d.text((pad, 6), title, fill=(20, 20, 20), font=font(title_fs))
    y = pad + th
    f = font(fs)
    for r, h in zip(rows, hs):
        x = pad
        for l, im in r:
            d.text((x, y), l, fill=(30, 30, 30), font=f)
            out.paste(Image.fromarray(im), (x, y + label_h))
            x += im.shape[1] + pad
        y += h
    return out
EOF
echo ok

# output:
ok

cd ~/work/pixelart-studies/books/genge; ls refs/Rough_Tor_Bodmin_Moor_granite_tor/; python3 -c "
from PIL import Image; import glob
for f in glob.glob('refs/Rough_Tor_Bodmin_Moor_granite_tor/*.jpg'): print(f, Image.open(f).size)"

# output:
_sheet.jpg
commons_Roughtor_Bodmin_Moor_2_jpg_76d7aed8.jpg
openverse_3b5b1a69_4bf4_4ef0_864d_db763bf0eef8_7b7ffd6e.jpg
openverse_Bodmin_Moor_Z17305_85dc9a65.jpg
sources.tsv
refs/Rough_Tor_Bodmin_Moor_granite_tor/commons_Roughtor_Bodmin_Moor_2_jpg_76d7aed8.jpg (1280, 853)
refs/Rough_Tor_Bodmin_Moor_granite_tor/_sheet.jpg (1440, 376)
refs/Rough_Tor_Bodmin_Moor_granite_tor/openverse_3b5b1a69_4bf4_4ef0_864d_db763bf0eef8_7b7ffd6e.jpg (1023, 678)
refs/Rough_Tor_Bodmin_Moor_granite_tor/openverse_Bodmin_Moor_Z17305_85dc9a65.jpg (1024, 768)

cd ~/work/pixelart-studies/books/genge; mkdir -p stages/e1; cat > src/e1_mastercopy_downsample.py <<'EOF'
"""E1a. Master copy by honest downsampling: what does a pixel grid do to Genge's Rough Tor drawing?
The tor in Fig 2.9 spans ~12 m (Rough Tor stacks are 4-8 m high); photo stack ~9 m across."""
import numpy as np, sys
sys.path.insert(0, 'src')
from util import load, area_down, up, sheet

SCALES = [('close 0.1 m/px', 0.10), ('landscape 0.35 m/px', 0.35), ('region 1 m/px', 1.0)]
fig = load('plates/crops/fig2_9b_roughtor_ao.png', gray=True)[60:, :]    # drop the '(b)' label strip
lin = load('plates/crops/fig2_5_roughtor_lines.png', gray=True)[90:470, 0:640]
ph = load('refs/Rough_Tor_Bodmin_Moor_granite_tor/commons_Roughtor_Bodmin_Moor_2_jpg_76d7aed8.jpg')[60:853, 20:1280]
rows = []
for name, src, span in [('Fig 2.5 lines', lin, 12.0), ('Fig 2.9b AO', fig, 12.0), ('photo', ph, 9.0)]:
    r = []
    for lab, mpp in SCALES:
        w = max(8, int(round(span / mpp)))
        d = area_down(src, w)
        k = max(1, int(round(360 / d.shape[1])))
        r.append((f'{name} | {lab} ({d.shape[1]}x{d.shape[0]}) x{k}', up(d, k)))
    rows.append(r)
sheet(rows, title='E1a  Rough Tor: Genge drawings vs photo, box-downsampled to our scales').save('stages/e1/e1a_downsample.png')
EOF
python3 src/e1_mastercopy_downsample.py && echo done

# output:
done

"""engine.py — a heightfield in any frame, ray-marched to a G-buffer.

Every subject in this study is a heightfield h(i, j) on a regular grid (spacing ds metres) laid in a
frame: origin O, in-plane axes E1 (along j) and E2 (along i), and normal N (the direction h grows).

  - a cliff face or a tor: E1 = east, E2 = up, N = toward the viewer; h = how far the rock stands out
  - terrain: E1 = east, E2 = north, N = up; h = elevation
Shading terms (sky occlusion, sun horizon, normals) are computed once in texture space, where they are
cheap and exact for a heightfield; the camera then ray-marches the heightfield and the G-buffer holds,
per screen pixel, the texel it hit. Painting (steps, lines, dither) happens in screen space, so every
decision is made at the pixel it lands on.
"""
import numpy as np
from numba import njit, prange


class Frame:
    def __init__(self, O, E1, E2, N):
        self.O, self.E1, self.E2, self.N = (np.asarray(v, np.float64) for v in (O, E1, E2, N))

    def to_local(self, P):
        d = P - self.O
        return np.stack([d @ self.E1, d @ self.E2, d @ self.N], -1)

    def to_world(self, L):
        return self.O + L[..., :1] * self.E1 + L[..., 1:2] * self.E2 + L[..., 2:3] * self.N


def normals(h, ds):
    gy, gx = np.gradient(h, ds)           # gy along i (E2), gx along j (E1)
    n = np.stack([-gx, -gy, np.ones_like(h)], -1)
    return n / np.linalg.norm(n, axis=-1, keepdims=True)   # local frame (E1, E2, N)


@njit(parallel=True, cache=True)
def _horizon(h, ds, dirs, radius, out):
    H, W = h.shape
    nd = dirs.shape[0]
    steps = int(radius / ds)
    for i in prange(H):
        for j in range(W):
            acc = 0.0
            for k in range(nd):
                dx, dy = dirs[k, 0], dirs[k, 1]
                best = 0.0
                s = 1
                while s <= steps:
                    x = j + dx * s
                    y = i + dy * s
                    xi, yi = int(x + 0.5), int(y + 0.5)
                    if xi < 0 or yi < 0 or xi >= W or yi >= H:
                        break
                    dist = s * ds
                    t = (h[yi, xi] - h[i, j]) / dist
                    if t > best:
                        best = t
                    s += 1 if s < 8 else 2
                acc += best / np.sqrt(1.0 + best * best)   # sin(horizon elevation)
            out[i, j] = acc / nd


def sky_occlusion(h, ds, radius=2.0, nd=12):
    """Mean sine of the horizon angle over nd azimuths (0 = open sky, ->1 = deep slot)."""
    a = np.arange(nd) * 2 * np.pi / nd
    dirs = np.stack([np.cos(a), np.sin(a)], -1)
    out = np.zeros_like(h, dtype=np.float64)
    _horizon(h.astype(np.float64), ds, dirs, radius, out)
    return out


@njit(parallel=True, cache=True)
def _sun_shadow(h, ds, lx, ly, lz, maxd, out):
    H, W = h.shape
    horiz = np.sqrt(lx * lx + ly * ly)
    if horiz < 1e-6:
        return
    dx, dy = lx / horiz, ly / horiz
    slope = lz / horiz
    steps = int(maxd / ds)
    for i in prange(H):
        for j in range(W):
            s = 1
            while s <= steps:
                x = j + dx * s
                y = i + dy * s
                xi, yi = int(x + 0.5), int(y + 0.5)
                if xi < 0 or yi < 0 or xi >= W or yi >= H:
                    break
                if h[yi, xi] > h[i, j] + slope * s * ds + 0.02:
                    out[i, j] = 1.0
                    break
                s += 1


def sun_shadow(h, ds, L_local, maxd=60.0):
    out = np.zeros_like(h, dtype=np.float64)
    if L_local[2] <= 0:
        out[:] = 1.0
        return out
    _sun_shadow(h.astype(np.float64), ds, float(L_local[0]), float(L_local[1]), float(L_local[2]), maxd, out)
    return out


@njit(parallel=True, cache=True)
def _march(h, ds, O, E1, E2, N, eye, dirs, ortho, fwd, tmax, hi, hj, hd):
    H, W = h.shape
    R = dirs.shape[0]
    hmax = h.max()
    hmin = h.min()
    for r in prange(R):
        # ray origin + direction in local coordinates
        if ortho:
            p = eye[r]
            d = fwd
        else:
            p = eye[0]
            d = dirs[r]
        px = (p[0] - O[0]) * E1[0] + (p[1] - O[1]) * E1[1] + (p[2] - O[2]) * E1[2]
        py = (p[0] - O[0]) * E2[0] + (p[1] - O[1]) * E2[1] + (p[2] - O[2]) * E2[2]
        pz = (p[0] - O[0]) * N[0] + (p[1] - O[1]) * N[1] + (p[2] - O[2]) * N[2]
        dx = d[0] * E1[0] + d[1] * E1[1] + d[2] * E1[2]
        dy = d[0] * E2[0] + d[1] * E2[1] + d[2] * E2[2]
        dz = d[0] * N[0] + d[1] * N[1] + d[2] * N[2]
        hi[r] = -1.0
        # skip to the slab hmin..hmax
        t = 0.0
        if pz > hmax and dz < 0:
            t = (pz - hmax) / (-dz) - ds
            if t < 0:
                t = 0.0
        step = ds * 0.5
        prev_above = True
        while t < tmax:
            x = px + dx * t
            y = py + dy * t
            z = pz + dz * t
            if dz < 0 and z < hmin - 1.0:
                break
            fj = x / ds
            fi = y / ds
            if fj < 0 or fi < 0 or fj > W - 1.001 or fi > H - 1.001:
                if (dx > 0 and fj > W) or (dx < 0 and fj < 0) or (dy > 0 and fi > H) or (dy < 0 and fi < 0):
                    if z < hmax:
                        break
                t += step * 4
                continue
            j0, i0 = int(fj), int(fi)
            aj, ai = fj - j0, fi - i0
            hh = (h[i0, j0] * (1 - aj) * (1 - ai) + h[i0, j0 + 1] * aj * (1 - ai)
                  + h[i0 + 1, j0] * (1 - aj) * ai + h[i0 + 1, j0 + 1] * aj * ai)
            if z <= hh:
                # bisection refine
                a, b = t - step, t
                for _ in range(10):
                    m = 0.5 * (a + b)
                    x2 = px + dx * m
                    y2 = py + dy * m
                    z2 = pz + dz * m
                    fj2, fi2 = x2 / ds, y2 / ds
                    j2 = min(max(int(fj2), 0), W - 2)
                    i2 = min(max(int(fi2), 0), H - 2)
                    aj2, ai2 = fj2 - j2, fi2 - i2
                    h2 = (h[i2, j2] * (1 - aj2) * (1 - ai2) + h[i2, j2 + 1] * aj2 * (1 - ai2)
                          + h[i2 + 1, j2] * (1 - aj2) * ai2 + h[i2 + 1, j2 + 1] * aj2 * ai2)
                    if z2 <= h2:
                        b = m
                    else:
                        a = m
                hi[r] = py / ds + dy * b / ds
                hj[r] = px / ds + dx * b / ds
                hd[r] = b
                break
            t += step


class Camera:
    """Perspective (fov_deg) or orthographic (ortho_w metres wide) camera."""
    def __init__(self, eye, target, W, H, fov_deg=40.0, ortho_w=None, up=(0, 0, 1)):
        self.eye = np.asarray(eye, float)
        self.W, self.H = W, H
        f = np.asarray(target, float) - self.eye
        self.fwd = f / np.linalg.norm(f)
        r = np.cross(self.fwd, up)
        if np.linalg.norm(r) < 1e-6:
            r = np.cross(self.fwd, (0, 1, 0))
        self.right = r / np.linalg.norm(r)
        self.upv = np.cross(self.right, self.fwd)
        self.fov, self.ortho_w = fov_deg, ortho_w

    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, ys)
        if self.ortho_w:
            hw = self.ortho_w / 2
            hh = hw * H / W
            eye = self.eye + X[..., None] * hw * self.right + Y[..., None] * hh * self.upv
            return eye.reshape(-1, 3), np.tile(self.fwd, (W * H, 1)), True
        t = np.tan(np.radians(self.fov) / 2)
        d = self.fwd + X[..., None] * t * self.right + Y[..., None] * t * H / W * self.upv
        d /= np.linalg.norm(d, axis=-1, keepdims=True)
        return np.tile(self.eye, (1, 1)), d.reshape(-1, 3), False


def march(h, ds, frame, cam, tmax=5000.0):
    eye, dirs, ortho = cam.rays()
    R = cam.W * cam.H
    hi = np.full(R, -1.0)
    hj = np.zeros(R)
    hd = np.zeros(R)
    _march(h.astype(np.float64), float(ds), frame.O, frame.E1, frame.E2, frame.N,
           np.ascontiguousarray(eye, dtype=np.float64), np.ascontiguousarray(dirs, dtype=np.float64),
           ortho, cam.fwd.astype(np.float64), tmax, hi, hj, hd)
    hit = hi >= 0
    return (hi.reshape(cam.H, cam.W), hj.reshape(cam.H, cam.W), hd.reshape(cam.H, cam.W),
            hit.reshape(cam.H, cam.W))


class GBuf:
    """Per-pixel lookups of texture-space fields at the texel each ray hit."""
    def __init__(self, hi, hj, depth, hit, ds):
        self.hit, self.depth, self.ds = hit, depth, ds
        self.I = np.clip(np.round(hi).astype(int), 0, None)
        self.J = np.clip(np.round(hj).astype(int), 0, None)
        self.fi, self.fj = hi, hj

    def take(self, tex, fill=0):
        H, W = tex.shape[:2]
        I = np.clip(self.I, 0, H - 1)
        J = np.clip(self.J, 0, W - 1)
        out = tex[I, J]
        if out.ndim == 2:
            out = np.where(self.hit, out, fill)
        else:
            out = np.where(self.hit[..., None], out, fill)
        return out


def combine(bufs):
    """Nearest-hit merge of several (GBuf, tag) — e.g. a cliff face plus the beach in front."""
    best = None
    for k, g in enumerate(bufs):
        d = np.where(g.hit, g.depth, np.inf)
        if best is None:
            best, owner = d, np.where(g.hit, k, -1)
        else:
            m = d < best
            best = np.where(m, d, best)
            owner = np.where(m, k, owner)
    return owner, best

"""palette.py — ramps of 8 steps per material: 0-3 shadow family (sky-lit, cooler), 4-7 lit family.

Targets are CIELAB (L*, C*, hue°) for the lit face (step 6) and the shadow face (step 2), measured
from photos by the rock study (elements/rock/LESSONS.md 1.2-1.6) or by me (noted). Everything else is
interpolated in Lab, so a palette swap relights an index canvas (studio convention: palette as light).
Genge's colour rule (p. 35, p. 154): flat base per unit, and *exaggerate* the difference between units
— `unit_ramps(base, n, spread)` makes n sibling ramps whose hue/value differences are pushed apart.
"""
import numpy as np


def lab2rgb(L, a, b):
    fy = (L + 16) / 116
    fx = fy + a / 500
    fz = fy - b / 200
    def finv(t):
        return np.where(t > 6 / 29, t ** 3, 3 * (6 / 29) ** 2 * (t - 4 / 29))
    X, Y, Z = 0.95047 * finv(fx), finv(fy), 1.08883 * finv(fz)
    r = 3.2406 * X - 1.5372 * Y - 0.4986 * Z
    g = -0.9689 * X + 1.8758 * Y + 0.0415 * Z
    bb = 0.0557 * X - 0.2040 * Y + 1.0570 * Z
    rgb = np.stack([r, g, bb], -1)
    rgb = np.where(rgb <= 0.0031308, 12.92 * rgb, 1.055 * np.clip(rgb, 0, None) ** (1 / 2.4) - 0.055)
    return np.clip(rgb, 0, 1)


def lch(L, C, h):
    return np.array([L, C * np.cos(np.radians(h)), C * np.sin(np.radians(h))])


def rgb2lab(rgb):
    rgb = np.asarray(rgb, float)
    lin = np.where(rgb <= 0.04045, rgb / 12.92, ((rgb + 0.055) / 1.055) ** 2.4)
    X = (0.4124 * lin[..., 0] + 0.3576 * lin[..., 1] + 0.1805 * lin[..., 2]) / 0.95047
    Y = 0.2126 * lin[..., 0] + 0.7152 * lin[..., 1] + 0.0722 * lin[..., 2]
    Z = (0.0193 * lin[..., 0] + 0.1192 * lin[..., 1] + 0.9505 * lin[..., 2]) / 1.08883
    f = lambda t: np.where(t > (6 / 29) ** 3, np.cbrt(t), t / (3 * (6 / 29) ** 2) + 4 / 29)
    fx, fy, fz = f(X), f(Y), f(Z)
    return np.stack([116 * fy - 16, 500 * (fx - fy), 200 * (fy - fz)], -1)


# lit (step 6) and shadow (step 2) targets, L*, C*, hue
MATS = {
    # Rough Tor / fresh-ish Cornish granite, lichen grey (my measure of the Commons photo: lit L 68, shadow L 30)
    'tor_granite':   dict(lit=(70, 6, 85), sh=(30, 4, 250)),
    # Enchantments granite: old wall lit 46 / shadow 22; fresh talus 67 (rock study 1.2)
    'granite_old':   dict(lit=(50, 4, 80), sh=(22, 4, 250)),
    'granite_fresh': dict(lit=(70, 3, 80), sh=(34, 3, 250)),
    # Columbia River Basalt: column faces 38-56 low sun, entablature 29 C15 h67 (1.3)
    'basalt':        dict(lit=(50, 14, 68), sh=(20, 8, 40)),
    'basalt_top':    dict(lit=(56, 22, 55), sh=(26, 12, 40)),    # red weathered flow top / bench
    # Olympic coast sandstone/mudstone: lit 21-24, shadow 14, olive tops 41 C6 h75 (1.4)
    'coast_rock':    dict(lit=(30, 5, 75), sh=(15, 4, 250)),
    # San Juan / Sucia sandstone: warm olive-tan lit 30 pockets 23 C12 h84 (1.5) — lifted for the lit face
    'sandstone':     dict(lit=(58, 16, 75), sh=(28, 10, 60)),
    'greenstone':    dict(lit=(52, 5, 150), sh=(18, 4, 220)),
    # Genge master copies: Ladram / West Bay palettes from the plates' own coloured pencil (my picks)
    'ladram_red':    dict(lit=(52, 38, 32), sh=(26, 22, 20)),
    'bridport_sand': dict(lit=(74, 38, 80), sh=(42, 24, 60)),
    # ground / vegetation used under the rock
    'grass_dry':     dict(lit=(66, 30, 88), sh=(34, 16, 110)),
    'grass_green':   dict(lit=(52, 30, 115), sh=(24, 16, 150)),
    'beach_sand':    dict(lit=(62, 10, 80), sh=(34, 6, 230)),
    'scree':         dict(lit=(58, 8, 70), sh=(28, 5, 240)),
}

SKY = {'day': (np.array([150, 185, 225]), np.array([205, 222, 238])),
       'overcast': (np.array([190, 196, 204]), np.array([222, 226, 230]))}


def ramp(mat, shift=(0.0, 0.0, 0.0), hour='day'):
    """8 sRGB colours (uint8). shift = (dL, dC, dh) applied to both anchors (unit exaggeration)."""
    m = MATS[mat] if isinstance(mat, str) else mat
    dL, dC, dh = shift
    lit = lch(m['lit'][0] + dL, max(0, m['lit'][1] + dC), m['lit'][2] + dh)
    sh = lch(m['sh'][0] + dL * 0.6, max(0, m['sh'][1] + dC * 0.6), m['sh'][2] + dh)
    if hour == 'overcast':   # compress families together: overcast light is sky-view, not sun
        lit = 0.55 * lit + 0.45 * sh + np.array([4, 0, 0])
    # steps: shadow 0..3 around sh (step 2), lit 4..7 around lit (step 6)
    Ls = []
    for s in range(8):
        if s < 4:
            base = sh.copy()
            base[0] = sh[0] + (s - 2) * 5.5
            if s == 0:
                base[0] = sh[0] - 13; base[1:] *= 0.7
        else:
            base = lit.copy()
            base[0] = lit[0] + (s - 6) * 7.0
            if s == 7:
                base[0] = lit[0] + 8; base[1:] *= 0.85
        Ls.append(base)
    Ls = np.array(Ls)
    return (lab2rgb(Ls[:, 0], Ls[:, 1], Ls[:, 2]) * 255 + 0.5).astype(np.uint8)


def unit_ramps(mat, n, spread=(6.0, 4.0, 10.0), seed=0):
    """n sibling ramps for n units of one rock, pushed apart (Genge: exaggerate unit differences)."""
    rng = np.random.default_rng(seed)
    out = []
    for k in range(n):
        s = np.array([rng.uniform(-1, 1) * spread[0], rng.uniform(-1, 1) * spread[1], rng.uniform(-1, 1) * spread[2]])
        out.append(ramp(mat, tuple(s)))
    return out


def sky_rgb(H, W, kind='day'):
    top, bot = SKY[kind]
    t = np.linspace(0, 1, H)[:, None, None]
    return np.broadcast_to(top * (1 - t) + bot * t, (H, W, 3)).astype(np.uint8)

"""geology.py — rock built from Genge's structural statistics, as heightfields with id fields.

Every builder returns a dict:
  h      heightfield (m) in its frame (face frame: protrusion toward viewer; terrain: elevation)
  ds     grid spacing (m)
  rock   bool mask of rock texels (else ground / background)
  bed    int id of the bed / sheet / flow each texel belongs to (-1 none)
  block  int id of the joint block (unique across the subject) (-1 none)
  unit   int lithological unit (indexes a palette list)
  comp   competence 0..1 per texel (1 = protrudes, ledge-former)
  bg     bool: background (sky) texels for a face; ground texels are rock=False, bg=False
Statistics follow MODEL.md §5: beds vary in thickness vertically and laterally, thin upward where
told to, pinch out; competent beds protrude; joints are few, perpendicular to bedding, and stop
(T-junction) at bed boundaries.
"""
import numpy as np
from scipy import ndimage as ndi


def _rng(seed):
    return np.random.default_rng(seed)


def lognormal(rng, mean, sigma, n):
    mu = np.log(mean) - sigma ** 2 / 2
    return rng.lognormal(mu, sigma, n)


def smooth_noise(n, scale_px, rng, amp=1.0):
    x = rng.standard_normal(n + 4 * int(scale_px) + 8)
    y = ndi.gaussian_filter1d(x, max(scale_px, 0.5))
    y = y[2 * int(scale_px) + 4: 2 * int(scale_px) + 4 + n]
    return amp * y / (y.std() + 1e-9)


def noise2(shape, scale_px, rng):
    a = ndi.gaussian_filter(rng.standard_normal(shape), scale_px)
    return a / (a.std() + 1e-9)


# ---------------------------------------------------------------- granite tor (Rough Tor, Fig 2.5/2.9)
def tor(width=12.0, height=7.0, ds=0.05, seed=3, sheet_mean=0.85, sheet_sigma=0.35,
        joint_mean=2.0, joint_sigma=0.45, round_v=0.45, undercut=1.6, round_h=0.25,
        gap=0.10, recess=0.35, taper=0.35, ground_h=0.6):
    """Stacked sheeting 'pillows' cut by vertical joints. round_v: fraction of each sheet's half-height
    that curves back (deep horizontal slots); undercut >1 makes the lower edge curve back more than the
    upper (the overhang Genge shades: 'shadow indicates top of overhang')."""
    rng = _rng(seed)
    W, H = int(width / ds), int((height + ground_h) / ds)
    v = (np.arange(H)[::-1] + 0.5) * ds - ground_h        # up, metres (row 0 = top of image)
    u = (np.arange(W) + 0.5) * ds
    h = np.full((H, W), -4.0)
    bed = -np.ones((H, W), int)
    block = -np.ones((H, W), int)
    # sheets from the ground up
    tops = [0.0]
    while tops[-1] < height:
        tops.append(tops[-1] + float(lognormal(rng, sheet_mean, sheet_sigma, 1)[0]))
    nb = 0
    cx = width / 2 + rng.uniform(-0.1, 0.1) * width
    for k in range(len(tops) - 1):
        v0, v1 = tops[k], min(tops[k + 1], height + 0.4)
        if v0 > height:
            break
        frac = v0 / height
        half = width * (0.5 - taper * frac ** 1.3 * 0.9) + rng.uniform(-0.8, 0.8)
        x0 = max(0.2, cx - half + rng.uniform(-0.7, 0.7) + frac * rng.uniform(-1.5, 1.5))
        x1 = min(width - 0.2, cx + half + rng.uniform(-0.7, 0.7) + frac * rng.uniform(-1.5, 1.5))
        if x1 - x0 < 1.0:
            break
        # joints in this sheet: offset per sheet (T-junctions at sheet boundaries)
        xs = [x0]
        while xs[-1] < x1:
            xs.append(xs[-1] + float(lognormal(rng, joint_mean, joint_sigma, 1)[0]))
        xs[-1] = x1
        if len(xs) > 2 and xs[-1] - xs[-2] < 0.6:
            xs.pop(-2)
        rows = (v >= v0) & (v < v1)
        tv = v1 - v0
        vv = (v[rows] - v0) / tv * 2 - 1                  # -1 bottom .. +1 top
        # pillow profile across height: superellipse, lower edge rounds harder (undercut)
        rv = np.where(vv < 0, round_v * undercut, round_v)
        prof_v = -np.clip(rv * tv * 0.5 * (1 - np.sqrt(np.clip(1 - np.abs(vv) ** 2.2, 0, 1))), 0, None)
        for b in range(len(xs) - 1):
            a0, a1 = xs[b], xs[b + 1]
            cols = (u >= a0) & (u < a1)
            if cols.sum() < 2:
                continue
            tw = a1 - a0
            uu = (u[cols] - a0) / tw * 2 - 1
            prof_h = -round_h * min(tw, 2.5) * 0.5 * (1 - np.sqrt(np.clip(1 - np.abs(uu) ** 3, 0, 1)))
            P = rng.uniform(-recess, recess) + 0.15 * np.sin(k * 1.7)
            tilt = rng.uniform(-0.08, 0.08)
            blockh = P + prof_v[:, None] + prof_h[None, :] + tilt * uu[None, :]
            # joint aperture: a narrow slot at the block's sides
            edge = np.minimum(u[cols] - a0, a1 - u[cols])
            slot = np.clip(1 - edge / max(gap, 1e-3), 0, 1) ** 2 * 0.6
            blockh = blockh - slot[None, :]
            sub = np.ix_(rows, cols)
            h[sub] = blockh
            bed[sub] = k
            block[sub] = nb
            nb += 1
    rock = bed >= 0
    # small surface relief: pitting / weathering hollows (texture, low amplitude)
    h = h + np.where(rock, 0.02 * noise2(h.shape, 2.0, rng), 0)
    # ground strip (grass) in front of the base, with a zigzag top
    gtop = 0.05 + 0.25 * np.abs(smooth_noise(W, 3, rng)) + 0.15 * rng.random(W)
    gmask = v[:, None] < gtop[None, :]
    h = np.where(gmask, 0.8, h)
    rock &= ~gmask
    bg = ~rock & ~gmask
    return dict(h=h, ds=ds, rock=rock, bed=np.where(rock, bed, -1), block=np.where(rock, block, -1),
                unit=np.where(rock, 0, -1), comp=np.where(rock, 1.0, 0.0), bg=bg, ground=gmask,
                u=u, v=v, width=width, height=height + ground_h)


# ---------------------------------------------------------------- bedded cliff face (West Bay, Olympic)
def beds_sequence(height, rng, bed_mean=0.6, sigma=0.5, thin_up=0.0, comp_frac=0.35, packets=None):
    """Bed tops from the base up. thin_up: fractional thinning of mean thickness from base to top
    (West Bay: 'cemented layers become more closely spaced and thinner towards the top', p. 98)."""
    tops, comp = [0.0], []
    while tops[-1] < height:
        f = tops[-1] / height
        m = bed_mean * (1 - thin_up * f)
        t = float(lognormal(rng, m, sigma, 1)[0])
        tops.append(tops[-1] + t)
        comp.append(1.0 if rng.random() < comp_frac else 0.0)
    return np.array(tops), np.array(comp)


def bedded_face(width=40.0, height=30.0, ds=0.1, seed=1, bed_mean=0.8, sigma=0.55, thin_up=0.5,
                comp_frac=0.3, protrude=0.35, soft_recess=0.25, dip_deg=0.0, taper=0.35, pinch=0.25,
                joint_mean=3.0, joint_sigma=0.5, joint_keep=0.5, joint_gap=0.08, bluffs=0,
                bluff_setback=4.0, furrows=0.15, skyline_rough=1.0, units=None, base_ground=1.5):
    """A cliff face of beds. Each bed boundary is a wandering line (lateral thickness change: taper);
    some boundaries die out (pinch). Competent beds protrude with a rounded lip and a shadowed soffit.
    Joints perpendicular to bedding, T-terminated at bed boundaries, only joint_keep of the candidate
    joints drawn through more than one bed. bluffs: the face is split into segments set back by
    bluff_setback (Genge §6.4: local horizons)."""
    rng = _rng(seed)
    W, H = int(width / ds), int((height + base_ground) / ds)
    u = (np.arange(W) + 0.5) * ds
    v = (np.arange(H)[::-1] + 0.5) * ds - base_ground
    U, V = np.meshgrid(u, v)
    tops, comp = beds_sequence(height + 2, rng, bed_mean, sigma, thin_up, comp_frac)
    nbnd = len(tops)
    # boundary k wanders laterally (thickness changes along the bed), plus dip
    Vb = np.zeros((nbnd, W))
    dip = np.tan(np.radians(dip_deg))
    for k in range(nbnd):
        wander = smooth_noise(W, 60 / ds * 0.1, rng, 1.0) * taper * (tops[min(k + 1, nbnd - 1)] - tops[max(k - 1, 0)]) * 0.25
        Vb[k] = tops[k] + dip * (u - width / 2) + wander
    Vb = np.maximum.accumulate(Vb, axis=0)              # boundaries never cross
    # bed index per texel
    bed = np.zeros((H, W), int)
    for k in range(1, nbnd):
        bed += (V >= Vb[k][None, :]).astype(int)
    # pinch-outs: some boundaries fade (the two beds merge over a stretch)
    alive = np.ones((nbnd, W), bool)
    for k in range(1, nbnd):
        if rng.random() < pinch:
            c = rng.uniform(0, width)
            L = rng.uniform(0.15, 0.5) * width
            alive[k] = np.abs(u - c) > L / 2 if rng.random() < 0.5 else (u - c) * rng.choice([-1, 1]) < 0
    # profile within each bed: position 0..1 from base to top
    b0 = np.take_along_axis(Vb, np.clip(bed, 0, nbnd - 1), 0)
    b1 = np.take_along_axis(Vb, np.clip(bed + 1, 0, nbnd - 1), 0)
    t = np.clip((V - b0) / np.maximum(b1 - b0, 1e-3), 0, 1)
    cb = comp[np.clip(bed, 0, len(comp) - 1)]
    # competent: protrude with rounded lip and undercut base; soft: recess, weathered into a groove
    lip = np.sin(np.pi * np.clip(t, 0, 1)) ** 0.35
    h_comp = protrude * lip - 0.15 * protrude * (1 - t) ** 6        # undercut near the bed's base
    h_soft = -soft_recess * np.sin(np.pi * t) ** 0.8 * 0.8
    h = np.where(cb > 0.5, h_comp, h_soft)
    # per-bed set back / stand out a little, varying along strike
    setb = np.array([smooth_noise(W, 8 / ds, rng, 0.08) for _ in range(nbnd)])
    h = h + np.take_along_axis(setb, np.clip(bed, 0, nbnd - 1), 0)
    # joints: candidate positions per bed (perpendicular to bedding = vertical here), T-terminated
    block = np.zeros((H, W), int)
    jmask = np.zeros((H, W), bool)
    nb = 0
    through = []   # joints that run through several beds (the few Genge would draw)
    for k in range(nbnd):
        xs = np.cumsum(lognormal(rng, joint_mean, joint_sigma, int(width / joint_mean * 3) + 4)) - rng.uniform(0, joint_mean)
        xs = xs[(xs > 0) & (xs < width)]
        rows = bed == k
        if not rows.any():
            continue
        edges = np.concatenate([[0], xs, [width]])
        idx = np.searchsorted(edges, u) - 1
        block[rows] = (nb + idx[None, :].repeat(H, 0))[rows]
        nb += len(edges)
        for x in xs:
            j = int(x / ds)
            jmask[rows & (np.abs(U - x) < joint_gap * 0.5 + ds * 0.5)] = True
    # open joints are slots
    jd = ndi.distance_transform_edt(~jmask) * ds
    h = h - 0.35 * np.clip(1 - jd / (joint_gap * 1.5 + 1e-3), 0, 1)
    # bluffs: segments set back; a bluff edge makes a local horizon when viewed obliquely
    seg = np.zeros(W, int)
    if bluffs > 0:
        cuts = np.sort(rng.uniform(0.12, 0.88, bluffs)) * width
        seg = np.searchsorted(cuts, u)
        off = np.cumsum(rng.uniform(0.5, 1.0, bluffs + 1) * bluff_setback)
        off -= off[0]
        # rounded bluff noses
        nose = np.zeros(W)
        for c in cuts:
            nose += -1.5 * np.exp(-((u - c) / 1.2) ** 2)
        h = h - off[seg][None, :] + nose[None, :] * 0.0
        h = h + np.where(True, 0, 0)
    # vertical furrows (weathering, p. 98): shallow grooves that break the lit faces
    if furrows > 0:
        fz = smooth_noise(W, 0.6 / ds, rng, 1.0)
        fz = np.clip(fz, 0, None) ** 2
        h = h - furrows * 0.25 * fz[None, :] * (0.6 + 0.4 * noise2((H, W), 6, rng))
    # skyline: top of face rough, stepped over bluffs
    top = height - np.abs(smooth_noise(W, 3 / ds, rng, skyline_rough)) - (seg * rng.uniform(0.2, 1.2))
    bg = V > top[None, :]
    ground = V < 0.02 + 0.3 * np.abs(smooth_noise(W, 2 / ds, rng, 1))[None, :]
    h = np.where(ground, 1.2 + 0.2 * noise2((H, W), 3, rng), h)
    h = np.where(bg, -40.0, h)
    rock = ~bg & ~ground
    unit = np.zeros((H, W), int)
    if units is not None:        # units: list of (from_height, unit_id)
        for z0, uid in units:
            unit[V >= z0] = uid
    return dict(h=h, ds=ds, rock=rock, bed=np.where(rock, bed, -1), block=np.where(rock, block, -1),
                unit=np.where(rock, unit, -1), comp=np.where(rock, cb, 0.0), bg=bg, ground=ground,
                alive=alive, Vb=Vb, u=u, v=v, width=width, height=height + base_ground, seg=seg)

"""genge.py — Genge's drawing rules as a pixel painter.

Input: a screen-space G-buffer dict `g` (see render.gbuffer). Output: an index canvas
(step 0..7 per pixel, -1 for non-rock) plus a colour image through per-unit palettes.

Each rule is a parameter, so any one can be switched off to see what it carries:

  shadows   'photo' = sun + all cast shadows (what a camera records)
            'genge' = no cast shadows at all (Genge §2.9.3: 'should not be included in a sketch')
            'mass'  = cast shadows kept only where they are big masses (>= mass_frac of the rock);
                      small ones dropped, because they obscure structure (the pixel translation)
  ao_inset  step drop in inset areas (Genge's first, grey, layer)            'inset areas shadowed'
  ao_deep   occlusion above which a pixel is forced to the darkest step     'darkest shadow in most inset areas'
  soffit    1-px darkest run under a protruding bed (the 'thick line' on    'shadow indicates top of overhang'
            competent beds with shadowed overhangs, p. 106)
  lip       1-px highlight on the top edge of a protruding lit bed          'highlights and lowlights along the top
                                                                             and bottom of beds' (p. 279)
  bed_lines impersistent 1-step-darker runs on bed boundaries, run length   'discontinuous lines … stroke lengths that
            scaled to the bed's thickness in pixels                          reflect the scale of the objects' (p. 28)
  joints    joint lines, but only the longest, and never more joint pixels  'not obscure the important features by
            than joint_ratio x bedding-line pixels                           adding too many' (p. 61; Fig 4.6)
  contact   darkest 1-px pocket where rock meets the ground in front        'thick line emphasizes foreground'
  skyline   clean silhouette; 1 step darker edge only where rock and sky    'thick line emphasizes horizon'
            are close in value
  streak    texture streaked along the bedding (value noise elongated       Fig 2.8c, Fig 16.7c
            along the bed), amplitude in steps
"""
import numpy as np
from scipy import ndimage as ndi

DEFAULT = dict(shadows='mass', mass_frac=0.06, ao_inset=1.2, ao_lo=0.18, ao_hi=0.45, ao_deep=0.55,
               soffit=True, lip=True, bed_lines=True, bed_run=(0.9, 0.35), joints=True, joint_ratio=0.35,
               joint_min_px=4, contact=True, skyline=True, streak=0.35, streak_len=6.0, seam=0.25,
               up=(0, 1, 0), lit_gamma=0.8, seed=0)

BAYER2 = np.array([[0.0, 0.5], [0.75, 0.25]])


def _bayer(H, W):
    return np.tile(BAYER2, (H // 2 + 1, W // 2 + 1))[:H, :W]


def _hash(*a):
    x = np.uint64(1469598103934665603)
    out = np.zeros(np.broadcast(*a).shape, np.uint64) + x
    for v in a:
        out = (out ^ np.asarray(v).astype(np.uint64)) * np.uint64(1099511628211)
    return (out % np.uint64(1000003)).astype(np.float64) / 1000003.0


def families(g, L, p):
    N = g['N']
    ndl = (N * np.asarray(L)[None, None, :]).sum(-1)
    rock = g['rock']
    self_sh = ndl <= 0.04
    cast = g['shadow'] > 0.5
    if p['shadows'] == 'photo':
        sh = self_sh | cast
    elif p['shadows'] == 'genge':
        sh = self_sh
    else:   # 'mass': keep only big connected cast-shadow masses
        lab, n = ndi.label(cast & rock)
        keep = np.zeros_like(cast)
        if n:
            sizes = ndi.sum(np.ones_like(lab), lab, np.arange(1, n + 1))
            big = np.flatnonzero(sizes >= p['mass_frac'] * max(rock.sum(), 1)) + 1
            keep = np.isin(lab, big)
        sh = self_sh | keep
    lit = rock & ~sh
    return lit, ndl


def value(g, L, p):
    lit, ndl = families(g, L, p)
    up = np.asarray(p['up'], float)
    sky = 0.5 + 0.5 * (g['N'] * up[None, None, :]).sum(-1)
    v_lit = 4.0 + 3.0 * np.clip(ndl, 0, 1) ** p['lit_gamma']
    v_sh = 1.1 + 2.1 * sky
    v = np.where(lit, v_lit, v_sh)
    ao = g['ao']
    # first layer (grey): inset areas lose up to ao_inset steps, staying in their family
    t = np.clip((ao - p['ao_lo']) / max(p['ao_hi'] - p['ao_lo'], 1e-3), 0, 1)
    floor = np.where(lit, 4.0, 0.6)
    v = np.maximum(v - p['ao_inset'] * t, floor)
    return v, lit


def streaks(shape, along, length, seed):
    """Value noise elongated along a direction field (angle radians per pixel, or a constant)."""
    rng = np.random.default_rng(seed)
    n = rng.standard_normal(shape)
    if np.isscalar(along):
        ang = np.degrees(along)
        k = np.zeros((1, int(length) * 2 + 1)); k[:] = 1
        n2 = ndi.convolve(n, k / k.sum(), mode='wrap')
        if abs(ang) > 1:
            n2 = ndi.rotate(n2, ang, reshape=False, mode='wrap')
    else:
        n2 = ndi.uniform_filter1d(n, int(length) * 2 + 1, axis=1)
    n2 = ndi.gaussian_filter(n2, 0.5)
    return n2 / (n2.std() + 1e-9)


def boundary(ids, mask, axis):
    """Pixels whose id differs from the neighbour below (axis=0) or to the right (axis=1)."""
    b = np.zeros_like(mask)
    if axis == 0:
        d = (ids[:-1] != ids[1:]) & mask[:-1] & mask[1:]
        b[:-1] |= d
    else:
        d = (ids[:, :-1] != ids[:, 1:]) & mask[:, :-1] & mask[:, 1:]
        b[:, :-1] |= d
    return b


def impersistent(line, run_px, gap_frac, seed):
    """Break a 1-px line mask into runs: walk connected pieces, keep dashes of ~run_px with gaps."""
    lab, n = ndi.label(line, structure=np.ones((3, 3)))
    if n == 0:
        return line
    ys, xs = np.nonzero(line)
    lid = lab[ys, xs]
    order = np.lexsort((ys, xs, lid))
    ys, xs, lid = ys[order], xs[order], lid[order]
    out = np.zeros_like(line)
    rng = np.random.default_rng(seed)
    start = 0
    for k in range(1, n + 1):
        end = np.searchsorted(lid, k, side='right')
        m = end - start
        pos = 0
        while pos < m:
            r = max(2, int(rng.lognormal(np.log(max(run_px, 2)), 0.4)))
            gp = max(1, int(r * gap_frac * rng.uniform(0.5, 1.5)))
            sel = slice(start + pos, start + min(pos + r, m))
            out[ys[sel], xs[sel]] = True
            pos += r + gp
        start = end
    return out


def paint(g, L, params=None):
    p = dict(DEFAULT)
    if params:
        p.update(params)
    H, W = g['rock'].shape
    rock = g['rock']
    v, lit = value(g, L, p)
    # texture streaked along the bedding (horizontal here; faces are drawn with bedding ~ horizontal)
    if p['streak'] > 0:
        s = streaks((H, W), p.get('bed_angle', 0.0), p['streak_len'], p['seed'] + 7)
        v = v + p['streak'] * 0.5 * s
        v = np.where(lit, np.clip(v, 4.0, 7.4), np.clip(v, 0.6, 3.49))
    # quantize: flat steps, 2x2 dither only in a narrow seam, never across the family bar
    B = _bayer(H, W)
    frac = v - np.floor(v)
    step = np.floor(v).astype(int)
    seam = np.abs(frac - 0.5) < p['seam']
    step = np.where(seam, np.floor(v - 0.5 + B + 0.0).astype(int), np.round(v).astype(int))
    step = np.where(lit, np.clip(step, 4, 7), np.clip(step, 0, 3))
    ao = g['ao']
    # second layer (black): the most inset areas are the darkest step, whatever the light
    if p['ao_deep'] < 1:
        deep = rock & (ao > p['ao_deep'])
        step = np.where(deep, np.where(ao > p['ao_deep'] + 0.12, 0, 1), step)
    lines = {}
    bed = g['bed']
    comp = g['comp']
    below_bed = np.roll(bed, -1, 0)
    bnd = boundary(bed, rock, 0)                     # pixel above a bed boundary
    bnd_below = np.roll(bnd, 1, 0) & rock           # pixel below it
    # soffit: the pixel under a protruding bed's base = darkest (thick line on competent beds)
    if p['soffit']:
        comp_above = np.roll(comp, 1, 0) > 0.5
        sof = bnd_below & comp_above & (comp < 0.5)
        step = np.where(sof, np.minimum(step, 1 if p['shadows'] != 'photo' else 0), step)
        step = np.where(sof & ~lit, 0, step)
        lines['soffit'] = sof
    # lip: top edge of a protruding bed catches light
    if p['lip']:
        comp_here = comp > 0.5
        lipm = bnd_below & comp_here & lit & (np.roll(comp, 1, 0) < 0.5)
        step = np.where(lipm, np.minimum(step + 1, 7), step)
        lines['lip'] = lipm
    # impersistent bed lines on the remaining boundaries
    bl_px = 0
    if p['bed_lines']:
        other = bnd & ~lines.get('soffit', np.zeros_like(bnd)) & ~np.roll(lines.get('lip', np.zeros_like(bnd)), -1, 0)
        if 'alive' in g:
            other &= g['alive']
        th = g.get('bed_px', None)
        run = p['bed_run'][0] * (np.median(th[rock]) if th is not None and rock.any() else 6.0)
        other = impersistent(other, max(3, run * 2.0), p['bed_run'][1] * 2, p['seed'])
        step = np.where(other, np.maximum(step - 1, np.where(lit, 4, 0)), step)
        lines['bed'] = other
        bl_px = other.sum() + lines.get('soffit', np.zeros(1)).sum()
    # joints: vertical boundaries between blocks inside one bed, longest first, capped
    if p['joints']:
        jb = boundary(g['block'], rock, 1) & ~bnd
        lab, n = ndi.label(jb, structure=np.array([[0, 1, 0], [0, 1, 0], [0, 1, 0]]))
        keep = np.zeros_like(jb)
        if n:
            sizes = ndi.sum(np.ones_like(lab), lab, np.arange(1, n + 1))
            order = np.argsort(-sizes)
            budget = p['joint_ratio'] * max(bl_px, 1) if p['bed_lines'] else np.inf
            used = 0
            for k in order:
                if sizes[k] < p['joint_min_px'] or used + sizes[k] > budget:
                    continue
                keep |= lab == k + 1
                used += sizes[k]
        step = np.where(keep, np.maximum(step - 1, 0), step)
        lines['joint'] = keep
    # contact with the ground in front: darkest pocket (thick foreground line)
    if p['contact'] and 'ground' in g:
        gd = g['ground']
        con = rock & np.roll(gd, -1, 0) & ~np.roll(rock, -1, 0)
        con |= rock & ndi.binary_dilation(gd, np.ones((1, 3))) & ~gd & np.roll(gd, -1, 0)
        step = np.where(con, 0, step)
        lines['contact'] = con
    step = np.where(rock, step, -1)
    return step, lit, lines, p


def colorize(step, g, ramps, ground_ramp=None, sky=None, ground_step=None, skyline=True):
    H, W = step.shape
    out = np.zeros((H, W, 3), np.uint8)
    if sky is not None:
        out[:] = sky
    unit = g['unit']
    for k, r in enumerate(ramps):
        m = (unit == k) & (step >= 0)
        out[m] = r[step[m]]
    if ground_ramp is not None and ground_step is not None:
        m = g['ground'] & g['hit']
        out[m] = ground_ramp[ground_step[m]]
    if skyline and sky is not None:
        # skyline: clean edge; darken the edge pixel 1 step only where rock is close to the sky in value
        rock = step >= 0
        edge = rock & ndi.binary_dilation(g['bg'], np.ones((3, 3))) & ~g['bg']
        lum = out.astype(float) @ np.array([0.3, 0.59, 0.11])
        skyl = np.asarray(sky, float).reshape(-1, 3).mean(0) @ np.array([0.3, 0.59, 0.11]) if np.ndim(sky) <= 1 \
            else (np.asarray(sky, float) @ np.array([0.3, 0.59, 0.11]))
        close = np.abs(lum - skyl) < 40
        for k, r in enumerate(ramps):
            m = edge & close & (unit == k)
            out[m] = r[np.maximum(step[m] - 1, 0)]
    return out


def ground_steps(g, L, p, rng_seed=0):
    """Simple ground (grass) steps: lit by N.L, darker under the rock's shadow, zigzag top kept."""
    N = g['N']
    ndl = (N * np.asarray(L)[None, None, :]).sum(-1)
    v = 4.5 + 2.0 * np.clip(ndl, 0, 1)
    rng = np.random.default_rng(rng_seed)
    v = v + 0.8 * rng.standard_normal(ndl.shape) * 0.5
    if p['shadows'] != 'genge':
        v = np.where(g['shadow'] > 0.5, 2.2, v)
    return np.clip(np.round(v), 0, 7).astype(int)

"""render.py — geology dict + frame + camera + light -> screen G-buffer dict for genge.paint."""
import numpy as np
import engine as E


def run_length_axis0(ids):
    """Per texel: length (texels) of the vertical run of equal ids it sits in."""
    H, W = ids.shape
    out = np.zeros((H, W), float)
    for j in range(W):
        col = ids[:, j]
        ch = np.flatnonzero(np.diff(col) != 0) + 1
        st = np.concatenate([[0], ch]); en = np.concatenate([ch, [H]])
        for a, b in zip(st, en):
            out[a:b, j] = b - a
    return out


class Scene:
    """Texture-space terms are cached per (geology, light, ao radius)."""
    def __init__(self, geo, frame, L_world, ao_radius=1.5, up_world=(0, 0, 1)):
        self.geo, self.frame = geo, frame
        h, ds = geo['h'], geo['ds']
        self.N = E.normals(h, ds)
        L = np.asarray(L_world, float); L /= np.linalg.norm(L)
        self.L_world = L
        self.L_local = np.array([L @ frame.E1, L @ frame.E2, L @ frame.N])
        self.up_local = np.array([np.dot(up_world, frame.E1), np.dot(up_world, frame.E2), np.dot(up_world, frame.N)])
        # heightfield rows run top->bottom but E2 points up for faces: flip i so that i grows along E2
        self.hf = h[::-1].copy()
        self.ao = E.sky_occlusion(self.hf, ds, ao_radius)[::-1]
        self.sh = E.sun_shadow(self.hf, ds, self.L_local)[::-1]
        self.bedlen = run_length_axis0(geo['bed']) * ds if 'bed' in geo else None

    def gbuffer(self, cam, mpp=None):
        geo = self.geo
        hi, hj, dep, hit = E.march(self.hf, geo['ds'], self.frame, cam)
        H = geo['h'].shape[0]
        g = E.GBuf(H - 1 - hi, hj, dep, hit, geo['ds'])
        g.hit = hit
        out = dict(hit=hit, depth=dep)
        N = self.N.copy()
        N[..., 1] *= -1            # normals() was computed with i growing downward; E2 is up
        out['N'] = g.take(N, 0.0)
        out['ao'] = g.take(self.ao, 0.0)
        out['shadow'] = g.take(self.sh, 0.0)
        for k in ('rock', 'ground', 'bg'):
            out[k] = g.take(geo[k], False) & hit
        out['bg'] = out['bg'] | ~hit
        for k in ('bed', 'block', 'unit'):
            out[k] = g.take(geo[k], -1)
        out['comp'] = g.take(geo['comp'], 0.0)
        out['h'] = g.take(geo['h'], 0.0)
        if 'alive_tex' in geo:
            out['alive'] = g.take(geo['alive_tex'], True)
        if self.bedlen is not None and mpp:
            out['bed_px'] = g.take(self.bedlen, 0.0) / mpp
        out['L_local'] = self.L_local
        out['up_local'] = self.up_local
        return out


def face_frame(width, height):
    """Face frame: E1 = east (x), E2 = up (z), N = toward viewer (-y). Origin at the face's lower left."""
    return E.Frame(O=(0, 0, 0), E1=(1, 0, 0), E2=(0, 0, 1), N=(0, -1, 0))


def front_camera(geo, mpp, elev_deg=0.0, az_deg=0.0, pad=0.0):
    """Orthographic camera looking at the face from the front (optionally from below/above, or aside)."""
    W_m, H_m = geo['width'], geo['height']
    Wpx = int(round((W_m + 2 * pad) / mpp))
    Hpx = int(round((H_m + 2 * pad) / mpp))
    ctr = np.array([W_m / 2, 0.0, H_m / 2 - (geo.get('ground_h', 0) or 0)])
    e, a = np.radians(elev_deg), np.radians(az_deg)
    d = np.array([np.sin(a) * np.cos(e), np.cos(a) * np.cos(e), -np.sin(e)])   # direction of view (into +y)
    eye = ctr - d * 200
    base_v = geo['v'].min()
    ctr[2] = base_v + H_m / 2
    eye = ctr - d * 200
    return E.Camera(eye, ctr, Wpx, Hpx, ortho_w=W_m + 2 * pad)

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s = open('render.py').read()
a = s.index('def face_frame')
s = s[:a] + '''def face_frame(geo):
    """Face frame: E1 = east (x), E2 = up (z), N = toward viewer (-y). Origin at the face's lower left."""
    return E.Frame(O=(0, 0, float(geo['v'].min())), E1=(1, 0, 0), E2=(0, 0, 1), N=(0, -1, 0))


def front_camera(geo, mpp, elev_deg=0.0, az_deg=0.0, pad=0.0):
    """Orthographic camera looking at the face from the front (optionally from below/above, or aside)."""
    W_m, H_m = geo['width'], geo['height']
    Wpx = int(round((W_m + 2 * pad) / mpp))
    Hpx = int(round((H_m + 2 * pad) / mpp))
    ctr = np.array([W_m / 2, 0.0, float(geo['v'].min()) + H_m / 2])
    e, a = np.radians(elev_deg), np.radians(az_deg)
    d = np.array([np.sin(a) * np.cos(e), np.cos(a) * np.cos(e), -np.sin(e)])
    return E.Camera(ctr - d * 200, ctr, Wpx, Hpx, ortho_w=W_m + 2 * pad)
'''
open('render.py','w').write(s)
EOF
echo ok

# output:
ok

cd ~/work/pixelart-studies/books/genge; cat > src/e1_tor.py <<'EOF'
"""E1b. Rough Tor: photo-light vs Genge rules, at close / landscape / region."""
import sys, numpy as np
sys.path.insert(0, 'src')
import geology, render, genge, palette
from util import load, area_down, up, sheet

geo = geology.tor(seed=int(sys.argv[1]) if len(sys.argv) > 1 else 3)
fr = render.face_frame(geo)
L = (-0.55, -0.55, 0.63)
sc = render.Scene(geo, fr, L, ao_radius=1.2)
ramp = palette.ramp('tor_granite')
gr = palette.ramp('grass_dry')
ph = load('refs/Rough_Tor_Bodmin_Moor_granite_tor/commons_Roughtor_Bodmin_Moor_2_jpg_76d7aed8.jpg')[60:853, 20:1280]
fig = load('plates/crops/fig2_9b_roughtor_ao.png')[60:, :]
rows = []
for lab, mpp in [('close 0.1', 0.1), ('landscape 0.35', 0.35), ('region 1.0', 1.0)]:
    cam = render.front_camera(geo, mpp)
    g = sc.gbuffer(cam, mpp)
    sky = palette.sky_rgb(cam.H, cam.W)
    r = []
    k = max(1, int(round(360 / cam.W)))
    r.append(('photo', up(area_down(ph, int(9.0 / mpp)), k)))
    r.append(('Genge 2.9b', up(area_down(fig, int(12.0 / mpp)), k)))
    for name, prm in [('photo light', dict(shadows='photo', ao_inset=0, ao_deep=2, soffit=False, lip=False,
                                           bed_lines=False, joints=False, contact=False, streak=0.2)),
                      ('Genge rules', dict(shadows='mass'))]:
        step, lit, lines, p = genge.paint(g, sc.L_local, prm)
        gs = genge.ground_steps(g, sc.L_local, p)
        img = genge.colorize(step, g, [ramp], gr, sky, gs, skyline=prm.get('skyline', True))
        r.append((f'{name} {lab} m/px', up(img, k)))
    rows.append(r)
sheet(rows, title='E1b  Rough Tor — photo | Genge Fig 2.9b | photo-light render | Genge-rules render').save('stages/e1/e1b_tor.png')
EOF
time nice -n 10 python3 src/e1_tor.py

# output:
nice -n 10 python3 src/e1_tor.py  5.45s user 0.59s system 21% cpu 27.579 total

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s = open('geology.py').read()
a = s.index('# ---------------------------------------------------------------- granite tor')
b = s.index('# ---------------------------------------------------------------- bedded cliff face')
new = '''# ---------------------------------------------------------------- granite tor (Rough Tor, Fig 2.5/2.9)
def tor(width=12.0, height=6.5, ds=0.05, seed=3, sheet_mean=0.95, sheet_sigma=0.4, bend=0.25,
        joint_mean=2.6, joint_sigma=0.6, round_v=0.95, undercut=1.35, end_round=0.9,
        gap=0.12, recess=0.45, taper=0.3, stacks=2, ground_h=0.6):
    """Stacked sheeting 'pillows' cut by vertical joints (Rough Tor). Sheet boundaries are curved
    (sheeting joints follow the old dome: bend), blocks are lens-shaped: round_v of each sheet's half
    height curves back (deep slots), the lower edge more (undercut: 'shadow indicates top of
    overhang'); block ends round back in plan (end_round). The outline is one or more stacks."""
    rng = _rng(seed)
    W, H = int(width / ds), int((height + ground_h) / ds)
    v = (np.arange(H)[::-1] + 0.5) * ds - ground_h
    u = (np.arange(W) + 0.5) * ds
    U, V = np.meshgrid(u, v)
    tops = [0.0]
    while tops[-1] < height + 1:
        tops.append(tops[-1] + float(lognormal(rng, sheet_mean, sheet_sigma, 1)[0]))
    nbnd = len(tops)
    Vb = np.zeros((nbnd, W))
    for k in range(nbnd):
        Vb[k] = tops[k] + bend * smooth_noise(W, 2.5 / ds, rng, 1.0) * (0.3 + 0.7 * k / nbnd)
    Vb = np.maximum.accumulate(Vb, axis=0)
    Vb[:, :] = np.maximum(Vb, Vb[max(0, 0)] * 0)
    bed = np.zeros((H, W), int)
    for k in range(1, nbnd):
        bed += (V >= Vb[k][None, :]).astype(int)
    b0 = np.take_along_axis(Vb, np.clip(bed, 0, nbnd - 1), 0)
    b1 = np.take_along_axis(Vb, np.clip(bed + 1, 0, nbnd - 1), 0)
    tv = np.maximum(b1 - b0, 1e-3)
    t = np.clip((V - b0) / tv, 0, 1) * 2 - 1                       # -1 bottom .. +1 top of sheet
    rv = np.where(t < 0, round_v * undercut, round_v)
    prof_v = -rv * tv * 0.5 * (1 - np.sqrt(np.clip(1 - np.abs(t) ** 2.0, 0, 1)))
    # stacks: centres and half-widths; each sheet's extent within its stack shrinks upward
    cents = np.sort(rng.uniform(0.25, 0.75, stacks)) * width
    hw0 = rng.uniform(0.35, 0.5, stacks) * width / max(1, stacks) * 1.6
    hts = rng.uniform(0.6, 1.0, stacks) * height; hts[rng.integers(stacks)] = height
    h = np.full((H, W), -4.0)
    block = -np.ones((H, W), int)
    inside = np.zeros((H, W), bool)
    nb = 0
    for k in range(nbnd - 1):
        rows_k = bed == k
        for s in range(stacks):
            if tops[k] > hts[s]:
                continue
            f = tops[k] / hts[s]
            half = hw0[s] * (1 - taper * f ** 1.5) + rng.uniform(-0.5, 0.5)
            c = cents[s] + rng.uniform(-0.6, 0.6)
            x0, x1 = max(0.3, c - half), min(width - 0.3, c + half)
            if x1 - x0 < 0.8:
                continue
            xs = [x0]
            while xs[-1] < x1:
                xs.append(xs[-1] + float(lognormal(rng, joint_mean, joint_sigma, 1)[0]))
            xs[-1] = x1
            if len(xs) > 2 and xs[-1] - xs[-2] < 0.5:
                xs.pop(-2)
            for bi in range(len(xs) - 1):
                a0, a1 = xs[bi], xs[bi + 1]
                cols = (u >= a0) & (u < a1)
                m = rows_k & cols[None, :]
                if m.sum() < 4:
                    continue
                tw = a1 - a0
                uu = (U[m] - a0) / tw * 2 - 1
                # ends at the stack edge round back strongly (lens), interior joints less
                er = np.where(uu < 0, end_round if bi == 0 else 0.35, end_round if bi == len(xs) - 2 else 0.35)
                prof_h = -er * min(tw, 2.0) * 0.5 * (1 - np.sqrt(np.clip(1 - np.abs(uu) ** 2.5, 0, 1)))
                P = rng.uniform(-recess, recess) + 0.4 * (1 - f)
                hh = P + prof_v[m] + prof_h + rng.uniform(-0.1, 0.1) * uu
                # the lens also thins toward its ends (sheets pinch at the stack edge)
                pinch = (np.abs(uu) ** 6) * (0.5 if bi in (0, len(xs) - 2) else 0.0)
                keep = np.abs(t[m]) < 1 - pinch
                mm = np.zeros_like(m); mm[m] = keep
                h[mm] = np.maximum(h[mm], hh[keep])
                block[mm] = nb
                inside |= mm
                nb += 1
    rock = inside
    h = h + np.where(rock, 0.025 * noise2(h.shape, 2.0, rng), 0)
    gtop = 0.05 + 0.25 * np.abs(smooth_noise(W, 3, rng)) + 0.15 * rng.random(W)
    gmask = V < gtop[None, :]
    h = np.where(gmask, 0.8, h)
    rock &= ~gmask
    bg = ~rock & ~gmask
    return dict(h=h, ds=ds, rock=rock, bed=np.where(rock, bed, -1), block=np.where(rock, block, -1),
                unit=np.where(rock, 0, -1), comp=np.where(rock, 1.0, 0.0), bg=bg, ground=gmask,
                u=u, v=v, width=width, height=height + ground_h)


'''
s = s[:a] + new + s[b:]
open('geology.py','w').write(s)
EOF
cd ..; python3 - <<'EOF'
s=open('src/genge.py').read()
s=s.replace("""    if p['contact'] and 'ground' in g:""","""    if p['contact'] and 'ground' in g and p.get('mpp', 0.1) <= p.get('contact_max_mpp', 0.4):""")
open('src/genge.py','w').write(s)
s=open('src/e1_tor.py').read()
s=s.replace("('Genge rules', dict(shadows='mass'))","('Genge rules', dict(shadows='mass', mpp=mpp))")
open('src/e1_tor.py','w').write(s)
EOF
nice -n 10 python3 src/e1_tor.py 3

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s = open('render.py').read()
s = s.replace("""    def gbuffer(self, cam, mpp=None):
        geo = self.geo""", """    def gbuffer(self, cam, mpp=None, ss=1):
        \"\"\"ss > 1: supersample ss x ss rays per pixel; floats are area-averaged (an honest pixel sees
        the mean of what falls in it: a 0.1 m joint at 0.35 m/px is a grey hint, not a black line),
        booleans by majority, ids by the centre sample.\"\"\"
        if ss > 1:
            import copy
            c2 = copy.copy(cam); c2.W, c2.H = cam.W * ss, cam.H * ss
            G = self.gbuffer(c2, mpp / ss if mpp else None, 1)
            out = {}
            c = ss // 2
            for k, a in G.items():
                if not isinstance(a, np.ndarray) or a.ndim < 2 or a.shape[0] != c2.H:
                    out[k] = a; continue
                sh = (cam.H, ss, cam.W, ss) + a.shape[2:]
                if a.dtype == bool:
                    out[k] = a.reshape(sh).mean((1, 3)) >= 0.5
                elif np.issubdtype(a.dtype, np.integer):
                    out[k] = a[c::ss, c::ss]
                else:
                    out[k] = a.reshape(sh).mean((1, 3))
            if 'N' in out:
                out['N'] = out['N'] / np.maximum(np.linalg.norm(out['N'], axis=-1, keepdims=True), 1e-6)
            if 'bed_px' in out:
                out['bed_px'] = out['bed_px'] / ss
            return out
        geo = self.geo""")
open('render.py','w').write(s)
s = open('e1_tor.py').read()
s = s.replace("g = sc.gbuffer(cam, mpp)", "g = sc.gbuffer(cam, mpp, ss=int(sys.argv[2]) if len(sys.argv) > 2 else 3)")
s = s.replace("geo = geology.tor(seed=int(sys.argv[1]) if len(sys.argv) > 1 else 3)",
 "geo = geology.tor(seed=int(sys.argv[1]) if len(sys.argv) > 1 else 3, sheet_mean=0.75, joint_mean=2.2, recess=0.5, stacks=3, height=7.0)")
open('e1_tor.py','w').write(s)
EOF
cd ..; nice -n 10 python3 src/e1_tor.py 5 3

# output:
[This command modified 1 file you've previously read: ../books/genge/src/render.py. Call Read before editing.]

cd ~/work/pixelart-studies/books/genge; ls refs/West_Bay_Bridport_cliffs/ refs/Kalaloch_sandstone_cliff refs/Sucia_Island_sandstone_honeycomb

# output:
refs/Kalaloch_sandstone_cliff:
_sheet.jpg
openverse_Mobile_Homes_7f28157b.jpg
openverse_Sandstone_rocks_tilted_almost_90_degrees_Beach_3_9137b1ad.jpg
openverse_rock_8cf8d475.jpg
sources.tsv

refs/Sucia_Island_sandstone_honeycomb:
_sheet.jpg
commons_Folded_Bedrock_Sucia_Island_and_Patos_Island_jpg_1f01218f.jpg
commons_Sucia_in_Summer_153693_jpg_c244fb78.jpg
openverse_20100804_DSCF0398_jpg_e0b256e4.jpg
openverse_Layers_of_Texture_b25a30d6.jpg
openverse_Sea_Lions_with_Mt_Baker_Sucia_Island_Kayak_Trip_53645ff3.jpg
openverse_Still_Life_on_Sucia_Island_31b532cb.jpg
sources.tsv

refs/West_Bay_Bridport_cliffs/:
_sheet.jpg
commons_Cliffs_West_Bay_Bridport_geograph_org_uk_5494083_9741edc9.jpg
commons_EastCliff_WestBay_jpg_c7f93602.jpg
commons_West_Bay_The_beach_and_Bridport_Sands_cliffs_geo_108bc900.jpg
commons_West_Bay_The_beach_and_Bridport_Sands_cliffs_geo_836a9cc4.jpg
openverse_16_023_West_Bay_Station_one_of_five_looking_down_670f2a82.jpg
openverse_West_Bay_Bridport_9ae7c59a.jpg
openverse_West_Bay_East_Beach_Beach_Safety_sign_6d614994.jpg
openverse_West_Bay_East_Beach_b6991a6b.jpg
sources.tsv

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s = open('geology.py').read()
a = s.index("    # bed index per texel\n    bed = np.zeros((H, W), int)")
b = s.index("    cb = comp[np.clip(bed, 0, len(comp) - 1)]")
new = '''    # pinch-outs: some boundaries die out along strike (the two beds merge over a stretch)
    alive = np.ones((nbnd, W), bool)
    for k in range(1, nbnd):
        if rng.random() < pinch:
            c = rng.uniform(0, width)
            L = rng.uniform(0.15, 0.5) * width
            alive[k] = np.abs(u - c) > L / 2 if rng.random() < 0.5 else (u - c) * rng.choice([-1, 1]) < 0
    alive[0] = True
    # effective beds: count alive boundaries below; b0/b1 = nearest alive boundaries below/above
    bed = np.zeros((H, W), int)
    b0 = np.full((H, W), -1e9); b1 = np.full((H, W), 1e9); kk = np.zeros((H, W), int)
    for k in range(nbnd):
        below = (V >= Vb[k][None, :]) & alive[k][None, :]
        bed += below
        upd = below & (Vb[k][None, :] > b0)
        b0 = np.where(upd, Vb[k][None, :], b0); kk = np.where(upd, k, kk)
        above = (V < Vb[k][None, :]) & alive[k][None, :] & (Vb[k][None, :] < b1)
        b1 = np.where(above, Vb[k][None, :], b1)
    b0 = np.where(b0 < -1e8, V - 1, b0); b1 = np.where(b1 > 1e8, V + 1, b1)
    t = np.clip((V - b0) / np.maximum(b1 - b0, 1e-3), 0, 1)
    bed = kk
'''
s = s[:a] + new + s[b:]
# remove the old pinch block that came after (it recomputed alive) and old b0/b1/t lines
old_pinch_start = s.index("    # pinch-outs: some boundaries fade")
old_pinch_end = s.index("    cb = comp[np.clip(bed, 0, len(comp) - 1)]")
s = s[:old_pinch_start] + s[old_pinch_end:]
s = s.replace("""    return dict(h=h, ds=ds, rock=rock, bed=np.where(rock, bed, -1), block=np.where(rock, block, -1),
                unit=np.where(rock, unit, -1), comp=np.where(rock, cb, 0.0), bg=bg, ground=ground,
                alive=alive, Vb=Vb""", """    return dict(h=h, ds=ds, rock=rock, bed=np.where(rock, bed, -1), block=np.where(rock, block, -1),
                unit=np.where(rock, unit, -1), comp=np.where(rock, cb, 0.0), bg=bg, ground=ground,
                alive=alive, Vb=Vb""")
open('geology.py','w').write(s)
EOF
grep -n "pinch\|b0 =\|alive" geology.py | head -30

# output:
Traceback (most recent call last):
  File "<stdin>", line 28, in <module>
ValueError: substring not found
13:told to, pinch out; competent beds protrude; joints are few, perpendicular to bedding, and stop
66:    b0 = np.take_along_axis(Vb, np.clip(bed, 0, nbnd - 1), 0)
110:                # the lens also thins toward its ends (sheets pinch at the stack edge)
111:                pinch = (np.abs(uu) ** 6) * (0.5 if bi in (0, len(xs) - 2) else 0.0)
112:                keep = np.abs(t[m]) < 1 - pinch
145:                comp_frac=0.3, protrude=0.35, soft_recess=0.25, dip_deg=0.0, taper=0.35, pinch=0.25,
149:    some boundaries die out (pinch). Competent beds protrude with a rounded lip and a shadowed soffit.
171:    # pinch-outs: some boundaries fade (the two beds merge over a stretch)
172:    alive = np.ones((nbnd, W), bool)
174:        if rng.random() < pinch:
177:            alive[k] = np.abs(u - c) > L / 2 if rng.random() < 0.5 else (u - c) * rng.choice([-1, 1]) < 0
179:    b0 = np.take_along_axis(Vb, np.clip(bed, 0, nbnd - 1), 0)
243:                alive=alive, Vb=Vb, u=u, v=v, width=width, height=height + base_ground, seg=seg)

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s = open('geology.py').read()
a = s.index("    # bed index per texel\n    bed = np.zeros((H, W), int)\n    for k in range(1, nbnd):")
b = s.index("    cb = comp[np.clip(bed, 0, len(comp) - 1)]")
new = '''    # pinch-outs: some boundaries die out along strike (the two beds merge over a stretch)
    alive = np.ones((nbnd, W), bool)
    for k in range(1, nbnd):
        if rng.random() < pinch:
            c = rng.uniform(0, width)
            L = rng.uniform(0.15, 0.5) * width
            alive[k] = np.abs(u - c) > L / 2 if rng.random() < 0.5 else (u - c) * rng.choice([-1, 1]) < 0
    # effective beds: nearest alive boundaries below (b0, index kk) and above (b1)
    b0 = np.full((H, W), -1e9); b1 = np.full((H, W), 1e9); bed = np.zeros((H, W), int)
    for k in range(nbnd):
        Vk = Vb[k][None, :]; ak = alive[k][None, :]
        upd = (V >= Vk) & ak & (Vk > b0)
        b0 = np.where(upd, Vk, b0); bed = np.where(upd, k, bed)
        dn = (V < Vk) & ak & (Vk < b1)
        b1 = np.where(dn, Vk, b1)
    b0 = np.where(b0 < -1e8, V - 1, b0); b1 = np.where(b1 > 1e8, V + 1, b1)
    t = np.clip((V - b0) / np.maximum(b1 - b0, 1e-3), 0, 1)
'''
s = s[:a] + new + s[b:]
# joints: master (through-going) joints persist over runs of beds
a = s.index("    # joints: candidate positions per bed")
b = s.index("    # open joints are slots")
new = '''    # joints: perpendicular to bedding (vertical here). Most stop at bed boundaries (T-junctions);
    # a few 'master' joints run through several consecutive beds (the ones Genge would draw).
    block = np.zeros((H, W), int)
    jmask = np.zeros((H, W), bool)
    nb = 0
    masters = []
    for _ in range(int(width / (joint_mean * 3)) + 1):
        k0 = rng.integers(0, nbnd); run = int(rng.integers(2, 8))
        masters.append((rng.uniform(0, width), k0, k0 + run))
    for k in range(nbnd):
        rows = bed == k
        if not rows.any():
            continue
        xs = np.cumsum(lognormal(rng, joint_mean, joint_sigma, int(width / joint_mean * 3) + 4)) - rng.uniform(0, joint_mean)
        xs = xs[(xs > 0) & (xs < width)]
        xs = xs[rng.random(len(xs)) < joint_keep * 2]
        xs = np.sort(np.concatenate([xs, [x for x, a0, a1 in masters if a0 <= k < a1]]))
        edges = np.concatenate([[0], xs, [width]])
        idx = np.searchsorted(edges, u) - 1
        block[rows] = (nb + idx[None, :].repeat(H, 0))[rows]
        nb += len(edges)
        for x in xs:
            jmask[rows & (np.abs(U - x) < joint_gap * 0.5 + ds * 0.5)] = True
'''
s = s[:a] + new + s[b:]
s = s.replace("""                alive=alive, Vb=Vb, u=u, v=v, width=width, height=height + base_ground, seg=seg)""",
"""                alive=alive, Vb=Vb, u=u, v=v, width=width, height=height + base_ground, seg=seg,
                alive_tex=np.ones((H, W), bool))""")
open('geology.py','w').write(s)
EOF
python3 -c "import sys; sys.path.insert(0,'.'); import geology; g=geology.bedded_face(); print(g['h'].shape, g['bed'].max())"

# output:
(315, 400) 57

cd ~/work/pixelart-studies/books/genge; mkdir -p stages/e3; cat > src/e3_beds.py <<'EOF'
"""E3. Bedding truth on a bedded cliff (West Bay, Bridport Sands; Genge Fig 6.6-6.7).
Rows: scales. Columns: photo | Genge 6.7a | schematic (Fig 4.7 mistake) | joints emphasized (Fig 4.6
mistake) | photo light (no rules) | Genge rules."""
import sys, numpy as np
sys.path.insert(0, 'src')
import geology, render, genge, palette
from util import load, area_down, up, sheet

seed = int(sys.argv[1]) if len(sys.argv) > 1 else 2
TRUE = dict(width=40, height=34, ds=0.05, seed=seed, bed_mean=1.3, sigma=0.55, thin_up=0.65, comp_frac=0.3,
            protrude=0.45, soft_recess=0.3, taper=0.45, pinch=0.3, joint_mean=2.5, furrows=0.35)
SCHEM = dict(TRUE, sigma=0.02, thin_up=0.0, taper=0.0, pinch=0.0, comp_frac=0.5)
L = (-0.4, -0.7, 0.6)
ramp = [palette.ramp('bridport_sand')]
gr = palette.ramp('beach_sand')
ph_close = load('refs/West_Bay_Bridport_cliffs/commons_West_Bay_The_beach_and_Bridport_Sands_cliffs_geo_836a9cc4.jpg')
ph_far = load('refs/West_Bay_Bridport_cliffs/commons_EastCliff_WestBay_jpg_c7f93602.jpg')
fig = load('plates/crops/fig6_7a_westbay.png')

variants = [('schematic', SCHEM, dict(bed_run=(50, 0.01), joint_ratio=1e9, joint_min_px=1)),
            ('joints emph.', TRUE, dict(joint_ratio=1e9, joint_min_px=1)),
            ('photo light', TRUE, dict(shadows='photo', ao_inset=0, ao_deep=2, soffit=False, lip=False,
                                       bed_lines=False, joints=False, contact=False, streak=0.2)),
            ('Genge rules', TRUE, dict())]
scenes = {}
rows = []
for lab, mpp, crop in [('close 0.1', 0.1, (14, 26, 12, 9)), ('landscape 0.35', 0.35, None), ('region 1.0', 1.0, None)]:
    r = []
    H0 = ph_close.shape[0]
    if crop:
        pc = ph_close[int(H0*0.15):int(H0*0.62), int(ph_close.shape[1]*0.35):int(ph_close.shape[1]*0.75)]
        r.append(('photo', up(area_down(pc, 120), 3)))
        r.append(('Genge 6.7a', up(area_down(fig[60:420, 80:560], 120), 3)))
    else:
        pf = ph_far[int(ph_far.shape[0]*0.18):int(ph_far.shape[0]*0.62), int(ph_far.shape[1]*0.05):int(ph_far.shape[1]*0.98)]
        w = int(round(40 / mpp * 3)) if mpp < 0.5 else int(round(120 / mpp))
        k = max(1, int(round(360 / (pf.shape[1] * (w / pf.shape[1])))))
        pd = area_down(pf, w)
        kk = max(1, int(round(360 / pd.shape[1])))
        r.append(('photo (whole cliff)', up(pd, kk)))
        fd = area_down(fig[60:700, 0:640], int(40 / mpp))
        r.append(('Genge 6.7a', up(fd, max(1, int(round(360 / fd.shape[1]))))))
    for name, gp, pp in variants:
        key = id(gp)
        if key not in scenes:
            geo = geology.bedded_face(**gp)
            scenes[key] = (geo, render.Scene(geo, render.face_frame(geo), L, ao_radius=1.5))
        geo, sc = scenes[key]
        cam = render.front_camera(geo, mpp)
        if crop:
            # a 12 x 9 m window at close range
            x0, z0, w, hgt = crop
            import engine as E
            ctr = np.array([x0 + w / 2, 0, z0 + hgt / 2])
            cam = E.Camera(ctr - np.array([0, 1, 0]) * 200, ctr, int(w / mpp), int(hgt / mpp), ortho_w=w)
        g = sc.gbuffer(cam, mpp, ss=3)
        prm = dict(pp); prm['mpp'] = mpp
        step, lit, lines, p = genge.paint(g, sc.L_local, prm)
        gs = genge.ground_steps(g, sc.L_local, p)
        sky = palette.sky_rgb(cam.H, cam.W)
        img = genge.colorize(step, g, ramp, gr, sky, gs)
        k = max(1, int(round(360 / cam.W)))
        r.append((f'{name}', up(img, k)))
    rows.append(r)
sheet(rows, title='E3  bedded cliff (West Bay): photo | Genge 6.7a | schematic | joints emphasized | photo light | Genge rules').save('stages/e3/e3_beds.png')
EOF
nice -n 10 python3 src/e3_beds.py 2

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s=open('geology.py').read()
s=s.replace("""def beds_sequence(height, rng, bed_mean=0.6, sigma=0.5, thin_up=0.0, comp_frac=0.35, packets=None):""",
"""def beds_sequence(height, rng, bed_mean=0.6, sigma=0.5, thin_up=0.0, comp_frac=0.35, packets=None, comp_mean=None):""")
s=s.replace("""        t = float(lognormal(rng, m, sigma, 1)[0])
        tops.append(tops[-1] + t)
        comp.append(1.0 if rng.random() < comp_frac else 0.0)""","""        c = 1.0 if rng.random() < comp_frac else 0.0
        mm = (comp_mean * (1 - 0.5 * thin_up * f)) if (c and comp_mean) else m
        t = float(lognormal(rng, mm, sigma, 1)[0])
        tops.append(tops[-1] + t)
        comp.append(c)""")
s=s.replace("""                bluff_setback=4.0, furrows=0.15, skyline_rough=1.0, units=None, base_ground=1.5):""",
"""                bluff_setback=4.0, furrows=0.15, skyline_rough=1.0, units=None, base_ground=1.5,
                comp_mean=None, furrow_depth=0.6, furrow_spacing=1.5, buttress=0.0, buttress_spacing=10.0):""")
s=s.replace("""    tops, comp = beds_sequence(height + 2, rng, bed_mean, sigma, thin_up, comp_frac)""",
"""    tops, comp = beds_sequence(height + 2, rng, bed_mean, sigma, thin_up, comp_frac, comp_mean=comp_mean)""")
s=s.replace("""    if furrows > 0:
        fz = smooth_noise(W, 0.6 / ds, rng, 1.0)
        fz = np.clip(fz, 0, None) ** 2
        h = h - furrows * 0.25 * fz[None, :] * (0.6 + 0.4 * noise2((H, W), 6, rng))""",
"""    if furrows > 0:
        # grooves: rounded, spacing ~furrow_spacing, each dying out at its own height (they start
        # below cemented ledges and run down the soft sand)
        fz = smooth_noise(W, furrow_spacing / ds * 0.35, rng, 1.0)
        fz = np.clip(fz, 0, None) ** 1.5
        reach = height * (0.3 + 0.7 * np.clip(0.5 + 0.5 * smooth_noise(W, furrow_spacing / ds, rng, 1.0), 0, 1))
        fade = np.clip((reach[None, :] - V) / 3.0, 0, 1)
        h = h - furrows * furrow_depth * fz[None, :] * fade * np.where(cb > 0.5, 0.5, 1.0)
    if buttress > 0:
        bz = smooth_noise(W, buttress_spacing / ds * 0.3, rng, 1.0)
        h = h + buttress * bz[None, :]""")
open('geology.py','w').write(s)
s=open('e3_beds.py').read()
s=s.replace("""TRUE = dict(width=40, height=34, ds=0.05, seed=seed, bed_mean=1.3, sigma=0.55, thin_up=0.65, comp_frac=0.3,
            protrude=0.45, soft_recess=0.3, taper=0.45, pinch=0.3, joint_mean=2.5, furrows=0.35)""",
"""TRUE = dict(width=40, height=34, ds=0.05, seed=seed, bed_mean=1.4, sigma=0.5, thin_up=0.6, comp_frac=0.4,
            comp_mean=0.3, protrude=0.35, soft_recess=0.06, taper=0.45, pinch=0.35, joint_mean=6.0,
            joint_keep=0.15, furrows=1.0, furrow_depth=0.7, furrow_spacing=2.0, buttress=1.2, buttress_spacing=12)""")
open('e3_beds.py','w').write(s)
EOF
cd ..; nice -n 10 python3 src/e3_beds.py 2

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge/src; cat >> geology.py <<'EOF'


# ---------------------------------------------------------------- basalt canyon (Yakima; Genge Ch. 7, p. 229)
def canyon(size=(600.0, 400.0), ds=1.0, seed=4, depth=220.0, flow_mean=28.0, flow_sigma=0.35,
           rim_frac=0.35, rim_steep=0.8, gully_spacing=45.0, gully_depth=25.0, river_w=30.0):
    """A river canyon cut into near-horizontal lava flows. Each flow has a competent lower part
    (colonnade: forms a steep rimrock riser, rim_frac of the flow) and a weak rubbly top (a bench).
    Beds are horizontal, so every flow boundary follows the contours and Vs up each gully (Genge
    p. 229: 'Only a horizontal boundary will run parallel to contours'). Frame: E1 east, E2 north, N up.
    The river runs along E1 at y = size[1]*0.3; the far wall rises to the north."""
    rng = _rng(seed)
    Wm, Hm = size
    W, H = int(Wm / ds), int(Hm / ds)
    x = (np.arange(W) + 0.5) * ds
    y = (np.arange(H)[::-1] + 0.5) * ds        # row 0 = north (far)
    X, Y = np.meshgrid(x, y)
    yr = Hm * 0.3 + 12 * smooth_noise(W, 60 / ds, rng, 1.0)[None, :]
    d = np.abs(Y - yr) - river_w / 2
    # smooth canyon profile: convex-concave wall, plateau at top
    prof = depth * np.clip(d / (Hm * 0.55), 0, 1) ** 0.9
    # gullies: incisions spaced along the wall, running straight down it
    gx = smooth_noise(W, gully_spacing / ds * 0.25, rng, 1.0)
    g = np.clip(gx, 0, None) ** 1.3
    wall = np.clip(d / (Hm * 0.55), 0, 1)
    Z = prof - gully_depth * g[None, :] * np.sin(np.pi * np.clip(wall * 1.15, 0, 1)) ** 0.7
    Z += 3 * noise2((H, W), 6 / ds, rng)
    Z = np.maximum(Z, 0)
    # flows: horizontal boundaries from the river up, thickness lognormal
    tops = [-5.0]
    while tops[-1] < depth + 20:
        tops.append(tops[-1] + float(lognormal(rng, flow_mean, flow_sigma, 1)[0]))
    tops = np.array(tops)
    k = np.clip(np.searchsorted(tops, Z) - 1, 0, len(tops) - 2)
    z0, z1 = tops[k], tops[k + 1]
    t = (Z - z0) / (z1 - z0)
    # terrace transform: the competent lower rim_frac of each flow becomes a steep riser, the rest a bench
    rf = rim_frac * (0.8 + 0.4 * rng.random(len(tops)))[k]
    tt = np.where(t < rf, (t / rf) * (1 - rim_steep * 0 ) , 1.0)
    # remap elevations: riser rises the whole flow's height quickly; bench nearly flat
    riser = np.clip(t / rf, 0, 1)
    bench = np.clip((t - rf) / (1 - rf), 0, 1)
    Zt = z0 + (z1 - z0) * (0.85 * riser + 0.15 * bench)
    Zt = np.where(wall > 0.02, Zt, Z)
    Zs = ndi.gaussian_filter(Zt, 0.6)
    comp = (t < rf) & (wall > 0.05)
    # exposure vs cover: steep ground is bare rock (rimrock), gentle ground grassed, scree below risers
    gy, gxx = np.gradient(Zs, ds)
    slope = np.degrees(np.arctan(np.hypot(gxx, gy)))
    rock = (slope > 38) & (wall > 0.03)
    rock = ndi.binary_opening(rock, np.ones((2, 2)))
    scree = (slope > 26) & ~rock & (wall > 0.03)
    water = d < 0
    unit = np.where(rock, 0, np.where(scree, 1, 2))
    unit = np.where(water, 3, unit)
    return dict(h=Zs, ds=ds, rock=rock, bed=np.where(rock, k, -1), block=np.where(rock, k, -1),
                unit=unit, comp=comp.astype(float), bg=np.zeros((H, W), bool), ground=~rock,
                flow=k, flow_t=t, slope=slope, water=water, scree=scree, x=x, y=y, width=Wm, height=Hm,
                tops=tops)
EOF
echo ok

# output:
ok

"""landscape.py — Genge's landscape rules (Ch. 7, Ch. 15, p. 229) as a terrain painter.

Rules, each a switch:
  horizons   local horizons (a near ridge in front of farther ground) get a clean edge: the far
             ground above the edge is hazed one step lighter, the near ridge-top pixel keeps its
             value (no dither) — Genge: 'horizons first'; multiple horizons 'divide the view into
             regions of different scale and detail' (§7.1)
  rimrock    competent flow bases (colonnade risers) are the 'heavier lines' of protruding beds:
             lit lip on top, darkest pocket at the foot where the riser meets the bench/scree (p. 106)
  traces     flow boundaries under grass: impersistent 1-step lines ('faint bedding trace on
             landscape … short parallel discontinuous lines', p. 106). Because flows are horizontal
             they follow contours and V up the gullies by construction (p. 229)
  gullies    valley axes drawn as impersistent darker runs ('minor gulleys provide volume', Fig 7.2b)
  scree      fall-line streaks on scree ('slope lines on scree slope', Fig 8.8b)
  haze       value range compressed toward the sky with distance (the aerial study's rule)
"""
import numpy as np
from scipy import ndimage as ndi
import genge

DEF = dict(horizons=True, rimrock=True, traces=True, gullies=True, scree=True, haze=0.5,
           shadows='mass', mass_frac=0.02, ao_inset=1.0, ao_lo=0.12, ao_hi=0.4, ao_deep=2.0,
           trace_run=5, gully_curv=0.9, horizon_jump=0.12, seed=0, lod_min_px=2.0)


def paint(g, L, sky_rgb, ramps, params=None):
    p = dict(DEF)
    if params:
        p.update(params)
    H, W = g['hit'].shape
    N = g['N']
    L = np.asarray(L, float)
    ndl = (N * L[None, None, :]).sum(-1)
    hit = g['hit']
    cast = g['shadow'] > 0.5
    self_sh = ndl <= 0.03
    if p['shadows'] == 'photo':
        sh = self_sh | cast
    elif p['shadows'] == 'genge':
        sh = self_sh
    else:
        lab, n = ndi.label(cast & hit)
        keep = np.zeros_like(cast)
        if n:
            sizes = ndi.sum(np.ones_like(lab), lab, np.arange(1, n + 1))
            keep = np.isin(lab, np.flatnonzero(sizes >= p['mass_frac'] * hit.sum()) + 1)
        sh = self_sh | keep
    lit = hit & ~sh
    sky = 0.5 + 0.5 * N[..., 2]
    v = np.where(lit, 4.0 + 3.0 * np.clip(ndl, 0, 1) ** 0.8, 1.1 + 2.0 * sky)
    t = np.clip((g['ao'] - p['ao_lo']) / (p['ao_hi'] - p['ao_lo']), 0, 1)
    v = np.maximum(v - p['ao_inset'] * t, np.where(lit, 4.0, 0.6))
    unit = g['unit']
    rng = np.random.default_rng(p['seed'])
    if p['scree']:
        s = genge.streaks((H, W), 0.0, 1, p['seed'] + 3)
        vs = ndi.uniform_filter1d(rng.standard_normal((H, W)), 5, axis=0)
        vs = vs / (vs.std() + 1e-9)
        v = np.where(unit == 1, v + 0.45 * vs, v)
    # grass: soft clumped texture
    gtex = ndi.gaussian_filter(rng.standard_normal((H, W)), 0.7)
    v = np.where(unit == 2, v + 0.35 * gtex / (gtex.std() + 1e-9), v)
    v = np.where(lit, np.clip(v, 4, 7.4), np.clip(v, 0.6, 3.49))
    B = genge._bayer(H, W)
    frac = v - np.floor(v)
    seam = np.abs(frac - 0.5) < 0.2
    step = np.where(seam, np.floor(v - 0.5 + B).astype(int), np.round(v).astype(int))
    step = np.where(lit, np.clip(step, 4, 7), np.clip(step, 0, 3))
    flow = g['bed_all']
    comp = g['comp'] > 0.5
    rockm = g['rock']
    lines = {}
    below = np.roll(flow, -1, 0)
    fb = (flow != below) & hit & np.roll(hit, -1, 0)
    # flow thickness in px: skip line classes whose spacing is under lod_min_px (Genge p. 98: in the
    # distance only a proportion of the layers can be drawn)
    fpx = g.get('flow_px', np.full((H, W), 99.0))
    big = fpx >= p['lod_min_px'] * 2
    if p['rimrock']:
        # foot of a riser: the pixel just below a competent-rock pixel that is not competent rock
        foot = np.roll(comp & rockm, 1, 0) & ~(comp & rockm) & hit & big
        step = np.where(foot, np.minimum(step, 1), step)
        lipm = comp & rockm & ~np.roll(comp & rockm, 1, 0) & lit & big
        step = np.where(lipm, np.minimum(step + 1, 7), step)
        lines['foot'], lines['lip'] = foot, lipm
    if p['traces']:
        tr = fb & ~rockm & big
        tr = genge.impersistent(tr, p['trace_run'], 0.6, p['seed'])
        step = np.where(tr, np.maximum(step - 1, np.where(lit, 4, 0)), step)
        lines['trace'] = tr
    if p['gullies'] and 'curv' in g:
        gm = (g['curv'] > p['gully_curv']) & hit & ~rockm
        gm = gm & (ndi.uniform_filter(gm.astype(float), 3) > 0.25)
        try:
            from skimage.morphology import skeletonize
            gm = skeletonize(gm)
        except Exception:
            pass
        gm = genge.impersistent(gm, 6, 0.4, p['seed'] + 1)
        step = np.where(gm, np.maximum(step - 1, 0), step)
        lines['gully'] = gm
    # colour
    out = np.array(sky_rgb, np.uint8).copy()
    for k, r in enumerate(ramps):
        m = (unit == k) & hit
        out[m] = r[np.clip(step[m], 0, 7)]
    # haze by depth (value range compressed toward the sky)
    if p['haze'] > 0:
        d = g['depth']
        dn = np.clip((d - np.percentile(d[hit], 5)) / (np.percentile(d[hit], 99) - np.percentile(d[hit], 5) + 1e-6), 0, 1)
        a = (p['haze'] * dn ** 1.3)[..., None]
        skyc = np.asarray(sky_rgb, float)
        hz = out * (1 - a) + skyc * a
        # quantize the haze into 4 bins so it stays palette-like
        out = np.where(hit[..., None], hz, out).astype(np.uint8)
    if p['horizons']:
        d = g['depth']
        up = np.roll(d, 1, 0)
        edge = hit & np.roll(hit, 1, 0) & (up > d * (1 + p['horizon_jump']))
        # far ground just above the edge gets lighter (separation), near edge pixel darkened 1 step
        far = np.roll(edge, -1, 0)
        out = np.where(far[..., None], np.clip(out.astype(int) + 14, 0, 255), out).astype(np.uint8)
        lines['horizon'] = edge
    return out, step, lines


def curvature(h, ds, sigma=2.0):
    hs = ndi.gaussian_filter(h, sigma)
    return ndi.laplace(hs) / ds ** 2 * sigma ** 2

cd ~/work/pixelart-studies/books/genge/src; python3 - <<'EOF'
s=open('render.py').read()
s=s.replace("""        out['L_local'] = self.L_local""","""        for k in geo.get('extra', []):
            out[k] = g.take(geo[k], -1 if np.issubdtype(np.asarray(geo[k]).dtype, np.integer) else 0.0)
        out['L_local'] = self.L_local""")
open('render.py','w').write(s)
EOF
cd ..; mkdir -p stages/e6; cat > src/e6_canyon.py <<'EOF'
"""E6. Landscape rules on a Yakima-style basalt canyon: photo light vs Genge landscape rules,
landscape perspective + 30-degree iso."""
import sys, numpy as np
sys.path.insert(0, 'src')
import geology, render, landscape, palette, engine as E
from util import load, area_down, up, sheet

geo = geology.canyon(seed=int(sys.argv[1]) if len(sys.argv) > 1 else 4)
geo['curv'] = landscape.curvature(geo['h'], geo['ds'], 3.0)
geo['bed_all'] = geo['flow']
geo['flow_m'] = np.diff(geo['tops'])[geo['flow']]
geo['extra'] = ['curv', 'bed_all', 'flow_m']
fr = E.Frame(O=(0, 0, 0), E1=(1, 0, 0), E2=(0, 1, 0), N=(0, 0, 1))
L = (-0.5, -0.35, 0.55)
sc = render.Scene(geo, fr, L, ao_radius=25.0)
ramps = [palette.ramp('basalt'), palette.ramp('scree'), palette.ramp('grass_dry'), palette.ramp('greenstone', (-10, 5, -60))]
views = {
 'landscape': E.Camera((300, -90, 160), (300, 300, 90), 320, 180, fov_deg=55),
 'iso 30deg': E.Camera(np.array((300, 200, 60)) + 800 * np.array([-0.61, -0.61, 0.5]), (300, 200, 60), 320, 200, ortho_w=380),
}
ph = load('refs/Yakima_River_Canyon_basalt_cliffs/openverse_Yakima_River_Canyon_Scenic_Byway_b438417e.jpg')
rows = []
for name, cam in views.items():
    g = sc.gbuffer(cam, 1.0, ss=2)
    if cam.ortho_w:
        mpp = cam.ortho_w / cam.W
        pixm = np.full(g['depth'].shape, mpp)
    else:
        pixm = g['depth'] * 2 * np.tan(np.radians(cam.fov) / 2) / cam.W
    g['flow_px'] = g['flow_m'] / np.maximum(pixm, 1e-3)
    sky = palette.sky_rgb(cam.H, cam.W)
    r = [('photo', up(area_down(ph, 320), 2))] if name == 'landscape' else [('', np.full((10, 10, 3), 245, np.uint8))]
    for lab, prm in [('photo light', dict(shadows='photo', horizons=False, rimrock=False, traces=False, gullies=False, ao_inset=0)),
                     ('Genge landscape rules', dict())]:
        img, step, lines = landscape.paint(g, sc.L_local, sky, ramps, prm)
        r.append((f'{lab} | {name}', up(img, 2)))
    rows.append(r)
sheet(rows, title='E6  basalt canyon: horizontal flows, rimrock, traces, gullies, local horizons').save('stages/e6/e6_canyon.png')
EOF
nice -n 10 python3 src/e6_canyon.py

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/geology.py').read()
s=s.replace("Zt = z0 + (z1 - z0) * (0.85 * riser + 0.15 * bench)","Zt = z0 + (z1 - z0) * (rim_steep * riser + (1 - rim_steep) * bench)")
s=s.replace("rim_frac=0.35, rim_steep=0.8,","rim_frac=0.18, rim_steep=0.4,")
open('src/geology.py','w').write(s)
s=open('src/e6_canyon.py').read()
s=s.replace("geo = geology.canyon(seed=int(sys.argv[1]) if len(sys.argv) > 1 else 4)",
 "geo = geology.canyon(size=(1400.0, 1000.0), ds=2.0, depth=320.0, gully_spacing=90, gully_depth=35, river_w=45, seed=int(sys.argv[1]) if len(sys.argv) > 1 else 4)")
s=s.replace("sc = render.Scene(geo, fr, L, ao_radius=25.0)","sc = render.Scene(geo, fr, L, ao_radius=30.0)")
s=s.replace("""'landscape': E.Camera((300, -90, 160), (300, 300, 90), 320, 180, fov_deg=55),
 'iso 30deg': E.Camera(np.array((300, 200, 60)) + 800 * np.array([-0.61, -0.61, 0.5]), (300, 200, 60), 320, 200, ortho_w=380),""",
"""'landscape': E.Camera((620, 60, 330), (760, 700, 150), 320, 180, fov_deg=60),
 'iso 30deg': E.Camera(np.array((700, 520, 150)) + 2000 * np.array([-0.61, -0.61, 0.5]), (700, 520, 150), 320, 200, ortho_w=640),""")
s=s.replace("('Genge landscape rules', dict())","('Genge landscape rules', dict(haze=0.3))")
s=s.replace("gullies=False, ao_inset=0)","gullies=False, ao_inset=0, haze=0.3)")
open('src/e6_canyon.py','w').write(s)
EOF
nice -n 10 python3 src/e6_canyon.py

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/geology.py').read()
s=s.replace("    Z += 3 * noise2((H, W), 6 / ds, rng)\n","    Z += 1.5 * noise2((H, W), 20 / ds, rng)\n")
s=s.replace("""    gy, gxx = np.gradient(Zs, ds)
    slope = np.degrees(np.arctan(np.hypot(gxx, gy)))
    rock = (slope > 38) & (wall > 0.03)
    rock = ndi.binary_opening(rock, np.ones((2, 2)))
    scree = (slope > 26) & ~rock & (wall > 0.03)""","""    gy, gxx = np.gradient(Zs, ds)
    slope = np.degrees(np.arctan(np.hypot(gxx, gy)))
    # exposure as masses (Genge §7.3: outline areas of exposure first): rimrock where the riser is
    # steep AND the flow's competent base, grown along strike; scree fans below rock; grass elsewhere
    ssm = ndi.gaussian_filter(slope, 1.0)
    rock = (ssm > 36) & comp & (wall > 0.03)
    rock = ndi.binary_closing(rock, np.ones((1, 5)))
    rock = ndi.binary_opening(rock, np.ones((2, 3)))
    lab, n = ndi.label(rock)
    if n:
        sz = ndi.sum(np.ones_like(lab), lab, np.arange(1, n + 1))
        rock = np.isin(lab, np.flatnonzero(sz >= 25) + 1)
    below = np.zeros_like(rock)
    for sft in range(1, 6):
        below |= np.roll(rock, sft, 0) | np.roll(rock, -sft, 0)
    scree = below & ~rock & (ssm > 24) & (wall > 0.03)
    scree = ndi.binary_opening(scree, np.ones((2, 2)))""")
open('src/geology.py','w').write(s)
s=open('src/e6_canyon.py').read()
s=s.replace("'landscape': E.Camera((620, 60, 330), (760, 700, 150), 320, 180, fov_deg=60),",
            "'landscape': E.Camera((200, 110, 190), (1000, 470, 70), 320, 180, fov_deg=60),")
open('src/e6_canyon.py','w').write(s)
s=open('src/landscape.py').read()
s=s.replace("""        tr = fb & ~rockm & big""","""        # only the prominent flows (thickest fraction), as Genge draws only the most prominent traces
        thick = g.get('flow_m', np.zeros((H, W)))
        thr = np.percentile(thick[hit], 100 * (1 - p.get('trace_frac', 0.4))) if hit.any() else 0
        tr = fb & ~rockm & big & (thick >= thr)""")
open('src/landscape.py','w').write(s)
EOF
nice -n 10 python3 src/e6_canyon.py

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/geology.py').read()
s=s.replace("    rock = (ssm > 36) & comp & (wall > 0.03)\n",
"""    exposed_flow = rng.random(len(tops)) < 0.45
    rock = (ssm > 36) & comp & (wall > 0.12) & exposed_flow[k]
""")
open('src/geology.py','w').write(s)
s=open('src/e6_canyon.py').read()
s=s.replace("gully_spacing=90, gully_depth=35","gully_spacing=70, gully_depth=60")
open('src/e6_canyon.py','w').write(s)
EOF
nice -n 10 python3 src/e6_canyon.py

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; sed -i '' 's/^L = (-0.5, -0.35, 0.55)/L = (0.75, -0.35, 0.45)/' src/e6_canyon.py && nice -n 10 python3 src/e6_canyon.py

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; sed -i '' "s/'grass_dry':     dict(lit=(66, 30, 88), sh=(34, 16, 110)),/'grass_dry':     dict(lit=(66, 30, 85), sh=(32, 14, 62)),/" src/palette.py; grep grass_dry src/palette.py; nice -n 10 python3 src/e6_canyon.py; mkdir -p stages/e2; cat > src/e2_light.py <<'EOF'
"""E2. Cast shadow vs ambient occlusion (Genge §2.9.3) under four lights, on the tor and the bedded
cliff. Columns: photo light (all cast shadows) | genge (none) | mass (big masses only) ."""
import sys, numpy as np
sys.path.insert(0, 'src')
import geology, render, genge, palette, engine as E
from util import up, sheet

LIGHTS = [('front-left high', (-0.55, -0.55, 0.63)), ('raking right low', (0.9, -0.25, 0.3)),
          ('back-lit (sun behind)', (0.3, 0.8, 0.5)), ('overcast', None)]
subjects = [('tor 0.1 m/px', geology.tor(seed=5, sheet_mean=0.75, joint_mean=2.2, recess=0.5, stacks=3, height=7.0),
             0.1, [palette.ramp('tor_granite')], palette.ramp('grass_dry')),
            ('bedded cliff 0.35 m/px', geology.bedded_face(width=40, height=34, ds=0.05, seed=2, bed_mean=1.4, sigma=0.5,
             thin_up=0.6, comp_frac=0.4, comp_mean=0.3, protrude=0.35, soft_recess=0.06, taper=0.45, pinch=0.35,
             joint_mean=6.0, joint_keep=0.15, furrows=1.0, furrow_depth=0.7, furrow_spacing=2.0, buttress=1.2,
             buttress_spacing=12), 0.35, [palette.ramp('bridport_sand')], palette.ramp('beach_sand'))]
rows = []
for sname, geo, mpp, ramps, gr in subjects:
    for lname, L in LIGHTS:
        ov = L is None
        Lw = (0.0, -0.3, 1.0) if ov else L
        sc = render.Scene(geo, render.face_frame(geo), Lw, ao_radius=1.5)
        cam = render.front_camera(geo, mpp)
        g = sc.gbuffer(cam, mpp, ss=3)
        if ov:   # overcast: no sun at all; the value is sky view (up-facing) + occlusion only
            g['shadow'][:] = 1.0
            g['N'] = g['N'].copy()
        r = []
        for mode in ['photo', 'genge', 'mass']:
            prm = dict(shadows=mode, mpp=mpp)
            if mode == 'photo':
                prm.update(ao_inset=0, ao_deep=2, soffit=False, lip=False, bed_lines=False, joints=False, contact=False)
            if ov:
                prm['shadows'] = 'photo'
            step, lit, lines, p = genge.paint(g, sc.L_local if not ov else np.array([0, 0, -1.0]), prm)
            gs = genge.ground_steps(g, sc.L_local, p)
            img = genge.colorize(step, g, [palette.ramp(list(palette.MATS)[0]) if False else ramps[0]], gr,
                                 palette.sky_rgb(cam.H, cam.W, 'overcast' if ov else 'day'), gs)
            k = max(1, int(round(300 / cam.W)))
            lab = {'photo': 'photo light', 'genge': 'Genge: AO only', 'mass': 'pixel: mass shadows + AO'}[mode]
            if ov:
                lab = {'photo': 'overcast, no rules', 'genge': 'overcast + rules', 'mass': 'overcast + rules'}[mode]
            r.append((f'{sname} | {lname} | {lab}', up(img, k)))
        rows.append(r)
sheet(rows, title='E2  cast shadow vs ambient occlusion, four lights', fs=12).save('stages/e2/e2_light.png')
EOF
nice -n 10 python3 src/e2_light.py

# output:
'grass_dry':     dict(lit=(66, 30, 85), sh=(32, 14, 62)),

cd ~/work/pixelart-studies/books/genge; mkdir -p stages/tests; ls refs/Sucia_Island_sandstone_honeycomb/; cat > src/t_subjects.py <<'EOF'
"""Tests on the real subjects: Ruby Beach sea stack (Olympic coast) and Sucia sandstone (San Juans).
photo | before (photo light, no rules) | after (Genge rules), at landscape and close."""
import sys, numpy as np
sys.path.insert(0, 'src')
import geology, render, genge, palette, engine as E
from util import load, area_down, up, sheet
from scipy import ndimage as ndi

def stack_mask(geo, seed=1):
    rng = np.random.default_rng(seed)
    u, v = geo['u'], geo['v']
    U, V = np.meshgrid(u, v)
    W = geo['width']
    # leaning tower: steep seaward (left) side, stepped landward side, ragged top (rock study's shape)
    c = W * 0.5 + (V / V.max()) * W * 0.08
    half = W * 0.36 * (1 - 0.55 * (V / V.max()) ** 1.4) + 1.2 * geology.smooth_noise(len(v), 20, rng, 1.0)[:, None]
    half = half + np.where(U > c, 1.5 * np.floor(V / 4.0) % 2, 0)
    top = V.max() - 3 - 4 * np.abs(geology.smooth_noise(len(u), 40, rng, 1.0))[None, :]
    return (np.abs(U - c) < half) & (V < top)

def run(name, geo, L, ramps, gr, photo, crops, extra=None, sky='day'):
    sc = render.Scene(geo, render.face_frame(geo), L, ao_radius=1.5)
    rows = []
    for lab, mpp in crops:
        cam = render.front_camera(geo, mpp)
        g = sc.gbuffer(cam, mpp, ss=3)
        k = max(1, int(round(300 / cam.W)))
        r = [('photo', up(area_down(photo, cam.W), k))]
        for nm, prm in [('before: photo light', dict(shadows='photo', ao_inset=0, ao_deep=2, soffit=False, lip=False,
                         bed_lines=False, joints=False, contact=False)), ('after: Genge rules', dict())]:
            prm = dict(prm, mpp=mpp, **(extra or {}))
            step, lit, lines, p = genge.paint(g, sc.L_local, prm)
            gs = genge.ground_steps(g, sc.L_local, p)
            img = genge.colorize(step, g, ramps, gr, palette.sky_rgb(cam.H, cam.W, sky), gs)
            r.append((f'{nm} | {lab}', up(img, k)))
        rows.append(r)
    sheet(rows, title=name).save(f'stages/tests/{name.split()[0]}.png')

# Ruby Beach stack: Hoh melange sandstone, beds dipping steeply, dark rock; coast ambient
geo = geology.bedded_face(width=30, height=34, ds=0.05, seed=7, bed_mean=1.1, sigma=0.6, thin_up=0.2, comp_frac=0.45,
                          comp_mean=0.8, protrude=0.6, soft_recess=0.35, dip_deg=-55, taper=0.5, pinch=0.3,
                          joint_mean=3.5, joint_keep=0.3, furrows=0.4, buttress=1.5, buttress_spacing=6, base_ground=1.0)
m = stack_mask(geo)
geo['bg'] = geo['bg'] | (~m & ~geo['ground'])
geo['rock'] = geo['rock'] & m
geo['h'] = np.where(geo['bg'], -40.0, geo['h'])
for kk in ('bed', 'block', 'unit'):
    geo[kk] = np.where(geo['rock'], geo[kk], -1)
ph = load('refs/Ruby_Beach_sea_stack/commons_Sea_stack_at_Ruby_Beach_jpg_e8a11f34.jpg')
H0, W0 = ph.shape[:2]
ph = ph[int(H0*0.08):int(H0*0.85), int(W0*0.18):int(W0*0.68)]
run('ruby_stack  Ruby Beach sea stack (Olympic coast)', geo, (-0.6, -0.5, 0.6), [palette.ramp('coast_rock')],
    palette.ramp('beach_sand'), ph, [('landscape 0.35', 0.35), ('close-ish 0.2', 0.2)], extra=dict(bed_angle=np.radians(-55)))

# Sucia sandstone: gently dipping tabular beds on a low shore bluff
geo = geology.bedded_face(width=24, height=8, ds=0.04, seed=11, bed_mean=0.35, sigma=0.5, thin_up=0.0, comp_frac=0.5,
                          comp_mean=0.4, protrude=0.25, soft_recess=0.12, dip_deg=12, taper=0.6, pinch=0.4,
                          joint_mean=2.0, joint_keep=0.3, furrows=0.3, base_ground=0.6)
ph = load('refs/Sucia_Island_sandstone_honeycomb/commons_Folded_Bedrock_Sucia_Island_and_Patos_Island_jpg_1f01218f.jpg')
run('sucia  Sucia Island sandstone (San Juans)', geo, (-0.5, -0.6, 0.6), [palette.ramp('sandstone')],
    palette.ramp('beach_sand'), ph, [('close 0.1', 0.1), ('landscape 0.35', 0.35)], extra=dict(bed_angle=np.radians(12)))
EOF
nice -n 10 python3 src/t_subjects.py

# output:
_sheet.jpg
commons_Folded_Bedrock_Sucia_Island_and_Patos_Island_jpg_1f01218f.jpg
commons_Sucia_in_Summer_153693_jpg_c244fb78.jpg
openverse_20100804_DSCF0398_jpg_e0b256e4.jpg
openverse_Layers_of_Texture_b25a30d6.jpg
openverse_Sea_Lions_with_Mt_Baker_Sucia_Island_Kayak_Trip_53645ff3.jpg
openverse_Still_Life_on_Sucia_Island_31b532cb.jpg
sources.tsv
~/work/pixelart-studies/books/genge/src/t_subjects.py:16: RuntimeWarning: invalid value encountered in power
  half = W * 0.36 * (1 - 0.55 * (V / V.max()) ** 1.4) + 1.2 * geology.smooth_noise(len(v), 20, rng, 1.0)[:, None]

cd ~/work/pixelart-studies/books/genge; python3 -c "
from PIL import Image
for f in ['openverse_Layers_of_Texture_b25a30d6.jpg','openverse_20100804_DSCF0398_jpg_e0b256e4.jpg','openverse_Still_Life_on_Sucia_Island_31b532cb.jpg']:
    im=Image.open('refs/Sucia_Island_sandstone_honeycomb/'+f); im.thumbnail((400,300)); im.save('(tmp))
" 2>/dev/null || python3 -c "
from PIL import Image
ims=[Image.open('refs/Sucia_Island_sandstone_honeycomb/'+f) for f in ['openverse_Layers_of_Texture_b25a30d6.jpg','openverse_20100804_DSCF0398_jpg_e0b256e4.jpg','openverse_Still_Life_on_Sucia_Island_31b532cb.jpg']]
s=Image.new('RGB',(1200,300))
for i,im in enumerate(ims): im.thumbnail((400,300)); s.paste(im,(i*400,0))
s.save('img/sucia_pick.png')"; ls img

# output:
contact_00.png
contact_01.png
contact_02.png
contact_03.png
contact_04.png
contact_05.png
contact_06.png
contact_07.png
contact_08.png
contact_09.png
contact_10.png
crops_check.png

cd ~/work/pixelart-studies/books/genge; python3 -c "
from PIL import Image
ims=[Image.open('refs/Sucia_Island_sandstone_honeycomb/'+f) for f in ['openverse_Layers_of_Texture_b25a30d6.jpg','openverse_20100804_DSCF0398_jpg_e0b256e4.jpg','openverse_Still_Life_on_Sucia_Island_31b532cb.jpg']]
s=Image.new('RGB',(1200,300))
for i,im in enumerate(ims): im.thumbnail((400,300)); s.paste(im,(i*400,0))
s.save('img/sucia_pick.png')"

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/t_subjects.py').read()
s=s.replace("half = half + np.where(U > c, 1.5 * np.floor(V / 4.0) % 2, 0)","half = half + np.where(U > c, 1.2 * np.floor(V / 5.0) - 0.25 * V, 0)")
s=s.replace("    return (np.abs(U - c) < half) & (V < top)","    m = (np.abs(U - c) < np.maximum(half, 0)) & (V < top)\n    return ndi.binary_opening(m, np.ones((15, 15)))")
s=s.replace("half = W * 0.36 * (1 - 0.55 * (V / V.max()) ** 1.4)","half = W * 0.36 * (1 - 0.55 * np.clip(V / V.max(), 0, 1) ** 1.4)")
s=s.replace("geo = geology.bedded_face(width=24, height=8, ds=0.04, seed=11, bed_mean=0.35, sigma=0.5, thin_up=0.0, comp_frac=0.5,\n                          comp_mean=0.4, protrude=0.25,",
 "geo = geology.bedded_face(width=24, height=8, ds=0.04, seed=11, bed_mean=1.0, sigma=0.5, thin_up=0.0, comp_frac=0.3,\n                          comp_mean=1.2, protrude=0.6,")
s=s.replace("commons_Folded_Bedrock_Sucia_Island_and_Patos_Island_jpg_1f01218f.jpg')","openverse_20100804_DSCF0398_jpg_e0b256e4.jpg')[30:560]")
open('src/t_subjects.py','w').write(s)
s=open('src/palette.py').read()
s=s.replace("'coast_rock':    dict(lit=(30, 5, 75), sh=(15, 4, 250)),","'coast_rock':    dict(lit=(40, 7, 80), sh=(17, 4, 250)),")
open('src/palette.py','w').write(s)
EOF
nice -n 10 python3 src/t_subjects.py 2>&1 | grep -v Warn

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/t_subjects.py').read()
s=s.replace("1.2 * geology.smooth_noise(len(v), 20, rng, 1.0)[:, None]","0.8 * geology.smooth_noise(len(v), 150, rng, 1.0)[:, None]")
s=s.replace("""        for nm, prm in [('before: photo light', dict(shadows='photo', ao_inset=0, ao_deep=2, soffit=False, lip=False,
                         bed_lines=False, joints=False, contact=False)), ('after: Genge rules', dict())]:""",
"""        vs = [('before: photo light', dict(shadows='photo', ao_inset=0, ao_deep=2, soffit=False, lip=False,
                         bed_lines=False, joints=False, contact=False)), ('after: Genge rules', dict())]
        if dark:
            vs.append(('after, dark-rock variant (lips, not lines)', dict(bed_lines=False, joint_ratio=0.1, ao_deep=0.7)))
        for nm, prm in vs:""")
s=s.replace("def run(name, geo, L, ramps, gr, photo, crops, extra=None, sky='day'):","def run(name, geo, L, ramps, gr, photo, crops, extra=None, sky='day', dark=False):")
s=s.replace("extra=dict(bed_angle=np.radians(-55)))","extra=dict(bed_angle=np.radians(-55)), dark=True)")
open('src/t_subjects.py','w').write(s)
EOF
nice -n 10 python3 src/t_subjects.py 2>&1 | grep -v Warn

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; cat > src/z8.py <<'EOF'
"""8x zoom: Rough Tor at 0.1 m/px — Genge Fig 2.9b | photo | Genge-rules render (same window)."""
import sys, numpy as np
sys.path.insert(0, 'src')
import geology, render, genge, palette
from util import load, area_down, up, sheet
geo = geology.tor(seed=5, sheet_mean=0.75, joint_mean=2.2, recess=0.5, stacks=3, height=7.0)
sc = render.Scene(geo, render.face_frame(geo), (-0.55, -0.55, 0.63), ao_radius=1.2)
cam = render.front_camera(geo, 0.1)
g = sc.gbuffer(cam, 0.1, ss=3)
out = {}
for nm, prm in [('photo light', dict(shadows='photo', ao_inset=0, ao_deep=2, soffit=False, lip=False, bed_lines=False,
                 joints=False, contact=False)), ('Genge rules', dict(mpp=0.1))]:
    step, lit, lines, p = genge.paint(g, sc.L_local, prm)
    out[nm] = genge.colorize(step, g, [palette.ramp('tor_granite')], palette.ramp('grass_dry'),
                            palette.sky_rgb(cam.H, cam.W), genge.ground_steps(g, sc.L_local, p))
ph = area_down(load('refs/Rough_Tor_Bodmin_Moor_granite_tor/commons_Roughtor_Bodmin_Moor_2_jpg_76d7aed8.jpg')[60:853, 20:1280], 90)
fig = area_down(load('plates/crops/fig2_9b_roughtor_ao.png')[60:, :], 120)
y0, x0 = 20, 30
rows = [[('Genge 2.9b (box-down)', up(fig[10:50, 30:80], 8)), ('photo (box-down)', up(ph[5:45, 20:70], 8))],
        [('photo light', up(out['photo light'][y0:y0+40, x0:x0+50], 8)), ('Genge rules', up(out['Genge rules'][y0:y0+40, x0:x0+50], 8))]]
sheet(rows, title='Rough Tor 0.1 m/px at 8x').save('sheets/Z_tor_8x.png')
EOF
nice -n 10 python3 src/z8.py

# output:
(Bash completed with no output)

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/genge.py').read()
s=s.replace("""    step = np.where(rock, step, -1)
    return step, lit, lines, p""","""    # despeckle: a lone pixel whose 4 neighbours all share one other step takes that step (no Morse dots)
    if p.get('despeckle', True):
        st = step.copy()
        nb = [np.roll(step, 1, 0), np.roll(step, -1, 0), np.roll(step, 1, 1), np.roll(step, -1, 1)]
        same = (nb[0] == nb[1]) & (nb[1] == nb[2]) & (nb[2] == nb[3]) & (nb[0] != step)
        keep_lines = np.zeros_like(rock)
        for m in lines.values():
            keep_lines |= m
        st = np.where(same & rock & ~keep_lines & np.roll(rock, 1, 0) & np.roll(rock, -1, 0), nb[0], st)
        step = st
    # holes: background pixels enclosed by rock (seen through a slot) become the darkest step
    bgm = g['bg'] & ~rock
    lab, n = ndi.label(bgm)
    if n:
        edge_ids = np.unique(np.concatenate([lab[0], lab[:, 0], lab[:, -1], lab[-1]]))
        sizes = ndi.sum(np.ones_like(lab), lab, np.arange(1, n + 1))
        small = np.isin(lab, [i + 1 for i in range(n) if sizes[i] < 0.002 * lab.size and (i + 1) not in edge_ids])
        rock = rock | small
        step = np.where(small, 0, step)
        g['bg'] = g['bg'] & ~small
        g['unit'] = np.where(small, g['unit'].max(), g['unit'])
    step = np.where(rock, step, -1)
    return step, lit, lines, p""")
open('src/genge.py','w').write(s)
EOF
nice -n 10 python3 src/z8.py && echo ok

# output:
ok

cd ~/work/pixelart-studies/books/genge; sed -i '' "s/    bgm = g\['bg'\] & ~rock/    bgm = ~rock \& ~(g['ground'] \& g['hit'])/" src/genge.py; grep -n "bgm =" src/genge.py; nice -n 10 python3 src/z8.py && echo ok

# output:
232:    bgm = ~rock & ~(g['ground'] & g['hit'])
ok

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/genge.py').read()
s=s.replace("""    unit = g['unit']
    for k, r in enumerate(ramps):
        m = (unit == k) & (step >= 0)""","""    unit = g['unit']
    if ((unit < 0) & (step >= 0)).any():      # centre-sample id missed the rock: take a neighbour's unit
        uu = unit.copy()
        for _ in range(3):
            bad = (uu < 0) & (step >= 0)
            if not bad.any():
                break
            nb = np.maximum.reduce([np.roll(uu, 1, 0), np.roll(uu, -1, 0), np.roll(uu, 1, 1), np.roll(uu, -1, 1)])
            uu = np.where(bad, nb, uu)
        unit = np.where((uu < 0) & (step >= 0), 0, uu)
    for k, r in enumerate(ramps):
        m = (unit == k) & (step >= 0)""")
open('src/genge.py','w').write(s)
EOF
nice -n 10 python3 src/z8.py && echo ok

# output:
ok

cd ~/work/pixelart-studies/books/genge; for s in "e1_mastercopy_downsample.py" "e1_tor.py 5 3" "e2_light.py" "e3_beds.py 2" "e6_canyon.py" "t_subjects.py"; do nice -n 10 python3 src/$s 2>&1 | grep -v Warn | tail -2; done; cp stages/e1/e1a_downsample.png sheets/A0_mastercopy_downsample.png; cp stages/e1/e1b_tor.png sheets/A_tor_scales.png; cp stages/e2/e2_light.png sheets/B_shadow_vs_AO.png; cp stages/e3/e3_beds.png sheets/C_bedded_cliff.png; cp stages/e6/e6_canyon.png sheets/D_basalt_canyon.png; cp stages/tests/ruby_stack.png sheets/E_ruby_stack.png; cp stages/tests/sucia.png sheets/F_sucia.png; ls sheets; du -sh .

# output:
python3: can't open file '~/work/pixelart-studies/books/genge/src/e1_tor.py 5 3': [Errno 2] No such file or directory
python3: can't open file '~/work/pixelart-studies/books/genge/src/e3_beds.py 2': [Errno 2] No such file or directory
A0_mastercopy_downsample.png
A_tor_scales.png
B_shadow_vs_AO.png
C_bedded_cliff.png
D_basalt_canyon.png
E_ruby_stack.png
F_sucia.png
Z_tor_8x.png
 95M    .

cd ~/work/pixelart-studies/books/genge; nice -n 10 python3 src/e1_tor.py 5 3 && nice -n 10 python3 src/e3_beds.py 2 && cp stages/e1/e1b_tor.png sheets/A_tor_scales.png && cp stages/e3/e3_beds.png sheets/C_bedded_cliff.png && echo ok

# output:
ok

cd ~/work/pixelart-studies/books/genge; python3 - <<'EOF'
s=open('src/genge.py').read()
s=s.replace("""    H, W = g['rock'].shape
    rock = g['rock']
    v, lit = value(g, L, p)""","""    H, W = g['rock'].shape
    rock = g['rock']
    # scale LOD for line classes (A0: Genge's lines die by 0.35 m/px; value masses carry farther):
    # drawn bed and joint lines only at close range; above line_lod_mpp, bedding is carried only by
    # the lip/soffit value pair of protruding beds at least lod_min_px thick
    mpp = p.get('mpp', 0.1)
    if mpp > p.get('line_lod_mpp', 0.2):
        p['bed_lines'] = False
        p['joints'] = False
    v, lit = value(g, L, p)""")
s=s.replace("""        sof = bnd_below & comp_above & (comp < 0.5)""","""        sof = bnd_below & comp_above & (comp < 0.5)
        if 'bed_px' in g and mpp > p.get('line_lod_mpp', 0.2):
            sof &= np.roll(g['bed_px'], 1, 0) >= p.get('lod_min_px', 2.0)""")
s=s.replace("""        lipm = bnd_below & comp_here & lit & (np.roll(comp, 1, 0) < 0.5)""","""        lipm = bnd_below & comp_here & lit & (np.roll(comp, 1, 0) < 0.5)
        if 'bed_px' in g and mpp > p.get('line_lod_mpp', 0.2):
            lipm &= g['bed_px'] >= p.get('lod_min_px', 2.0)""")
open('src/genge.py','w').write(s)
# West Bay: joint slots nearly closed; furrows carry
for f in ['src/e3_beds.py','src/e2_light.py']:
    s=open(f).read()
    s=s.replace("joint_mean=6.0,\n            joint_keep=0.15,","joint_mean=6.0,\n            joint_keep=0.15, joint_gap=0.02,")
    s=s.replace("joint_mean=6.0, joint_keep=0.15,","joint_mean=6.0, joint_keep=0.15, joint_gap=0.02,")
    open(f,'w').write(s)
# canyon: rimrock only on upper walls, stronger haze
s=open('src/geology.py').read()
s=s.replace("rock = (ssm > 36) & comp & (wall > 0.12) & exposed_flow[k]","rock = (ssm > 36) & comp & (Zs > 0.62 * depth) & exposed_flow[k]")
open('src/geology.py','w').write(s)
s=open('src/e6_canyon.py').read()
s=s.replace("haze=0.3)","haze=0.55)")
open('src/e6_canyon.py','w').write(s)
# tor: low wide left-heavy pile
for f in ['src/e1_tor.py','src/z8.py','src/e2_light.py']:
    s=open(f).read()
    s=s.replace("sheet_mean=0.75, joint_mean=2.2, recess=0.5, stacks=3, height=7.0","sheet_mean=0.7, joint_mean=3.0, recess=0.6, stacks=2, height=5.5, taper=0.2, seed_extra=0")
    s=s.replace(", seed_extra=0","")
    open(f,'w').write(s)
# ruby: broad triangle with stepped right side
s=open('src/t_subjects.py').read()
s=s.replace("half = W * 0.36 * (1 - 0.55 * np.clip(V / V.max(), 0, 1) ** 1.4)","half = W * 0.47 * (1 - 0.9 * np.clip(V / V.max(), 0, 1) ** 0.9)")
s=s.replace("c = W * 0.5 + (V / V.max()) * W * 0.08","c = W * 0.5 + np.clip(V / V.max(), 0, 1) * W * 0.12")
open('src/t_subjects.py','w').write(s)
EOF
grep -n "def tor" src/geology.py; grep -n "taper" src/e1_tor.py | head -3

# output:
42:def tor(width=12.0, height=6.5, ds=0.05, seed=3, sheet_mean=0.95, sheet_sigma=0.4, bend=0.25,
7:geo = geology.tor(seed=int(sys.argv[1]) if len(sys.argv) > 1 else 3, sheet_mean=0.7, joint_mean=3.0, recess=0.6, stacks=2, height=5.5, taper=0.2)

cd ~/work/pixelart-studies/books/genge; for s in "e1_tor.py 5 3" "e2_light.py" "e3_beds.py 2" "e6_canyon.py" "t_subjects.py" "z8.py"; do set -- $s; nice -n 10 python3 src/$1 $2 $3 2>&1 | grep -v Warn | tail -2; done; cp stages/e1/e1b_tor.png sheets/A_tor_scales.png; cp stages/e2/e2_light.png sheets/B_shadow_vs_AO.png; cp stages/e3/e3_beds.png sheets/C_bedded_cliff.png; cp stages/e6/e6_canyon.png sheets/D_basalt_canyon.png; cp stages/tests/ruby_stack.png sheets/E_ruby_stack.png; cp stages/tests/sucia.png sheets/F_sucia.png; echo done

# output:
python3: can't open file '~/work/pixelart-studies/books/genge/src/e1_tor.py 5 3': [Errno 2] No such file or directory
python3: can't open file '~/work/pixelart-studies/books/genge/src/e3_beds.py 2': [Errno 2] No such file or directory
done

cd ~/work/pixelart-studies/books/genge; nice -n 10 python3 src/e1_tor.py 5 3 && nice -n 10 python3 src/e3_beds.py 2 && cp stages/e1/e1b_tor.png sheets/A_tor_scales.png && cp stages/e3/e3_beds.png sheets/C_bedded_cliff.png && echo ok

# output:
ok

cd ~/work/pixelart-studies/books/genge; sed -i '' 's/furrows=1.0, furrow_depth=0.7, furrow_spacing=2.0/furrows=1.0, furrow_depth=0.3, furrow_spacing=3.0/' src/e3_beds.py src/e2_light.py; nice -n 10 python3 src/e3_beds.py 2 && nice -n 10 python3 src/e2_light.py && cp stages/e3/e3_beds.png sheets/C_bedded_cliff.png && cp stages/e2/e2_light.png sheets/B_shadow_vs_AO.png; echo ok

# output:
ok

# What Genge teaches, in pixels: lessons, tests, honest comparisons

The book's model, in its own terms and with page numbers, is in `MODEL.md`. This file covers what I
did with it: the translations, the tests, the dead ends, and what didn't carry.

**Where it stands, plainly.** Three of Genge's principles carry into our medium and changed the
pictures:

1. structure survives shadow: ambient occlusion instead of cast shadow;
2. drawn lines depend on scale;
3. geology decides the geometry (bedding statistics, thin competent beds, fluting).

Some rules only work at close range: his line hierarchy, the impersistent lines, and the joint
budget. Two don't carry at all, and I say why below. None of this is at Ferrari's 8× bar. The basalt
canyon and the sea stack are the weakest tests. A cold critic ranked the sheets
F > A0 > Z > A > E > B > C > D before the last fixes. After those fixes, C (West Bay) went from worst
to one of the better sheets at landscape scale.


cd ~/work/pixelart-studies/books/genge
nice -n 10 python3 src/e1_mastercopy_downsample.py   # sheets A0: plate vs photo, box-downsampled
nice -n 10 python3 src/e1_tor.py 5 3                 # A: Rough Tor, 3 scales (seed, supersample)
nice -n 10 python3 src/z8.py                         # Z: tor at 8x
nice -n 10 python3 src/e2_light.py                   # B: cast shadow vs AO, 4 lights, 2 subjects
nice -n 10 python3 src/e3_beds.py 2                  # C: bedded cliff, schematic vs truthful
nice -n 10 python3 src/e6_canyon.py                  # D: basalt canyon, landscape + 30 deg iso
nice -n 10 python3 src/t_subjects.py                 # E, F: Ruby Beach stack, Sucia sandstone


Everything is numpy plus numba on one core, taking 5 to 60 s per sheet, with peak memory under
1.5 GB. The folder is about 100 MB, most of it plate renders and reference photos.

---

## 0. Reading

I read the whole text, all 16 chapters and the appendix, and looked at about 50 figure pages at
150 dpi, 9 of them again at 300 dpi (`plates/`). The drawings changed what I took from the text in
three places.

- **Fig. 2.9 (Rough Tor)** is the whole shading chapter in one plate. Its labels are the rules:
  inset areas are shadowed, the shadow marks the top of an overhang, deep fractures are dark, and
  the darkest shadow sits in the most inset areas. It is also a granite tor, so it is directly our
  Enchantments and Skykomish rock.
- **Fig. 1.1 (Ladram Bay)** makes the book's case. In the photo the sea stack is a dark red lump;
  in the sketch it is a stack of beds and ledges. The Ruby Beach stacks have the same problem in
  pixels.
- **Fig. 15.4 (the landscape block diagram)** is our 30° iso. It is drawn with ridge lines, the
  lower edges of cliffs, short steep strokes on cliff faces, and darker tone at breaks of slope.

## 1. Master copies by honest downsampling (sheet A0) — the finding that organised everything

I box-downsampled Genge's own Rough Tor drawings (Figs 2.5 and 2.9b) to our three scales, and did
the same to a photo of the tor.

![A0](sheets/A0_mastercopy_downsample.png)

- **Close, 0.1 m/px:** his line drawing survives almost intact. It is a perfectly good pixel
  drawing.
- **Landscape, 0.35 m/px:** the lines turn into grey mush. The photo still reads, because what
  carries it is value masses: the dark crescent under each pillow of rock.
- **Region, 1 m/px:** the drawing is gone. The photo is a lumpy mass with a dark lower left.

**So Genge's line rules are close-range rules. His shading rule, the ambient-occlusion masses,
carries to landscape scale, and only silhouette and big masses carry to region scale.** This is also
why supersampling matters. A 0.1 m joint seen at 0.35 m/px should be a grey hint, not a black line,
so the G-buffer is rendered 3×3 and area-averaged before quantizing (`Scene.gbuffer(ss=3)`).

## 2. Translations, one principle at a time

The method is a heightfield in any frame (a cliff face, a tor, or terrain), built from the
geological statistics in `geology.py`. It is ray-marched to a G-buffer (`engine.py`, `render.py`)
and painted as an index canvas of steps 0–7: 0–3 is the shadow family and 4–7 the lit family, the
studio's convention (`genge.py`, `landscape.py`). Each palette ramp (`palette.py`) is built from
measured L\*C\*h values, so swapping the palette relights the canvas. Every rule below is one
parameter, so it can be switched off and compared.

| # | Genge (page) | pixel translation | verdict |
|---|---|---|---|
| 1 | Don't draw cast shadows, only ambient occlusion (p. 37, Fig 2.9) | `shadows='mass'`: keep cast shadows only as big masses (≥ 6% of the rock, or 2% in landscapes) and drop small ones. `ao_inset` takes up to about one step off inset areas without leaving the family (his grey layer). Occlusion above `ao_deep` goes to step 0 or 1 whatever the light (his black layer). | **Carries, strongest result** (sheet B). Back-lit and overcast, the photo-light cliff is a flat brown slab; with the rules, every bed and ledge stays readable. Under raking light, pure Genge (no cast shadows) loses a real value mass, which the critic caught. The pixel version keeps it and shows structure inside it. **Needs rethinking:** a game has sun, so the shadow stays as design and only the small, fragmenting shadows go. |
| 2 | Line width is emphasis: thick at the horizon, the foreground and protruding beds, thin and clean elsewhere (Fig 2.5; p. 106) | There are no lines in pixels, only value edges. The "thick" line becomes a pair: a lit lip on the top edge of a protruding bed (`lip`) and a darkest 1-px soffit under it (`soffit`). The foreground contact becomes a step-0 pocket (`contact`), used at close range only. At 1 m/px it made a black band and was cut. The skyline gets a clean, undithered edge, darkened one step only where rock and sky are close in value. | **Carries at close and landscape range** (Sucia, sheet F, was the critic's best cell). The lip and soffit pair survives downsampling where a dark line alone does not. |
| 3 | Impersistent lines, with stroke length matched to the object's scale (p. 28) | `bed_lines`: bed boundaries broken into dashes. Mean run is 1.8× the bed's thickness in pixels, with gaps. | **Close range only.** At 0.35 m/px and beyond, dashes turn into Morse code and barcode (sheet C, first round). So `line_lod_mpp=0.2` turns drawn lines off above 0.2 m/px, and bedding is then carried only by the lip and soffit pair of beds at least 2 px thick. That is Genge's own "only a proportion of the layers can be drawn" (p. 98). |
| 4 | Don't over-emphasise joints (Fig 4.6; p. 61) | `joints`: longest joint segments first, capped at `joint_ratio` × bedding-line pixels. In the geometry, most joints stop at bed boundaries (T-junctions) and a few master joints run through several beds. | Carries, but the geometry mattered more than the cap. Open joint slots of 0.1 m made a tartan grid at 0.35 m/px until supersampling and near-closed slots (`joint_gap`). Joints are also off above 0.2 m/px. |
| 5 | Beds vary in thickness vertically and laterally, pinch out, thin upward, competent beds protrude (pp. 98, 141–144; Figs 4.7, 5.12) | `bedded_face`: lognormal thicknesses, `thin_up`, laterally wandering boundaries (`taper`), boundaries that die out and merge the beds (`pinch`). Thin competent "doggers" (`comp_mean`) sit in weak sand, and weathering furrows run below the ledges. | **Carries, and it was the big lever** on West Bay. The first two rounds were brick walls because I had drawn every bed as competent and every joint as open. That is exactly the schematic-bedding mistake Genge names in Fig 4.7. |
| 6 | Furrows show the cliff surface (p. 98) | Shallow vertical grooves that fade out at random heights. | Carries once they are shallow (0.3 m). When they were deep (0.7 m), AO turned them into dark bars. His furrows are light marks, not slots. |
| 7 | Landscapes: horizons first, bedding traces follow topography, Vs up the gullies (Ch 7; p. 229) | `canyon`: horizontal flows. Each flow's competent base forms a riser and its weak top a bench, so every trace follows the contours and Vs by construction. `landscape.py` adds traces (prominent flows only), gully lines, a lip and foot on rimrock, haze by depth, and local-horizon separation. | **Partly.** The rule of Vs emerges for free, which is the right way round: geology first, drawing second. The picture is still weak. Gully and spur shadows run as a regular zebra, local horizons barely separate, and my first versions painted rimrock on every flow (stairs), which the critic rightly called a barcode. Now rimrock is exposed only near crests (`Zs > 0.62 depth`) and only on some flows. |
| 8 | Many small things: draw enough to show the texture (~30 of hundreds), slightly enlarged (p. 201) | Not built as a separate painter. The despeckle rule (no lone pixels) and the ≥ 2 px gate on lips and soffits are the pixel form of "draw enough, slightly enlarged". | Thin. Crystals and clasts at 0.1 m/px were not tested. |
| 9 | Colour: a flat base per unit, exaggerated between units, no blending (p. 35) | `palette.unit_ramps(mat, n, spread)` makes sibling ramps pushed apart in L, C and hue. Texture is a second layer streaked along the bedding (`streak`, `bed_angle`), like his Fig 2.8c and Fig 16.7c. | The ramps exist but no sheet shows a multi-unit rock. That test (Yakima flow tops red against colonnade grey) is still to do. The along-bedding streak works: on dipping rock it is what makes the dip read (sheet E). |

## 3. Tests on the real subjects (before | after beside the photo)

**Rough Tor, a granite tor (sheets A and Z), stands in for Enchantments and Skykomish granite at
close range.** The pillows are lens shaped, with a harder undercut below
(`undercut`: "shadow indicates top of overhang"), on curved sheeting joints.

![A](sheets/A_tor_scales.png)
![Z](sheets/Z_tor_8x.png)

- **Right:** at 0.1 m/px it reads as a stacked-pillow tor, with dark slots and lit convex tops.
- **Wrong:** the silhouette is a symmetric block where the photo is a low, left-heavy, overhanging
  pile. The soffits are even straight bands where the photo has tapering crescents. There is still
  a brick rhythm.
- The rules make little visible difference on the tor in front light, because the cast shadows here
  are the occlusion. That is an honest null.

**West Bay, the Bridport Sands cliff (sheet C), is the method test for the Olympic coast cliffs.**

![C](sheets/C_bedded_cliff.png)

- At landscape scale the truthful geometry plus the rules gives horizontal banding, soft fluting and
  thinning upward. That is close to the whole-cliff photo, and far from the schematic and
  joints-emphasised columns.
- At region scale all versions are near-flat. The critic judged the plain schematic closest. That
  is fine: at 1 m/px nothing but masses should survive.
- **Missing:** the grass cap. The gold face against the green slope is what makes the photo read.

**Ruby Beach sea stack (sheet E), on the Olympic coast.**

![E](sheets/E_ruby_stack.png)

- The after column shows the steep dip, which the photo hides in darkness, as Ladram Bay does. That
  is Genge's Fig 1.1 argument, working.
- On dark rock, though, drawn dark lines pile up into hatching, which is his own mistake from
  Fig 6.9. The dark-rock variant, with lips instead of lines, is better.
- **Rule:** on rock whose lit face is below L\* ~45, carry the bedding with light lips, not dark
  lines.
- The silhouette is still a thumb, not the photo's broad stepped triangle.

**Sucia sandstone in the San Juans (sheet F).** At close range the lip and soffit on the upper ledge
and the overhang shadow read as real 3D. At landscape scale the cell is a wash. The photo's bench
dips steadily at about 15° with one deep overhang; mine undulates.

![F](sheets/F_sucia.png)

**Basalt canyon, Yakima (sheet D).** This is the rule of Vs, contour-following traces, spur and
gully light, and haze. It reads as a dry canyon wall under raking light, but it is the weakest
sheet. See row 7 of the table.

![D](sheets/D_basalt_canyon.png)

**Not done:**

- the Dragontail real-terrain match;
- a voxsim G-buffer from a place lane;
- other students' scenes;
- a master copy of the Fig 15.4 block diagram.

The painter takes the same per-pixel inputs voxsim already has (normal, sky occlusion, sun shadow,
bed and block ids, depth), so a lane can try it on a G-buffer dump.

## 4. What didn't carry, and why

- **"Don't draw cast shadows" taken literally.** It is right for a notebook and wrong for a game.
  Shadow is our value design and our time of day. What carries is the reason behind it: structure
  has to stay readable inside the shadow.
- **Line hierarchy beyond close range.** Lines are sub-pixel at 0.35 m/px. Genge's own plates
  downsampled prove it (A0).
- **Dashed lines for interpretation (p. 28), and labels, scale bars and looking-towards bearings.**
  These are notebook conventions with no pixel meaning. The nearest in-game equivalent is a debug
  overlay.
- **Vertical exaggeration** is a mistake in sketches but normal in schematic block diagrams (p. 257).
  I kept true scale throughout.

## 5. Dead ends, in order

1. The tor as flat bricks: the pillows had no lens profile and no undercut.
2. Point-sampled rays: every 0.1 m joint became a black line, a tartan grid. Fixed by supersampling.
3. West Bay with every bed competent and every joint open: a brick wall, twice. Then deep furrows
   that read as dark bars.
4. A contact line at 1 m/px: a black band.
5. Canyon terraces taking 85% of each flow's relief: giant stairs. Then rimrock on every flow:
   barcode.
6. Dark bed lines on dark coast rock: hatching.
7. Isolated dither pixels and centre-sample id misses: Morse dots and blue holes. Fixed by
   `despeckle` and the hole fill.

## 6. Rules for other students and the voxsim lanes

1. **Geology before drawing.** Build the rock from its statistics:
   - beds lognormal, thinning upward where the rock does, wandering along strike, some pinching out;
   - thin competent beds in weak ones;
   - joints few, perpendicular to bedding, mostly stopping at bed boundaries.

   Most of what reads as "true" comes from here, not from paint.
2. **Honest pixels.** Supersample, then area-average before quantizing. A feature narrower than a
   pixel is a one-step hint, never a black line.
3. **Structure survives shadow.** Keep cast shadows only as big masses (≥ 2–6% of the object).
   Inside every family, inset areas lose about one step, and the most inset go to the darkest step.
   Check every scene back-lit and overcast: if it goes flat, this rule is missing.
4. **Lines by scale.**
   - At ≤ 0.2 m/px: impersistent bed lines (run ≈ 2× bed thickness), a lip and soffit on protruding
     beds, a few joints (≤ 0.35× bedding pixels), and a contact pocket.
   - Above that: only the lip and soffit on beds at least 2 px thick.
   - At region scale: only silhouette and masses.
5. **On dark rock (lit L\* below ~45), use light lips, not dark lines.**
6. **Landscapes:**
   - horizontal beds make contour-following traces that V up the gullies, so let the geometry draw
     them;
   - draw traces only for the prominent flows;
   - expose rimrock only near crests;
   - separate each local horizon with a lighter, cooler value.
7. **Colour:** a flat base per unit, with units pushed apart (`unit_ramps`), and texture streaked
   along the bedding (`bed_angle`), never stippled.

## 7. Code

In `src/`:

- `engine.py`: frames, the heightfield ray-march, horizon-based sky occlusion and sun shadow.
- `render.py`: `Scene` and `gbuffer(ss=)`, plus the face frame and camera.
- `geology.py`: `tor`, `bedded_face`, `canyon`, all parameters documented in the docstrings.
- `genge.py`: `paint(g, L, params)` → `(step, lit, lines, p)`, and `colorize`. The parameter table
  is in the module docstring.
- `landscape.py`: the terrain painter.
- `palette.py`: Lab-anchored ramps and `unit_ramps`.

The minimal port for a lane is `genge.paint` on its own G-buffer. It needs these keys:
`rock, ground, bg, hit, N, ao, shadow, bed, block, unit, comp, bed_px, mpp`.

## 8. Next, in order

1. Masses and silhouette first, checked at 1 m/px against the photo before any rule (the critic's
   first fix): a left-heavy tor, a broad stepped Ruby stack, the West Bay grass cap.
2. Crescent soffits of varying thickness, and irregular partial joints.
3. The canyon: branched, tapering gullies, and value steps between successive ridges.
4. A multi-unit colour test on Yakima flows.
5. A voxsim G-buffer from a place lane.
6. The Dragontail terrain with these rules.

cd ~/work/pixelart-studies/books/genge; rm -rf src/__pycache__; du -sh . plates refs stages sheets; ls sheets

# output:
93M .
 63M    plates
 18M    refs
1.4M    stages
1.2M    sheets
A0_mastercopy_downsample.png
A_tor_scales.png
B_shadow_vs_AO.png
C_bedded_cliff.png
D_basalt_canyon.png
E_ruby_stack.png
F_sucia.png
Z_tor_8x.png