GameCult/Ratatoskr
Watch
9
0
Fork
You've already forked Ratatoskr
 0
Code Issues Pull requests Projects Releases 3 Packages Wiki Activity Actions

ridge-bake: move the A+D baker, tests and mutation runner into scripts/

Browse source 
Moved from the mockup workspace: ridge.py, sun.py, places.json, tests. core.py is
dropped (ridge.py only imported its save_png, unused; the rest carried hard-coded
paths). places.json is read beside ridge.py; tiles and caches come from RIDGE_TILES
and RIDGE_WORK.

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
This commit is contained in:
Nibu 2026-11-13 23:19:01 +00:00
commit bffe988f05
7 changed files with 1647 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
scripts/ridge-bake/__pycache__/ridge.cpython-313.pyc Normal file
View file

Binary file not shown.

657
scripts/ridge-bake/places.json Normal file
View file

@ -0,0 +1,657 @@
{
"_doc": "Place presets for the A+D baker (ridge.py). eye = observer lat/lon; eye_h_m = eye above the local DEM surface (ground_override_m replaces the DEM ground elevation, used for Vienna to keep round 1's 183.3 m); look_at = target lat/lon whose bearing centres the 72 deg window (or az_deg); rmax_km = ring radius (decides the tile list: every 1x1 deg GLO-30 tile whose rectangle intersects the circle); floor_agl_m = authored eye height above the notional valley floor used only for the floor-plane colour gradient; s_ref_deg = front slope that counts as steep for the painter; tz_h = standard-time offset of the local clock on the sun date.",
"date": "2026-09-23",
"presets": [
{
"id": "vienna",
"name": "Vienna, looking SW to Rax/Schneeberg",
"kind": "far range over a flat foreground",
"lat": 48.1929826,
"lon": 16.3534237,
"eye_h_m": 30.0,
"ground_override_m": 183.3,
"az_deg": 214.0,
"rmax_km": 95,
"floor_agl_m": 30,
"s_ref_deg": 8,
"tz_h": 2,
"tz_name": "CEST",
"viewpoint_source": {
"status": "existing observer of the round-1 mockup (the scene camera used by the game), not photo-sourced; eye 30 m above ground as in round 1"
}
},
{
"id": "matterhorn",
"name": "Matterhorn from the Gornergrat (photo cell, 9.7 km)",
"kind": "lone horn, close, 27 deg of relief",
"lat": 45.983139,
"lon": 7.782466,
"eye_h_m": 1.6,
"look_at": [
45.97666,
7.65678
],
"rmax_km": 60,
"floor_agl_m": 992.3298828124998,
"s_ref_deg": 30,
"tz_h": 2,
"tz_name": "CEST",
"initial_guess_not_used": {
"lat": 45.9856,
"lon": 7.7452
},
"viewpoint_source": {
"status": "sourced",
"source": "Wikimedia Commons geosearch (namespace 6), geotagged photos",
"query": "geosearch ns=6 gslimit=500 gsradius=10000 around 25 centres (peak + 3 rings x 8 bearings, radii 3600-12000 m), then EXIF GPSImgDirection of every file in range via prop=imageinfo&iiprop=commonmetadata",
"files_found": 5003,
"in_range": 3241,
"with_exif_heading": 205,
"facing_the_peak": 52,
"pool": "photos whose camera heading is within 30 deg of the bearing to the summit",
"rule": "densest 1 km cell of photos 3000-12000 m from the summit whose EXIF camera heading faces the summit (within 30 deg) with the summit visible on the DEM (straight line, first 250 m ignored, 8 m clutter allowed); viewpoint = median coordinate of the cell",
"median_camera_heading_deg": 245.8,
"photos_in_cell": 7,
"cells_tried": [
{
"n_facing": 7,
"lat": 45.98314,
"lon": 7.78247,
"peak_visible": true
}
],
"distance_to_peak_m": 9749,
"bearing_to_peak_deg": 265.8,
"dem_ground_m": 3090.7,
"eye_height_rule": "DEM ground at the viewpoint + 1.6 m (the source gives no eye height)",
"sample_files": [
{
"title": "File:Gornergrat-04.jpg",
"lat": 45.983061,
"lon": 7.782031,
"heading": 256.0,
"url": "https://commons.wikimedia.org/wiki/File%3AGornergrat-04.jpg"
},
{
"title": "File:Gornergrat-03.jpg",
"lat": 45.983058,
"lon": 7.782038,
"heading": 245.8,
"url": "https://commons.wikimedia.org/wiki/File%3AGornergrat-03.jpg"
},
{
"title": "File:Gornergrat-07.jpg",
"lat": 45.982653,
"lon": 7.782466,
"heading": 250.2,
"url": "https://commons.wikimedia.org/wiki/File%3AGornergrat-07.jpg"
},
{
"title": "File:Gornergratsee 2021 10.jpg",
"lat": 45.986958,
"lon": 7.780303,
"heading": 243.4,
"url": "https://commons.wikimedia.org/wiki/File%3AGornergratsee_2021_10.jpg"
},
{
"title": "File:Early morning Eiger.jpg",
"lat": 45.983242,
"lon": 7.783361,
"heading": 254.7,
"url": "https://commons.wikimedia.org/wiki/File%3AEarly_morning_Eiger.jpg"
},
{
"title": "File:Klein Matterhorn mit Wolken an Theodulhorn und Furgghorn \u00fcber dem Oberen Theodulgletscher- a view from Gornergrat 2021 02 (cropped).jpg",
"lat": 45.983139,
"lon": 7.784172,
"heading": 242.1,
"url": "https://commons.wikimedia.org/wiki/File%3AKlein_Matterhorn_mit_Wolken_an_Theodulhorn_und_Furgghorn_%C3%BCber_dem_Oberen_Theodulgletscher-_a_view_from_Gornergrat_2021_02_%28cropped%29.jpg"
},
{
"title": "File:Breithorn and neighbours from Gornergrat 2021 02.jpg",
"lat": 45.983139,
"lon": 7.784172,
"heading": 242.1,
"url": "https://commons.wikimedia.org/wiki/File%3ABreithorn_and_neighbours_from_Gornergrat_2021_02.jpg"
}
],
"exif_headings": [
{
"title": "File:Gornergrat-04.jpg",
"heading": 256.0
},
{
"title": "File:Gornergrat-03.jpg",
"heading": 245.8
},
{
"title": "File:Gornergrat-07.jpg",
"heading": 250.2
},
{
"title": "File:Gornergratsee 2021 10.jpg",
"heading": 243.4
},
{
"title": "File:Early morning Eiger.jpg",
"heading": 254.7
},
{
"title": "File:Klein Matterhorn mit Wolken an Theodulhorn und Furgghorn \u00fcber dem Oberen Theodulgletscher- a view from Gornergrat 2021 02 (cropped).jpg",
"heading": 242.1
},
{
"title": "File:Breithorn and neighbours from Gornergrat 2021 02.jpg",
"heading": 242.1
}
],
"n_with_exif_heading": 7,
"summit_dem_m": 4324.9,
"summit_nominal": [
45.97639,
7.65833
]
}
},
{
"id": "fuji",
"name": "Fuji from the Shiraito Falls side, 13.8 km WSW (photo cell)",
"kind": "lone cone, smooth",
"lat": 35.312467222222494,
"lon": 138.587223194445,
"eye_h_m": 1.6,
"look_at": [
35.36061,
138.72731
],
"rmax_km": 60,
"floor_agl_m": 138.6306213378906,
"s_ref_deg": 20,
"tz_h": 9,
"tz_name": "JST",
"initial_guess_not_used": {
"lat": 35.4994,
"lon": 138.7686
},
"viewpoint_source": {
"status": "sourced",
"source": "Wikimedia Commons geosearch (namespace 6), geotagged photos",
"query": "geosearch ns=6 gslimit=500 gsradius=10000 around 25 centres (peak + 3 rings x 8 bearings, radii 9000-30000 m), then EXIF GPSImgDirection of every file in range via prop=imageinfo&iiprop=commonmetadata",
"files_found": 6494,
"in_range": 4280,
"with_exif_heading": 588,
"facing_the_peak": 132,
"pool": "photos whose camera heading is within 30 deg of the bearing to the summit",
"rule": "densest 1 km cell of photos 8000-30000 m from the summit whose EXIF camera heading faces the summit (within 30 deg) with the summit visible on the DEM (straight line, first 250 m ignored, 8 m clutter allowed); viewpoint = median coordinate of the cell",
"median_camera_heading_deg": 64.4,
"photos_in_cell": 8,
"cells_tried": [
{
"n_facing": 15,
"lat": 35.36248,
"lon": 138.86287,
"peak_visible": false
},
{
"n_facing": 8,
"lat": 35.31247,
"lon": 138.58722,
"peak_visible": true
}
],
"distance_to_peak_m": 13808,
"bearing_to_peak_deg": 67.2,
"dem_ground_m": 486.9,
"eye_height_rule": "DEM ground at the viewpoint + 1.6 m (the source gives no eye height)",
"sample_files": [
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 03.JPG",
"lat": 35.31114,
"lon": 138.5886,
"heading": 48.0,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_03.JPG"
},
{
"title": "File:Mount Fuji, view from Shiraito no Taki.jpg",
"lat": 35.313589,
"lon": 138.588394,
"heading": 68.8,
"url": "https://commons.wikimedia.org/wiki/File%3AMount_Fuji%2C_view_from_Shiraito_no_Taki.jpg"
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 02.JPG",
"lat": 35.313248055556,
"lon": 138.58760805556,
"heading": 64.2,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_02.JPG"
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 07.JPG",
"lat": 35.312428055556,
"lon": 138.58723138889,
"heading": 88.8,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_07.JPG"
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 08.JPG",
"lat": 35.312418055556,
"lon": 138.587215,
"heading": 64.0,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_08.JPG"
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 04.JPG",
"lat": 35.312478055556,
"lon": 138.58717805556,
"heading": 62.2,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_04.JPG"
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 06.JPG",
"lat": 35.312456388889,
"lon": 138.58715805556,
"heading": 65.2,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_06.JPG"
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 05.JPG",
"lat": 35.312479722222,
"lon": 138.58714138889,
"heading": 64.5,
"url": "https://commons.wikimedia.org/wiki/File%3A%E5%AF%8C%E5%A3%AB%E5%B1%B1%EF%BC%88%E3%81%B5%E3%81%98%E3%81%95%E3%82%93%E3%80%81%E8%8B%B1%E8%AA%9E%EF%BC%9AMount_Fuji%EF%BC%89_05.JPG"
}
],
"exif_headings": [
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 03.JPG",
"heading": 48.0
},
{
"title": "File:Mount Fuji, view from Shiraito no Taki.jpg",
"heading": 68.8
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 02.JPG",
"heading": 64.2
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 07.JPG",
"heading": 88.8
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 08.JPG",
"heading": 64.0
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 04.JPG",
"heading": 62.2
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 06.JPG",
"heading": 65.2
},
{
"title": "File:\u5bcc\u58eb\u5c71\uff08\u3075\u3058\u3055\u3093\u3001\u82f1\u8a9e\uff1aMount Fuji\uff09 05.JPG",
"heading": 64.5
}
],
"n_with_exif_heading": 8,
"summit_dem_m": 3765.5,
"summit_nominal": [
35.36061,
138.72731
]
}
},
{
"id": "halfdome",
"name": "Half Dome from the Yosemite Valley floor (Sentinel Bridge photo cell)",
"kind": "granite wall across a valley",
"lat": 37.745833333333,
"lon": -119.592606,
"eye_h_m": 1.6,
"look_at": [
37.74435,
-119.53426
],
"rmax_km": 40,
"floor_agl_m": 5.5442382812500455,
"s_ref_deg": 30,
"tz_h": -7,
"tz_name": "PDT",
"initial_guess_not_used": {
"lat": 37.7305,
"lon": -119.5734
},
"viewpoint_source": {
"status": "sourced",
"source": "Wikimedia Commons geosearch (namespace 6), geotagged photos",
"query": "geosearch ns=6 gslimit=500 gsradius=10000 around 25 centres (peak + 3 rings x 8 bearings, radii 2700-9000 m), then EXIF GPSImgDirection of every file in range via prop=imageinfo&iiprop=commonmetadata",
"files_found": 2505,
"in_range": 1788,
"with_exif_heading": 181,
"facing_the_peak": 64,
"pool": "photos whose camera heading is within 30 deg of the bearing to the summit",
"rule": "densest 1 km cell of photos 2000-9000 m from the summit whose EXIF camera heading faces the summit (within 30 deg) with the summit visible on the DEM (straight line, first 250 m ignored, 8 m clutter allowed); viewpoint = median coordinate of the cell",
"median_camera_heading_deg": 91.6,
"photos_in_cell": 9,
"cells_tried": [
{
"n_facing": 9,
"lat": 37.73853,
"lon": -119.57124,
"peak_visible": false
},
{
"n_facing": 9,
"lat": 37.74583,
"lon": -119.59261,
"peak_visible": true
}
],
"distance_to_peak_m": 5138,
"bearing_to_peak_deg": 91.8,
"dem_ground_m": 1209.2,
"eye_height_rule": "DEM ground at the viewpoint + 1.6 m (the source gives no eye height)",
"sample_files": [
{
"title": "File:Yosemite National Park - 51700634719.jpg",
"lat": 37.746606,
"lon": -119.590317,
"heading": 87.0,
"url": "https://commons.wikimedia.org/wiki/File%3AYosemite_National_Park_-_51700634719.jpg"
},
{
"title": "File:Half Dome, Yosemite National Park - 53324601178.jpg",
"lat": 37.748264,
"lon": -119.592606,
"heading": 98.8,
"url": "https://commons.wikimedia.org/wiki/File%3AHalf_Dome%2C_Yosemite_National_Park_-_53324601178.jpg"
},
{
"title": "File:Half Dome, Yosemite National Park - 53324824955.jpg",
"lat": 37.748322,
"lon": -119.592736,
"heading": 90.1,
"url": "https://commons.wikimedia.org/wiki/File%3AHalf_Dome%2C_Yosemite_National_Park_-_53324824955.jpg"
},
{
"title": "File:2013-09-19 18 35 25 Half Dome and the Merced River in Yosemite National Park.JPG",
"lat": 37.7435,
"lon": -119.58966666667,
"heading": 91.6,
"url": "https://commons.wikimedia.org/wiki/File%3A2013-09-19_18_35_25_Half_Dome_and_the_Merced_River_in_Yosemite_National_Park.JPG"
},
{
"title": "File:Half Dome from the Yosemite Valley floor in mid-Spring.jpeg",
"lat": 37.745833333333,
"lon": -119.59233333333,
"heading": 92.0,
"url": "https://commons.wikimedia.org/wiki/File%3AHalf_Dome_from_the_Yosemite_Valley_floor_in_mid-Spring.jpeg"
},
{
"title": "File:Sentinel Bridge & Half Dome.jpg",
"lat": 37.743419,
"lon": -119.590562,
"heading": 86.8,
"url": "https://commons.wikimedia.org/wiki/File%3ASentinel_Bridge_%26_Half_Dome.jpg"
},
{
"title": "File:Half Dome pre Super moon 2.jpg",
"lat": 37.745769,
"lon": -119.592972,
"heading": 106.5,
"url": "https://commons.wikimedia.org/wiki/File%3AHalf_Dome_pre_Super_moon_2.jpg"
},
{
"title": "File:Half Dome pre Super moon.jpg",
"lat": 37.745747,
"lon": -119.593033,
"heading": 97.1,
"url": "https://commons.wikimedia.org/wiki/File%3AHalf_Dome_pre_Super_moon.jpg"
},
{
"title": "File:Yosemite National Park 003.jpg",
"lat": 37.750814,
"lon": -119.598664,
"heading": 80.1,
"url": "https://commons.wikimedia.org/wiki/File%3AYosemite_National_Park_003.jpg"
}
],
"exif_headings": [
{
"title": "File:Yosemite National Park - 51700634719.jpg",
"heading": 87.0
},
{
"title": "File:Half Dome, Yosemite National Park - 53324601178.jpg",
"heading": 98.8
},
{
"title": "File:Half Dome, Yosemite National Park - 53324824955.jpg",
"heading": 90.1
},
{
"title": "File:2013-09-19 18 35 25 Half Dome and the Merced River in Yosemite National Park.JPG",
"heading": 91.6
},
{
"title": "File:Half Dome from the Yosemite Valley floor in mid-Spring.jpeg",
"heading": 92.0
},
{
"title": "File:Sentinel Bridge & Half Dome.jpg",
"heading": 86.8
},
{
"title": "File:Half Dome pre Super moon 2.jpg",
"heading": 106.5
},
{
"title": "File:Half Dome pre Super moon.jpg",
"heading": 97.1
},
{
"title": "File:Yosemite National Park 003.jpg",
"heading": 80.1
}
],
"n_with_exif_heading": 9,
"summit_dem_m": 2669.1,
"summit_nominal": [
37.74597,
-119.5329
]
}
},
{
"id": "amadablam",
"name": "Ama Dablam and the Everest massif, 4 km WSW of Tengboche near Khumjung (photo cell)",
"kind": "far high range behind a near spire",
"lat": 27.822000000000003,
"lon": 86.722666666667,
"eye_h_m": 1.6,
"look_at": [
27.8617,
86.86079
],
"rmax_km": 60,
"floor_agl_m": 829.8232421875,
"s_ref_deg": 30,
"tz_h": 5.75,
"tz_name": "NPT",
"initial_guess_not_used": {
"lat": 27.8369,
"lon": 86.7643
},
"viewpoint_source": {
"status": "sourced",
"source": "Wikimedia Commons geosearch (namespace 6), geotagged photos",
"query": "geosearch ns=6 gslimit=500 gsradius=10000 around 25 centres (peak + 3 rings x 8 bearings, radii 4800-16000 m), then EXIF GPSImgDirection of every file in range via prop=imageinfo&iiprop=commonmetadata",
"files_found": 1398,
"in_range": 1143,
"with_exif_heading": 54,
"facing_the_peak": 13,
"pool": "photos whose camera heading is within 30 deg of the bearing to the summit",
"rule": "densest 1 km cell of photos 4000-16000 m from the summit whose EXIF camera heading faces the summit (within 30 deg) with the summit visible on the DEM (straight line, first 250 m ignored, 8 m clutter allowed); viewpoint = median coordinate of the cell",
"median_camera_heading_deg": 79.6,
"photos_in_cell": 2,
"cells_tried": [
{
"n_facing": 2,
"lat": 27.822,
"lon": 86.72267,
"peak_visible": true
}
],
"distance_to_peak_m": 14299,
"bearing_to_peak_deg": 72.0,
"dem_ground_m": 3782.0,
"eye_height_rule": "DEM ground at the viewpoint + 1.6 m (the source gives no eye height)",
"sample_files": [
{
"title": "File:Thamserku (6).JPG",
"lat": 27.822666666667,
"lon": 86.722666666667,
"heading": 73.6,
"url": "https://commons.wikimedia.org/wiki/File%3AThamserku_%286%29.JPG"
},
{
"title": "File:Thamserku (7).JPG",
"lat": 27.821333333333,
"lon": 86.722666666667,
"heading": 85.7,
"url": "https://commons.wikimedia.org/wiki/File%3AThamserku_%287%29.JPG"
}
],
"exif_headings": [
{
"title": "File:Thamserku (6).JPG",
"heading": 73.6
},
{
"title": "File:Thamserku (7).JPG",
"heading": 85.7
}
],
"n_with_exif_heading": 2,
"summit_dem_m": 6659.6,
"summit_nominal": [
27.8617,
86.8614
]
}
},
{
"id": "trecime",
"name": "Tre Cime di Lavaredo from the Dreizinnenhuette area, 2.4 km N (photo cell)",
"kind": "cirque wall, spires",
"lat": 46.636906,
"lon": 12.310064,
"eye_h_m": 1.6,
"look_at": [
46.61806,
12.29583
],
"rmax_km": 40,
"floor_agl_m": 735.6748046875,
"s_ref_deg": 30,
"tz_h": 2,
"tz_name": "CEST",
"initial_guess_not_used": {
"lat": 46.5833,
"lon": 12.2536
},
"viewpoint_source": {
"status": "sourced",
"source": "Wikimedia Commons geosearch (namespace 6), geotagged photos",
"query": "geosearch ns=6 gslimit=500 gsradius=10000 around 25 centres (peak + 3 rings x 8 bearings, radii 3000-10000 m), then EXIF GPSImgDirection of every file in range via prop=imageinfo&iiprop=commonmetadata",
"files_found": 2028,
"in_range": 763,
"with_exif_heading": 207,
"facing_the_peak": 33,
"pool": "photos whose camera heading is within 30 deg of the bearing to the summit",
"rule": "densest 1 km cell of photos 1500-10000 m from the summit whose EXIF camera heading faces the summit (within 30 deg) with the summit visible on the DEM (straight line, first 250 m ignored, 8 m clutter allowed); viewpoint = median coordinate of the cell",
"median_camera_heading_deg": 194.6,
"photos_in_cell": 5,
"cells_tried": [
{
"n_facing": 5,
"lat": 46.63691,
"lon": 12.31006,
"peak_visible": true
}
],
"distance_to_peak_m": 2363,
"bearing_to_peak_deg": 207.4,
"dem_ground_m": 2403.0,
"eye_height_rule": "DEM ground at the viewpoint + 1.6 m (the source gives no eye height)",
"sample_files": [
{
"title": "File:Cime e Nuvole.jpg",
"lat": 46.63605,
"lon": 12.310064,
"heading": 200.3,
"url": "https://commons.wikimedia.org/wiki/File%3ACime_e_Nuvole.jpg"
},
{
"title": "File:DreiZinnenP.jpg",
"lat": 46.636770555556,
"lon": 12.310305555556,
"heading": 180.1,
"url": "https://commons.wikimedia.org/wiki/File%3ADreiZinnenP.jpg"
},
{
"title": "File:Soffio di sguardo.jpg",
"lat": 46.636906,
"lon": 12.311142,
"heading": 194.6,
"url": "https://commons.wikimedia.org/wiki/File%3ASoffio_di_sguardo.jpg"
},
{
"title": "File:2405-DreizinnenH.jpg",
"lat": 46.637782777778,
"lon": 12.309966388889,
"heading": 178.4,
"url": "https://commons.wikimedia.org/wiki/File%3A2405-DreizinnenH.jpg"
},
{
"title": "File:Dolomiten-Fingerkraut in den Sextener Dolomiten - N\u00e4he Drei-Zinnen-H\u00fctte (9801216636).jpg",
"lat": 46.642352,
"lon": 12.30829,
"heading": 200.0,
"url": "https://commons.wikimedia.org/wiki/File%3ADolomiten-Fingerkraut_in_den_Sextener_Dolomiten_-_N%C3%A4he_Drei-Zinnen-H%C3%BCtte_%289801216636%29.jpg"
}
],
"exif_headings": [
{
"title": "File:Cime e Nuvole.jpg",
"heading": 200.3
},
{
"title": "File:DreiZinnenP.jpg",
"heading": 180.1
},
{
"title": "File:Soffio di sguardo.jpg",
"heading": 194.6
},
{
"title": "File:2405-DreizinnenH.jpg",
"heading": 178.4
},
{
"title": "File:Dolomiten-Fingerkraut in den Sextener Dolomiten - N\u00e4he Drei-Zinnen-H\u00fctte (9801216636).jpg",
"heading": 200.0
}
],
"n_with_exif_heading": 5,
"summit_dem_m": 2746.9,
"summit_nominal": [
46.6186,
12.3025
]
}
}
]
}

507
scripts/ridge-bake/ridge.py Normal file
View file

@ -0,0 +1,507 @@
"""Round-2 baker: A's painter baked from D's sweep, driven by a place preset from places.json (ray sampling, the round-1 recommendation).
Offline baker for the ridge backdrop (scripts/ridge-bake); export.py drives it.
load_preset(id) -> configure() : observer, tiles, mosaic origin, ring radius, window azimuth, painter constants from places.json
stage0(detail) : ray march of all 6400 azimuth columns (curvature+refraction), optional authored sub-DEM foreground detail, cached
sweep() : D's per-column scan -> spans (with own height/run = per-peak slope) -> wrapped chains -> 7-knot layer cap
set_view(slots) : authored Y mapping: angle -> row, ground plane at the bottom (see set_view)
bake(slots) : sun-independent per-layer arrays
relight(layers, az, el) : shading + atmosphere + floor + sky
Vertical exaggeration is AUTHORED and derived per view by set_view (printed on the page)."""
import sys, time, math, json
import os, hashlib
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent))
import numpy as np, tifffile
def _env_dir(name, what):
v = os.environ.get(name)
if not v: raise RuntimeError(f"{name} is not set: it must name {what} (no default location is assumed)")
return Path(v)
B = _env_dir('RIDGE_WORK', 'the work directory (the mosaic / ray caches)')
COP = _env_dir('RIDGE_TILES', 'the directory of Copernicus GLO-30 tiles (Copernicus_DSM_COG_10_*_DEM.tif)')
W = 6400 # ring columns: 360 deg at 17.78 px/deg
WIN = 1280 # window = 72 deg of the ring
HPX_DEG = W / 360.0
H = 600 # frame rows
COLW = math.radians(360.0 / W)
R_EARTH, KREF = 6371000.0, 0.13
PROM_M, NEAR_CAP, MINW, LINK_D, LINK_A, LOOK = 25.0, 4000.0, 50, 0.45, 0.25, 4
FOREGROUND_M = 5000.0 # a chain nearer than this that is narrower than MINW is kept only where no kept foreground chain covers its columns (the nearest layer must cover every column; a lone stub beside a long near ridge is an arch)
TAPER = 8 # minimum columns over which a chain end dies into its valley (D's 8; round 1 used a fixed 20 AND borrowed the far end's value: the spike)
TAPER_K, TAPER_MAX = 30.0, 72 # a taller end falls over more columns: ceil(30 columns per degree of drop), 8..72: at the authored 1.5-4x exaggeration an end leaves the ridge at about 45 deg or less, not a wall
def taper_len(drop_deg): return np.clip(np.ceil(np.asarray(drop_deg) * TAPER_K), TAPER, TAPER_MAX)
KNOTS = np.array([80, 48, 28, 16, 9, 5, 2.5]) * 1000.0
_r5 = ('#a9b9d3', '#8e9fc2', '#6f86a6', '#556f7d', '#3f5546')
def hx(h): return np.array([int(h[i:i+2], 16) for i in (1, 3, 5)]) / 255.0
def mix(a, b, t): return a + (b - a) * np.asarray(t)[..., None]
RAMP = [mix(hx(_r5[min(int(t), 3)]), hx(_r5[min(int(t), 3) + 1]), t - min(int(t), 3)) for t in np.linspace(0, 4, len(KNOTS))]
GROUND = hx('#4a5f45'); CONTOUR_NICE = (30.0, 50.0, 100.0, 150.0, 200.0, 300.0, 500.0, 1000.0); CONTOUR_MIN_PX = 9.0 # iso-height interval is the smallest nice value that keeps lines >= 9 px apart (p10 of the layer)
ROT = 0 # test hook: frame rotated by ROT columns; noise/azimuth stay tied to the absolute column
P = {} # configured preset
rowang = np.zeros(H); PXDEG_V = 100.0; TOP_DEG = 2.6; BOTTOM_DEG = -3.4; Hc = H
def load_preset(pid, path=Path(__file__).resolve().parent / 'places.json'):
d = json.load(open(path)); p = next(x for x in d['presets'] if x['id'] == pid); p = dict(p); p['date'] = d['date']; return p
# ====================================================================== configuration, DEM
def tiles_for(p):
lat, lon, r = p['lat'], p['lon'], p['rmax_km']; dl = r / 111.32; dn = r / (111.32 * math.cos(math.radians(lat))); out = []
for la in range(math.floor(lat - dl), math.floor(lat + dl) + 1):
for lo in range(math.floor(lon - dn), math.floor(lon + dn) + 1):
cy = min(max(lat, la), la + 1); cx = min(max(lon, lo), lo + 1)
if math.hypot((cy - lat) * 111.32, (cx - lon) * 111.32 * math.cos(math.radians(lat))) <= r: out.append((la, lo))
return out
def tile_path(la, lo): return COP / f"Copernicus_DSM_COG_10_{'N' if la >= 0 else 'S'}{abs(la):02d}_00_{'E' if lo >= 0 else 'W'}{abs(lo):03d}_00_DEM.tif"
def cache_key(*parts):
"""Short digest of everything a cached array was computed from. A cache file is read only under the key of its own inputs, so a changed preset field, tile list or baker constant never meets a stale file."""
h = hashlib.sha256()
for q in parts: h.update(q.tobytes() if isinstance(q, np.ndarray) else repr(q).encode()); h.update(b'|')
return h.hexdigest()[:12]
def mosaic_key(p): return cache_key('mosaic.v1', sorted(tiles_for(p)), 3600)
_mos = {}
def mosaic(p):
"""Tile grid covering the circle (bounding rows/cols), cached as .npy per preset. Missing tiles (open sea has no GLO-30 tile) are zero."""
key = (p['id'], mosaic_key(p))
if key in _mos: return _mos[key]
ts = tiles_for(p); las = [t[0] for t in ts]; los = [t[1] for t in ts]; lat_top = max(las) + 1; lon_left = min(los)
c = B / f"mosaic-{p['id']}-{key[1]}.npy"
if c.exists(): m = np.load(c, mmap_mode='r')
else:
rows = []
for la in range(max(las), min(las) - 1, -1):
rows.append(np.concatenate([tifffile.imread(tile_path(la, lo)) if tile_path(la, lo).exists() else np.zeros((3600, 3600), np.float32) for lo in range(min(los), max(los) + 1)], axis=1))
m = np.concatenate(rows, axis=0).astype(np.float32); np.save(c, m)
_mos[key] = (m, lat_top, lon_left); return _mos[key]
def bil(a, lat_top, lon_left, lat, lon):
y = (lat_top - lat) * 3600 - 0.5; x = (lon - lon_left) * 3600 - 0.5
y0 = np.clip(np.floor(y).astype(int), 0, a.shape[0] - 2); x0 = np.clip(np.floor(x).astype(int), 0, a.shape[1] - 2)
fy = np.clip(y - y0, 0, 1); fx = np.clip(x - x0, 0, 1)
return (a[y0, x0] * (1 - fx) * (1 - fy) + a[y0, x0 + 1] * fx * (1 - fy) + a[y0 + 1, x0] * (1 - fx) * fy + a[y0 + 1, x0 + 1] * fx * fy)
def configure(p):
"""Set the module state for a preset: observer, eye elevation (DEM ground + eye height, or ground_override), window azimuth."""
global P
P = p; m, lt, ll = mosaic(p)
ground = p['ground_override_m'] if 'ground_override_m' in p else float(bil(m, lt, ll, np.array(p['lat']), np.array(p['lon'])))
P['ground_m'] = ground; P['eye_m'] = ground + p['eye_h_m']
if 'look_at' in p:
la, lo = p['look_at']; P['az_center'] = (math.degrees(math.atan2((lo - p['lon']) * math.cos(math.radians(p['lat'])), la - p['lat'])) + 360) % 360
else: P['az_center'] = p['az_deg']
P['win0'] = int(round(P['az_center'] / 360 * W - WIN / 2)) % W
return P
def win(img): return np.take(img, (P['win0'] + np.arange(WIN)) % W, axis=1)
# ====================================================================== authored sub-DEM foreground detail (Aetheria / CultMath ridged multifractal)
def _frac(v): return v - np.floor(v)
def _hash2(x, y): return _frac(np.sin(x * 127.1 + y * 311.7) * 43758.5453) # cultmath_hash(vec2), as Ratatoskr noise.ts
def _vnoise2(x, y): # cultmath_value_noise: 4 hashed corners, smoothstep blend
cx, cy = np.floor(x), np.floor(y); u = x - cx; v = y - cy; u = u * u * (3 - 2 * u); v = v * v * (3 - 2 * v)
a, b, c, d = _hash2(cx, cy), _hash2(cx + 1, cy), _hash2(cx, cy + 1), _hash2(cx + 1, cy + 1)
return (a + (b - a) * u) * (1 - v) + (c + (d - c) * u) * v
LACUN, WEIGHT_GAIN, AMP_GAIN, CRAG = 2.07, 2.0, 0.5, 0.5 # ridgeWorld.ts constants
DETAIL = dict(cell_m=60.0, octaves=4, amp_m=24.0, fade_from_m=500.0, fade_to_m=4000.0, seed=17)
def detail_height(E, N, dist):
"""Ridged multifractal (ridgeWorld.ts `profile`, folded value noise weighted by the previous octave) in WORLD metres (east, north of the observer),
so it is automatically periodic around the ring. First octave wavelength = 60 m = two GLO-30 cells (nothing coarser than the DEM can resolve is
touched), then x2.07. Amplitude fades to zero between fade_from and fade_to metres of distance. Seeded, deterministic, authored."""
cfg = DETAIL; e = 0.03 + (1 - CRAG) * 0.4; top = math.sqrt(1 + e * e); s = 0.0; nrm = 0.0; amp = 1.0; f = 1.0 / cfg['cell_m']; w = 1.0
for k in range(cfg['octaves']):
v = 2 * _vnoise2(E * f + cfg['seed'] * 7.31, N * f + cfg['seed'] * 3.17 + k * 13) - 1
sig = (1 - np.sqrt(v * v + e * e) / top) ** 2 * w; w = np.minimum(1.0, sig * WEIGHT_GAIN); s = s + sig * amp; nrm += amp; amp *= AMP_GAIN; f *= LACUN
t = np.clip((dist - cfg['fade_from_m']) / (cfg['fade_to_m'] - cfg['fade_from_m']), 0, 1); fade = 1 - t * t * (3 - 2 * t)
return cfg['amp_m'] * fade * (s / nrm - 0.3)
# ====================================================================== stage 0: march
def schedule(detail):
d = [500.0]
while d[-1] < P['rmax_km'] * 1000:
d.append(d[-1] + min(250.0, max(10.0 if detail else 40.0, 0.004 * d[-1])))
return np.array(d)
def march(detail=False, curvature=True):
m, lt, ll = mosaic(P); az = np.radians((np.arange(W) + 0.5) / W * 360.0); d = schedule(detail); cosl = math.cos(math.radians(P['lat']))
A = np.zeros((len(d), W), np.float32); Hh = np.zeros_like(A); sa, ca = np.sin(az), np.cos(az)
for k, dk in enumerate(d):
lat = P['lat'] + dk * ca / 111320.0; lon = P['lon'] + dk * sa / (111320.0 * cosl)
h = bil(m, lt, ll, lat, lon)
if detail and dk < DETAIL['fade_to_m']: h = h + detail_height(dk * sa, dk * ca, dk)
drop = dk * dk / (2 * R_EARTH) * (1 - KREF) if curvature else 0.0
A[k] = np.degrees(np.arctan2(h - P['eye_m'] - drop, dk)); Hh[k] = h
return d, A, Hh
def rays_key(detail):
"""Inputs of the ray march: observer, eye elevation, range, the DEM mosaic, the sample schedule (hashed as the array it produces), the Earth / refraction constants, the ring width, and (detail) the foreground detail noise."""
return cache_key('rays.v1', P['lat'], P['lon'], P['rmax_km'], round(float(P['eye_m']), 3), mosaic_key(P), schedule(detail), R_EARTH, KREF, W, bool(detail),
(sorted(DETAIL.items()), LACUN, WEIGHT_GAIN, AMP_GAIN, CRAG) if detail else None)
def stage0(detail=False, force=False):
c = B / f"rays-{P['id']}{'-detail' if detail else ''}-{rays_key(detail)}.npz"
if c.exists() and not force:
z = np.load(c); return z['d'], z['ang'], z['h'], None
t0 = time.time(); mosaic(P); t1 = time.time(); d, A, Hh = march(detail); t2 = time.time(); np.savez(c, d=d, ang=A, h=Hh)
return d, A, Hh, dict(mosaic=t1 - t0, march=t2 - t1, samples=len(d))
# ====================================================================== periodic helpers
def blurw(a, s):
r = int(3 * s) + 1; k = np.exp(-0.5 * (np.arange(-r, r + 1) / s) ** 2); k /= k.sum()
a = np.asarray(a, np.float64); xp = np.concatenate([a[..., -r:], a, a[..., :r]], axis=-1); o = np.zeros_like(a)
for i, kv in enumerate(k): o += kv * xp[..., i:i + a.shape[-1]]
return o
def gradw(a): return (np.roll(a, -1) - np.roll(a, 1)) / 2.0
def _smooth(f): return f * f * (3 - 2 * f)
def pn2(cell, seed):
rng = np.random.default_rng(seed); n = max(int(round(W / cell)), 2); gh = int(H / cell) + 3
g = rng.random((gh, n)).astype(np.float32)
xs = ((np.arange(W) + ROT) % W) * (n / W); x0 = xs.astype(int); fx = _smooth(xs - x0).astype(np.float32); x1 = (x0 + 1) % n
ys = np.arange(H) / cell; y0 = ys.astype(int); fy = _smooth(ys - y0).astype(np.float32)
tp = g[y0][:, x0] * (1 - fx) + g[y0][:, x1] * fx; bt = g[y0 + 1][:, x0] * (1 - fx) + g[y0 + 1][:, x1] * fx
return tp * (1 - fy)[:, None] + bt * fy[:, None]
def pn1(cell, seed): return pn2(cell, seed)[0].astype(np.float64)
def abscol(): return (np.arange(W) + ROT) % W
# ====================================================================== sweep (D), wrapped
NH = 24 # samples of the face history kept per peak (valley -> crest, uniform in elevation angle)
def strip_history(a, dl, hh, k0, k1):
"""The marcher's record of one peak strip: the ray samples from the last local minimum (k0) to the peak (k1), reduced to what the eye sees (a sample hides
behind any earlier, higher one: running maximum of the angle) and resampled at NH uniform steps of elevation angle. Returns (distance[NH], height[NH]), index 0 = valley, NH-1 = crest."""
seg = np.asarray(a[k0:k1 + 1]); runmax = np.maximum.accumulate(seg); vis = np.concatenate([[True], runmax[1:] > runmax[:-1]])
va = runmax[vis]; vd = np.asarray(dl[k0:k1 + 1])[vis]; vh = np.asarray(hh[k0:k1 + 1])[vis]
if len(va) < 2 or va[-1] <= va[0]: return np.full(NH, vd[-1]), np.full(NH, vh[-1])
t = np.linspace(va[0], va[-1], NH); return np.interp(t, va, vd), np.interp(t, va, vh)
def extract(d, ang, hgt, prom=PROM_M):
"""Per column: (peak_ang, valley_ang, peak_d, peak_h, valley_h, valley_d) followed by the face history (NH distances, NH heights) of the strip valley -> peak.
A peak commits when the ray has fallen hy[k] below it and it stands `prom` METRES above the lowest ground the strip has seen (the lowest height, not the height of the lowest-ANGLE sample: that one is often near and high)."""
K, Wc = ang.shape; hy = np.degrees(prom / d).tolist(); dl = d.tolist(); cols = []
for x in range(Wc):
a = ang[:, x].tolist(); hh = hgt[:, x].tolist(); out = []
def span(pk, lo, pd, ph, vh, vd, k0, k1): out.append((pk, lo, pd, ph, vh, vd) + tuple(np.concatenate(strip_history(a, dl, hh, k0, k1))))
runmax = -90.0; mn = a[0]; mnh = hh[0]; mnd = dl[0]; mnk = 0; cand = a[0]; ch = hh[0]; cd = dl[0]; ck = 0; lowh = hh[0]; capped = False # NEAR_CAP is tried once, at the first sample past 4 km
for k in range(1, K):
v = a[k]; dk = dl[k]; hk = hh[k]; lowh = min(lowh, hk)
if dk >= NEAR_CAP and not capped and not out and runmax < -80:
capped = True; j = max(range(k), key=a.__getitem__)
if hh[j] - lowh >= prom: span(a[j], mn, dl[j], hh[j], mnh, mnd, min(mnk, j), j); runmax = a[j]; mn = cand = v; mnh = ch = hk; mnd = cd = dk; mnk = ck = k; lowh = hk
if v < mn: mn = cand = v; mnh = ch = hk; mnd = cd = dk; mnk = ck = k
elif v > cand: cand = v; ch = hk; cd = dk; ck = k
if cand - v >= hy[k] and ch - lowh >= prom:
if cand > runmax: span(cand, mn, cd, ch, mnh, mnd, mnk, ck); runmax = cand
mn = cand = v; mnh = ch = hk; mnd = cd = dk; mnk = ck = k; lowh = hk
if ch - lowh >= prom and cand > runmax: span(cand, mn, cd, ch, mnh, mnd, mnk, ck)
cols.append(out)
return cols
def chains(cols):
Wc = len(cols); parent = []; ids = []
for sp in cols: ids.append(list(range(len(parent), len(parent) + len(sp)))); parent.extend(range(len(parent), len(parent) + len(sp)))
def find(i):
while parent[i] != i: parent[i] = parent[parent[i]]; i = parent[i]
return i
for x, sp in enumerate(cols):
for i, s in zip(ids[x], sp):
p, dd = s[0], s[2]; best = None
for back in range(1, LOOK + 1):
xb = (x - back) % Wc
for j, s2 in zip(ids[xb], cols[xb]):
ld = abs(math.log(dd / s2[2]))
if ld < LINK_D and abs(p - s2[0]) < LINK_A * (1 + 0.5 * (back - 1)):
sc = ld / LINK_D + abs(p - s2[0]) / LINK_A + 0.3 * back
if best is None or sc < best[0]: best = (sc, j)
if best: parent[find(i)] = find(best[1])
out = {}
for x, sp in enumerate(cols):
for i, s in zip(ids[x], sp): out.setdefault(find(i), []).append((x,) + tuple(s))
return list(out.values())
def chain_layer(c):
"""Chain -> dense per-column arrays (deg / m). The chain is unrolled at its largest circular gap into an arc [0, x1] of local columns; every
column outside the arc takes the value of the NEAREST arc end (start end or finish end) and then dies toward that end's valley over TAPER columns."""
a = np.array(c); xs = a[:, 0].astype(int); u = np.unique(xs)
gaps = np.diff(np.concatenate([u, [u[0] + W]])); start = int(u[(gaps.argmax() + 1) % len(u)]); xl = (xs - start) % W
Pk, L, D, PH, VH, VD = (a[:, i] for i in range(1, 7)); HD = a[:, 7:7 + NH]; HH = a[:, 7 + NH:7 + 2 * NH]
run = np.maximum(np.abs(D - VD), 300.0); slope = np.degrees(np.arctan2(np.maximum(PH - VH, 0), run))
x1 = int(xl.max()); ii = np.arange(0, x1 + 1); res = {}; x = np.arange(W)
near_end = np.where(x - x1 <= W - x, x1, 0) # outside the arc: which end is nearer going around the ring
dist = np.where(x <= x1, 0, np.minimum(x - x1, W - x)).astype(float)
def dense(v):
full = np.full(W, np.nan); full[xl] = v; ok = ~np.isnan(full[ii]); full[ii] = np.interp(ii, ii[ok], full[ii][ok]); return np.where(x <= x1, full, full[near_end])
for k, v in dict(top=Pk, lo=L, ph=PH, vh=VH, slope=slope, pd=D, vd=VD).items(): res[k] = dense(v)
res['hd'] = np.stack([dense(HD[:, j]) for j in range(NH)]); res['hh'] = np.stack([dense(HH[:, j]) for j in range(NH)]) # the face history, one (NH, W) array each: column x carries the history of the peak that owns it
drop = res['top'] - np.minimum(res['lo'], res['top'] - 0.4); Tn = taper_len(drop) # end drop in degrees -> taper length at this column
t = np.clip(dist / Tn, 0, 1); t = t * t * (3 - 2 * t)
res['top'] = res['top'] - drop * t
res['wt'] = 1 - np.clip(dist / Tn, 0, 1)
for k in ('top', 'lo', 'ph', 'vh', 'slope', 'wt', 'pd', 'vd', 'hd', 'hh'): res[k] = np.roll(res[k], start, axis=-1)
res['dist'] = float(np.median(D)); res['n'] = len(c); res['x0'] = start; res['x1'] = (start + x1) % W; res['span'] = x1 + 1
res['rise_ang'] = float(np.median(Pk - L)); res['rise_h'] = float(np.median(PH - VH)); res['slope_med'] = float(np.median(slope))
return res
def col_gap(present):
"""Circular distance in columns from each column to the nearest present column (0 where present)."""
x = np.arange(3 * W); p = np.tile(present, 3); last = np.maximum.accumulate(np.where(p, x, -10 * W)); nxt = np.minimum.accumulate(np.where(p, x, 10 * W)[::-1])[::-1]
return np.minimum(x - last, nxt - x)[W:2 * W].astype(float)
def to_slots(layers):
slots = {}
for L in layers: slots.setdefault(int(np.argmin(np.abs(np.log(KNOTS) - math.log(L['dist'])))), []).append(L)
out = []
for s in sorted(slots):
g = slots[s]; top = np.stack([np.where(l['wt'] > 0, l['top'], -1e3 - col_gap(l['wt'] > 0)) for l in g]); j = top.argmax(0) # only chains that are present in a column compete; where none is, the chain whose extent is NEAREST stands (a geometric rule: the list order of chains depends on where column 0 is, and absent-column values feed the blurs at every chain end)
S = dict(slot=s, dist=float(KNOTS[s]), chains=len(g), dists=[l['dist'] for l in g])
for k in ('top', 'lo', 'ph', 'vh', 'slope', 'wt', 'pd', 'vd'): S[k] = np.take_along_axis(np.stack([l[k] for l in g]), j[None], 0)[0]
for k in ('hd', 'hh'): S[k] = np.take_along_axis(np.stack([l[k] for l in g]), j[None, None], 0)[0]
out.append(S)
return out
def sweep(d, ang, hgt, prom=PROM_M):
t0 = time.time(); cols = extract(d, ang, hgt, prom); t1 = time.time()
ch = chains(cols); near = lambda c: np.median([t[3] for t in c]) < FOREGROUND_M; keep = [c for c in ch if len(c) >= MINW]
def arc(c):
u = np.unique([t[0] for t in c]); g = np.diff(np.concatenate([u, [u[0] + W]])); i = g.argmax(); st = u[(i + 1) % len(u)]; return (st + np.arange((u[i] - st) % W + 1)) % W # the columns a chain spans, unrolled at its largest gap
cov = np.zeros(W, bool)
for c in keep:
if near(c): cov[(arc(c)[:, None] + np.arange(-MINW // 4, MINW // 4 + 1)[None, :]).ravel() % W] = True # dilated by MINW/4 columns: a stub that occludes a near ridge for a few columns sits inside it
for c in sorted((c for c in ch if len(c) < MINW and near(c)), key=lambda c: -len(c)): # widest first: a stub is judged against the stubs already kept
xs = arc(c)
if not cov[xs].any(): keep.append(c); cov[xs] = True
lay = [chain_layer(c) for c in keep]; slots = to_slots(lay); t2 = time.time()
return dict(cols=cols, nchains=len(ch), keep=lay, slots=slots, raw=sum(len(c) for c in cols), t_extract=t1 - t0, t_chain=t2 - t1)
# ====================================================================== the authored Y mapping
def set_view(slots, sky_row=0.04, near_row=0.66):
"""AUTHORED vertical mapping, linear in elevation angle: row = (TOP_DEG - angle) * PXDEG_V over the whole 600-row frame (cylindrical, horizontal 17.78 px/deg).
Two anchors: the highest ridge top in the window sits at `sky_row` (4%) of the frame height, the LOWEST crest of the nearest layer in the window sits at `near_row` (66%);
everything below is that layer's body, its foot, and the ground-plane band. PXDEG_V = (near_row - sky_row) * 600 / (highest top - lowest near crest), clamped to
[17.78, 178] px/deg = vertical exaggeration 1x..10x. TOP_DEG = highest top + sky_row*600/PXDEG_V; BOTTOM_DEG = TOP_DEG - 600/PXDEG_V; horizon (angle 0) at row TOP_DEG*PXDEG_V."""
global rowang, PXDEG_V, TOP_DEG, BOTTOM_DEG
cols = (P['win0'] + np.arange(WIN)) % W; sky = max(float(S['top'][cols][S['wt'][cols] > 0].max()) for S in slots if (S['wt'][cols] > 0).any())
near = slots[-1]; act = near['wt'][cols] > 0; c_low = float(blurw(near['top'], 3.0)[cols][act].min()) if act.any() else sky - 2.0
ppd = float(np.clip((near_row - sky_row) * H / max(sky - c_low, 0.05), HPX_DEG, 10 * HPX_DEG))
PXDEG_V = ppd; TOP_DEG = sky + sky_row * H / ppd; BOTTOM_DEG = TOP_DEG - H / ppd; rowang = TOP_DEG - (np.arange(H) + 0.5) / PXDEG_V
return dict(top_deg=TOP_DEG, bottom_deg=BOTTOM_DEG, px_per_deg=PXDEG_V, exaggeration=PXDEG_V / HPX_DEG, horizon_row=TOP_DEG * PXDEG_V, sky_deg=sky, near_crest_low_deg=c_low,
sky_row=sky_row, near_row=near_row)
# ====================================================================== bake (sun-independent), per layer
def azrad(): return (abscol() + 0.5) / W * 2 * math.pi
def face_field(hd, hh, ang, lo, pd, ph, vd, vh):
"""Per-pixel face from the marcher's history (hd, hh: NH distances / heights per column, valley -> crest). A pixel at elevation `rowang` of column x sits at fraction
f = (rowang - lo) / (ang - lo) of the strip; its distance and height are the history interpolated at f. Front slope s = dh/dd along the ray comes from the history over +-2 nodes
(the strip's mean slope where the window has no run); the lateral tilt is the height change between the neighbouring columns at the SAME pixel row, corrected for the radial
step (c - a*s) / (tangential length). n ~ (-s r - la t, 1). Returns Z, D, s, la, f as (H, W) arrays."""
f = np.clip((rowang[:, None] - lo[None, :]) / np.maximum(ang - lo, 1e-3)[None, :], 0, 1); pos = f * (NH - 1); i0 = np.minimum(pos.astype(int), NH - 2); t = pos - i0; xs = np.arange(W)[None, :]
lerp = lambda A: A[i0, xs] * (1 - t) + A[i0 + 1, xs] * t
lo_i, hi_i = np.clip(np.arange(NH) - 2, 0, NH - 1), np.clip(np.arange(NH) + 2, 0, NH - 1)
dd = hd[hi_i] - hd[lo_i]; mean_s = ((ph - vh) / np.maximum(np.abs(pd - vd), 300.0))[None, :]
sn = np.where(dd >= 20.0, (hh[hi_i] - hh[lo_i]) / np.maximum(dd, 20.0), mean_s)
Z, Dp, sl = lerp(hh), lerp(hd), lerp(sn)
a_ = np.roll(Dp, -1, axis=1) - np.roll(Dp, 1, axis=1); c_ = np.roll(Z, -1, axis=1) - np.roll(Z, 1, axis=1)
b_ = np.maximum(2 * COLW * Dp, 1.0); la = np.where(np.abs(a_) > 3 * b_, 0.0, np.clip((c_ - a_ * sl) / b_, -1.5, 1.5)) # neighbours at very different distances are an occlusion edge, not a surface: no tilt
return Z.astype(np.float32), Dp, blurw(sl, 2.0).astype(np.float32), blurw(la, 2.0).astype(np.float32), f
def bake_layer(S, rank, far):
b = S['slot']; base = RAMP[b]; reg = int(round((rank - 2) * 1.5)); R_ = lambda v: np.roll(v, reg, axis=-1); px = 1.0 / PXDEG_V
s_ref = P.get('s_ref_deg', 8.0)
ang = blurw(S['top'], 2.5) + (pn1(38, 100 + b) - 0.5) * 6 * px + (pn1(7, 200 + b) - 0.5) * 2 * px; ang = R_(ang) # edge wobble authored in pixels
lo = R_(blurw(S['lo'], 2.5)); ph = R_(blurw(S['ph'], 2.0)); vh = R_(blurw(S['vh'], 2.0))
sn = np.clip(R_(blurw(S['slope'], 3.0)) / s_ref, 0, 1); wt = R_(S['wt'])
hd = R_(blurw(S['hd'], 2.0)); hhist = R_(blurw(S['hh'], 2.0)); pd_ = R_(blurw(S['pd'], 2.0)); vd_ = R_(blurw(S['vd'], 2.0))
top = ang[None, :] - rowang[:, None]; toppx = (top * PXDEG_V).astype(np.float32)
m = np.clip(toppx / 1.6 + 0.5, 0, 1) * np.clip(blurw((wt > 0).astype(float), 1.0) * 1.5, 0, 1)[None, :] # a layer is drawn only where a chain exists (or tapers)
rise_px = np.clip((ang - blurw(lo, 12.0)) * PXDEG_V, 1, None) # foot length from a smoothed valley: no per-column curtains
foot = np.clip(0.85 * rise_px, 28 * np.clip(wt, 0, 1), 210).clip(1) # a tapering end carries no 28 px curtain: the body shrinks with the chain's weight, so a chain end leaves no stub
depthpx = np.clip(toppx / foot[None], 0, 1)
fade = np.clip((depthpx - 0.6) / 0.4, 0, 1); fade = fade * fade * (3 - 2 * fade); alpha = (m * (1 - fade)).astype(np.float32) # the body dissolves at its foot: what is behind (other layers, floor, sky) shows, so a column boundary leaves no wall down the frame
w = (depthpx ** 0.8 * 0.92).astype(np.float32) # weight of the foot colour (floor/haze mix, decided at relight)
face = np.broadcast_to((base * 1.12).astype(np.float32), (H, W, 3)).copy()
Z, Dp, sl, la, f = face_field(hd, hhist, ang, lo, pd_, ph, vd_, vh)
snp = np.clip(np.degrees(np.arctan(sl)) / s_ref, 0, 1) # per-pixel face steepness from the history (the column scalar `sn` stays for stroke weights only)
face *= (1.07 - 0.22 * snp)[..., None]
mot = 0.55 * pn2(22, 300 + b) + 0.45 * pn2(5, 400 + b)
face *= (0.90 + 0.18 * mot[..., None]) * (0.96 + 0.08 * pn2(1, 500 + b)[..., None])
echo = np.clip(((ang - 24 * px * (1 + rank * .3))[None, :] - rowang[:, None]) * PXDEG_V / 1.6 + .5, 0, 1)
face *= (1 - 0.13 * echo[..., None] * (1 - depthpx[..., None] * .5))
ink = base * 0.42 + np.array([0.02, 0.03, 0.08])
ew = (0.7 + 1.8 * snp) * (1 - 0.35 * far)
edge = np.clip(1 - np.abs(toppx - 0.6 - ew * 0.5) / ew, 0, 1) * (toppx > -0.5)
ea = 0.55 * (1 - 0.45 * far) * edge * (0.4 + 0.6 * snp)
dzdpx = np.maximum(np.abs(np.roll(Z, -1, axis=0) - np.roll(Z, 1, axis=0)) / 2.0, 1e-3); dzdpx[0] = dzdpx[1]; dzdpx[-1] = dzdpx[-2] # metres of true height per pixel row, per pixel
faceok = (f > 0.02) & (f < 0.98) & (wt[None, :] > 0.5)
cm = next((c for c in CONTOUR_NICE if c >= CONTOUR_MIN_PX * float(np.percentile(dzdpx[faceok], 75) if faceok.any() else 3.0)), CONTOUR_NICE[-1]) # metres per contour line: lines >= 9 px apart over three quarters of this layer's face
p_ = (Z / cm) % 1.0; dpx = (np.minimum(p_, 1 - p_) * cm / dzdpx).astype(np.float32)
face_mask = ((top > 0.003) & (rowang[:, None] > lo[None, :])).astype(np.float32)
ca = (0.95 * (1 - 0.45 * far) * (1 - depthpx) ** 2.0 * wt[None, :] * (0.6 + 0.4 * snp) * face_mask).astype(np.float32)
hatch = np.zeros((H, W), np.float32)
if far < 0.8:
rng = np.random.default_rng(700 + b); xs = np.cumsum(rng.integers(2, 7, W // 2)); xs = xs[xs < W]
up = gradw(ang) * PXDEG_V
for i, xa in enumerate(xs):
x = (xa - ROT) % W; u1, u2 = rng.random(2)
if u1 > 0.30 + 0.70 * sn[x] or wt[x] < 0.5: continue
Lpx = (8 + 50 * sn[x] ** 0.7) * (0.55 + 0.7 * u2) * (1 - 0.5 * far)
y0 = (rowang[0] - ang[x]) * PXDEG_V; yy = np.arange(H); t = (yy - y0 - 2) / max(Lpx, 1)
ok = (t >= 0) & (t < 1); xx = (x + np.round(-np.clip(up[x] * 0.25, -0.6, 0.6) * (yy - y0)).astype(int)) % W
hatch[yy[ok], xx[ok]] = np.maximum(hatch[yy[ok], xx[ok]], (1 - t[ok]) ** 1.5)
ha = np.clip(ea + hatch * 0.50 * (1 - 0.4 * far), 0, 1).astype(np.float32)
rb = (np.clip(1 - toppx / 7.0, 0, 1) * (toppx > -0.5)).astype(np.float32)
az = azrad(); dd = np.stack([S['vd'], 0.5 * (S['vd'] + S['pd']), S['pd']]); hh = np.stack([S['vh'], 0.5 * (S['vh'] + S['ph']), S['ph']]) + 3.0 # face samples: valley, middle, crest
spts = dict(lat=np.roll(P['lat'] + dd * np.cos(az) / 111320.0, reg, axis=1), lon=np.roll(P['lon'] + dd * np.sin(az) / (111320.0 * math.cos(math.radians(P['lat']))), reg, axis=1), h=np.roll(hh, reg, axis=1))
return dict(m=m.astype(np.float32), a=alpha, face=face, w=w, dpx=dpx, ca=ca, ha=ha, rb=rb, mot=mot.astype(np.float32), sl=sl, la=la, ink=ink, far=far, dist=S['dist'], rank=rank, ang=ang, lo=lo, spts=spts)
def bake(slots):
return [bake_layer(S, rank, float(np.clip((math.log(S['dist']) - math.log(2000)) / (math.log(70000) - math.log(2000)), 0, 1))) for rank, S in enumerate(slots)]
# ====================================================================== cast shadows (sun-dependent): DEM march toward the sun from each layer's face samples
SH_STEPS = 40.0 * 1.1 ** np.arange(0, 66) # 40 m .. 20 km, geometric
def cast_shadows(layers, az_s, el_s):
"""Per layer a (3, W) array in [0,1]: how much of the sun is blocked by the real DEM at that layer's valley / middle / crest sample. One vectorised march for all
layers. Soft: the blocked fraction is the terrain excess over the sun ray, over a penumbra of 6 m + 3% of the distance travelled. Sun at or below 0.3 deg: no direct light, no shadow."""
if el_s <= 0.3: return [np.zeros((3, W)) for _ in layers]
m, lt, ll = mosaic(P); lat = np.concatenate([Ly['spts']['lat'].ravel() for Ly in layers]); lon = np.concatenate([Ly['spts']['lon'].ravel() for Ly in layers]); h = np.concatenate([Ly['spts']['h'].ravel() for Ly in layers])
t = math.tan(math.radians(el_s)); cosl = math.cos(math.radians(P['lat'])); ca, sa = math.cos(math.radians(az_s)), math.sin(math.radians(az_s)); sh = np.zeros_like(h)
for s in SH_STEPS:
terr = bil(m, lt, ll, lat + s * ca / 111320.0, lon + s * sa / (111320.0 * cosl)) - s * s / (2 * R_EARTH)
sh = np.maximum(sh, np.clip((terr - (h + s * t)) / (6.0 + 0.03 * s), 0, 1))
return [a.reshape(3, W) for a in np.split(sh, len(layers))]
# ====================================================================== relight (sun-dependent)
TAU_Z = np.array([0.09, 0.14, 0.30]); BETA_R = np.array([0.0058, 0.0135, 0.0331]); BETA_M = 0.023; ATM = 0.15
ATM_FLOOR = 1.2 # authored: valley haze over the ground plane is heavier than the layer extinction ATM, so the floor thins into haze gradually
G_MIE, W_MIE = 0.7, 0.10; REF = dict(az=180.0, el=62.0); REF_AMB = np.array([0.62, 0.74, 1.0])
def airmass(el): e = max(el, 0.3); return 1.0 / (math.sin(math.radians(e)) + 0.50572 * (e + 6.07995) ** -1.6364)
def sun_colour(el): return np.exp(-TAU_Z * airmass(el))
def phase(c):
pr = 3.0 / (16 * math.pi) * (1 + c * c); pm = (1 - G_MIE ** 2) / (4 * math.pi * (1 + G_MIE ** 2 - 2 * G_MIE * c) ** 1.5)
return (1 - W_MIE) * pr + W_MIE * pm
def sun_vec(az, el): a, e = math.radians(az), math.radians(el); return np.array([math.sin(a) * math.cos(e), math.cos(a) * math.cos(e), math.sin(e)])
def EX_REF_NORM(Ec, I, amb_k): return float((0.45 * REF_AMB * amb_k + 0.75 * Ec * I * 0.7).mean()) / 0.87
def sun_state(az, el):
Ec = sun_colour(el) / sun_colour(REF['el']); I = (max(math.sin(math.radians(el)), 0.0) / math.sin(math.radians(REF['el']))) ** 0.5
tw = 0.30 * float(np.clip((el + 6) / 8, 0, 1)) ** 2; amb_k = 0.12 + 0.88 * float(np.clip((el + 3) / 45, 0, 1)) ** 0.7
ex = 0.45 + 0.55 * float(np.clip((el + 4) / 18, 0, 1)) ** 0.6
return dict(az=az, el=el, L=sun_vec(az, el), sun=Ec * I, sky_sun=Ec * (I + tw), amb=REF_AMB * amb_k, amb_k=amb_k, ex=ex / EX_REF_NORM(Ec, I, amb_k), ex_t=ex)
def _g0():
st = sun_state(REF['az'], REF['el']); d = 0.7; g = (0.45 * st['amb'] + 0.75 * st['sun'] * d) * (np.array([0.90, 0.92, 1.04]) + (np.array([1.03, 1.0, 0.95]) - np.array([0.90, 0.92, 1.04])) * d)
return g / g.mean()
G0 = _g0()
def col_scatter(st):
az = azrad(); c = math.cos(math.radians(st['el'])) * np.cos(az - math.radians(st['az']))
pr = np.array([phase(float(x)) for x in c]); g = np.clip((pr / phase(0.0)) ** 0.35, 0.5, 2.4)
rs = sun_state(REF['az'], REF['el']); ref = 0.8 * rs['sky_sun'] + 0.35 * rs['amb']; ref = ref / ref.mean()
col = 0.8 * st['sky_sun'] + 0.35 * st['amb']; col = col / col.mean()
return c, g[:, None] * (col / ref)[None, :]
def hue_gain(st, d):
lit_tint = np.array([1.03, 1.0, 0.95]); shd_tint = np.array([0.90, 0.92, 1.04])
tint = shd_tint + (lit_tint - shd_tint) * d
return (0.45 * st['amb'] + 0.75 * st['sun'] * d) * tint / G0 * st['ex']
def floor_img(st, hazec, fmot):
"""Ground-plane band: below the horizon a flat floor `floor_agl_m` below the eye; plane distance dg = h / tan(depression); colour = ground albedo lit by an
up-facing normal, thinned into haze by the same Beer-Lambert extinction as the layers. Above the horizon the same array is the haze (sky covers it)."""
a = rowang; dg = P['floor_agl_m'] / np.tan(np.radians(np.clip(-a, 0.002, 89.0)))
T = np.exp(-(BETA_R + BETA_M)[None, :] * (dg / 1000.0)[:, None] * ATM_FLOOR); T[a >= 0] = 0.0
ground = (GROUND * hue_gain(st, float(np.clip((math.sin(math.radians(st['el'])) + 0.2) / 1.2, 0, 1))))
return ground[None, None, :] * fmot[..., None] * T[:, None, :] + hazec[None, :, :] * (1 - T)[:, None, :]
def ndl_field(Ly, L):
"""n . L per pixel from the baked face: n ~ (-sl r - la t, 1), r / t the column's radial / tangential directions. Returns (ndl, |n|-normaliser)."""
az = azrad(); sa, ca_ = np.sin(az), np.cos(az); sl, la = Ly['sl'], Ly['la']
Lr = (L[0] * sa + L[1] * ca_)[None, :]; Lt = (L[0] * ca_ - L[1] * sa)[None, :]; nrm = np.sqrt(1 + sl * sl + la * la)
return (L[2] - sl * Lr - la * Lt) / nrm, nrm
RELIEF_OFF_DEG, RELIEF_EL = 45.0, 32.0 # authored relief light: RELIEF_OFF_DEG off the view axis, RELIEF_EL up
def win_centre_az(): return (P['win0'] + WIN / 2) / W * 360.0
def relief_vec(sun_az):
"""The authored relief light direction (toward the light). The real sun lights a face the viewer looks at only when it is off to one side; with the sun behind the viewer every face is
head-on and the mountain goes flat. Shading n.L (ndl_field, rim, spec) therefore uses this fixed direction: azimuth = the view centre az +/- RELIEF_OFF_DEG, on the left/right side of the view axis
that the real sun is on (so the two agree where they can), RELIEF_EL up. Sun colour/strength, haze, colour temperature and cast shadows keep the real sun:
the relief only decides where on the form the sun's light falls."""
ctr = win_centre_az(); side = 1.0 if math.sin(math.radians(sun_az - ctr)) >= 0 else -1.0
return sun_vec(ctr + side * RELIEF_OFF_DEG, RELIEF_EL)
def fresnel(sl, nrm): return np.clip(1 - sl / nrm, 0, 1) ** 3 # (1 - n . V)^3, V = -r: strongest on faces that lie away from the viewer
def spec_field(L, sl, la, nrm, ndl):
"""GGX-like highlight (alpha 0.3) of the face (sl, la, nrm as in ndl_field) for the light L: half vector between L and the direction back to the viewer (-r) in every column."""
az = azrad(); sa, ca_ = np.sin(az), np.cos(az)
Hv = L[None, :] + np.stack([-sa, -ca_, np.zeros(W)], 1); Hv /= np.linalg.norm(Hv, axis=1, keepdims=True)
Hr = (Hv[:, 0] * sa + Hv[:, 1] * ca_)[None, :]; Ht = (Hv[:, 0] * ca_ - Hv[:, 1] * sa)[None, :]
nh = np.clip((Hv[:, 2][None, :] - sl * Hr - la * Ht) / nrm, 0, 1); al = 0.30
D = (al * al) / (1 + nh * nh * (al * al - 1)) ** 2 * al * al
return 0.30 * D * np.clip(ndl, 0, 1)
def relight_layer(Ly, st, hazec, F, sh3=None):
"""Shading per PIXEL (ndl_field), not per column. n.L uses the relief light; the sun's colour, strength (st['sun'] scales every direct term, so night and dusk stay dark) and shadows are the real sun's."""
L = relief_vec(st['az']); az = azrad(); sa, ca_ = np.sin(az), np.cos(az); sl, la = Ly['sl'], Ly['la']
ndl, nrm = ndl_field(Ly, L); d = np.clip((ndl + 0.2) / 1.2, 0, 1)
gain = hue_gain(st, d[..., None])
rim_dir = (np.stack([sa, ca_], 1) * L[None, :2]).sum(1)[None, :]
fres = fresnel(sl, nrm)
rim = (0.12 + 1.0 * np.clip(rim_dir, 0, 1)) * fres
spec = spec_field(L, sl, la, nrm, ndl)
if sh3 is None: sp = np.zeros((1, 1), np.float32); face = Ly['face'] * gain.astype(np.float32)
else:
sh3 = blurw(sh3, 2.0); f = np.clip((rowang[:, None] - Ly['lo'][None, :]) / np.maximum(Ly['ang'] - Ly['lo'], 1e-3)[None, :], 0, 1) # position on the face: 0 valley .. 1 crest
sp = np.where(f < 0.5, sh3[0] + (sh3[1] - sh3[0]) * 2 * f, sh3[1] + (sh3[2] - sh3[1]) * 2 * (f - 0.5)).astype(np.float32)
g_amb = hue_gain(st, 0.0)[None, None, :]; face = Ly['face'] * (g_amb + (gain - g_amb) * (1 - sp)[..., None]).astype(np.float32) # shadow removes the direct share only
fa = 1 - math.exp(-Ly['dist'] / 20000.0) # valley mist grows with the layer's distance
foot = F * (1 - fa) + hazec[None] * fa
wf = Ly['w'] * (1 - 0.6 * fa) # far layers keep 60% of their own body colour at the foot
wf = 1 - (1 - wf) * np.clip(1 + rowang * PXDEG_V / 40.0, 0, 1)[:, None] # below the horizon every layer's foot has merged into the ground plane within 40 px: no per-slot colour seams down the floor
face = face * (1 - wf[..., None]) + foot * wf[..., None]
face += (Ly['rb'] * (rim * 0.55) * (1 - sp))[..., None] * (st['sun'] * 1.1 + 0.25)[None, None, :]
face += ((spec * (0.5 + Ly['mot'])) * (Ly['m'] > 0.5) * (1 - sp))[..., None] * st['sun'][None, None, :] * 0.5
lw = (0.70 + 0.90 * (1 - d)) * (1 - 0.30 * Ly['far'])
line = np.clip(1 - Ly['dpx'] / lw, 0, 1) * Ly['ca']
ia = np.clip(line + Ly['ha'], 0, 1)
ink = Ly['ink'][None, None, :] * (0.6 + 0.4 * gain.mean(-1))[..., None]
face = face * (1 - ia[..., None]) + ink * ia[..., None]
T = np.exp(-(BETA_R + BETA_M) * Ly['dist'] / 1000.0 * ATM)
return face * T[None, None, :] + hazec[None, :, :] * (1 - T)[None, None, :]
def sky_hz(st, paperf):
az = azrad(); c, hg = col_scatter(st)
hz = hx('#f5e2c8')[None, :] * hg * st['ex_t']
tint = (st['amb'] / st['amb'].mean()) / (REF_AMB / REF_AMB.mean()); zen = hx('#a9bcd6')[None, :] * tint[None, :] * st['ex_t']
u = np.clip(rowang / max(TOP_DEG, 2.5), 0, 1)[:, None] ** 0.7
img = zen[None] * u[..., None] + hz[None] * (1 - u)[..., None]
ev = np.radians(rowang)[:, None]; cs = np.cos(ev) * math.cos(math.radians(st['el'])) * np.cos(az - math.radians(st['az']))[None, :] + np.sin(ev) * math.sin(math.radians(st['el']))
ang = np.degrees(np.arccos(np.clip(cs, -1, 1)))
glow = 0.95 * np.exp(-(ang / 4.5) ** 2) + 0.55 * np.exp(-(ang / 22.0) ** 2)
img = img + glow[..., None] * (st['sky_sun'] * 0.8 * st['ex_t'])[None, None, :]
return img * paperf[..., None] / 0.93, hz
def relight(layers, az, el, paperf, grain, fmot, shadows=True):
st = sun_state(az, el); shs = cast_shadows(layers, az, el) if shadows else [None] * len(layers); sky, hz = sky_hz(st, paperf)
hazec = hz * (hx('#cdd7e7') / hx('#f5e2c8'))[None, :]
F = floor_img(st, hazec, fmot) * (paperf[..., None] / 0.93)
img = np.where((rowang >= 0)[:, None, None], sky, F).astype(np.float32)
for Ly, sh3 in zip(layers, shs):
face = relight_layer(Ly, st, hazec.astype(np.float32), F.astype(np.float32), sh3); face *= (paperf[..., None] / 0.93)
img = img * (1 - Ly['a'][..., None]) + face * Ly['a'][..., None]
return img * grain[..., None]
def statics():
paperf = 0.93 + 0.09 * (0.5 * pn2(2, 11) + 0.3 * pn2(9, 12) + 0.2 * pn2(40, 13)) + 0.03 * (pn2(1, 14) - .5)
return paperf, 0.95 + 0.07 * pn2(1, 900), 0.88 + 0.2 * (0.6 * pn2(30, 810) + 0.4 * pn2(7, 811))
def debug_spans(cols, win0=None, w=WIN):
win0 = P['win0'] if win0 is None else win0
img = np.ones((H, w, 3)) * 0.97; ramp = [hx(c) for c in ('#d03030', '#e8a020', '#60b040', '#2090b0', '#6030c0')]
order = sorted(((s[2], i, s[0], s[1]) for i in range(w) for s in cols[(win0 + i) % W]), key=lambda t: -t[0])
for dd, x, p, l in order:
t = np.clip((math.log(dd) - math.log(500)) / (math.log(95000) - math.log(500)), 0, 1) * 4
i = min(int(t), 3); col = mix(ramp[i], ramp[i + 1], t - i)
y0 = int(np.clip((TOP_DEG - p) * PXDEG_V, 0, H - 1)); y1 = int(np.clip((TOP_DEG - max(l, p - 3.0)) * PXDEG_V, 0, H - 1))
img[y0:y1 + 1, x] = col
return img
def apply_view(v):
"""re-apply a view returned by set_view (used by the rotation test and the before/after comparison so that frames line up exactly)"""
global rowang, PXDEG_V, TOP_DEG, BOTTOM_DEG
PXDEG_V = v['px_per_deg']; TOP_DEG = v['top_deg']; BOTTOM_DEG = v['bottom_deg']; rowang = TOP_DEG - (np.arange(H) + 0.5) / PXDEG_V

29
scripts/ridge-bake/sun.py Normal file
View file

@ -0,0 +1,29 @@
"""Solar position: NOAA general solar position formulas (Spencer/fractional-year low precision: equation of time and declination Fourier series,
hour angle, zenith/azimuth), https://gml.noaa.gov/grad/solcalc/solareqns.PDF . Accuracy ~0.1 deg, plenty for a relight. Apparent elevation adds
Saemundsson refraction (Meeus ch.16). UTC in, degrees out. Azimuth clockwise from north."""
import math
from datetime import datetime, timezone, timedelta
def solar(dt_utc, lat, lon):
doy = dt_utc.timetuple().tm_yday; hr = dt_utc.hour + dt_utc.minute / 60 + dt_utc.second / 3600
g = 2 * math.pi / 365 * (doy - 1 + (hr - 12) / 24)
eq = 229.18 * (0.000075 + 0.001868 * math.cos(g) - 0.032077 * math.sin(g) - 0.014615 * math.cos(2 * g) - 0.040849 * math.sin(2 * g))
de = (0.006918 - 0.399912 * math.cos(g) + 0.070257 * math.sin(g) - 0.006758 * math.cos(2 * g) + 0.000907 * math.sin(2 * g)
- 0.002697 * math.cos(3 * g) + 0.00148 * math.sin(3 * g))
tst = hr * 60 + eq + 4 * lon # true solar time, minutes
ha = math.radians(tst / 4 - 180); la = math.radians(lat)
sin_el = math.sin(la) * math.sin(de) + math.cos(la) * math.cos(de) * math.cos(ha)
el = math.degrees(math.asin(sin_el))
az = (math.degrees(math.atan2(math.sin(ha), math.cos(ha) * math.sin(la) - math.tan(de) * math.cos(la))) + 180) % 360
ref = 1.02 / math.tan(math.radians(el + 10.3 / (el + 5.11))) / 60 if el > -5 else 0.0
return az, el + ref, el, math.degrees(de), eq
def find_el(day, lat, lon, target, lo_h, hi_h):
"""bisect the UTC hour in [lo_h, hi_h] where apparent elevation crosses `target`."""
f = lambda h: solar(day + timedelta(hours=h), lat, lon)[1] - target
a, b = lo_h, hi_h; fa = f(a)
for _ in range(50):
m = (a + b) / 2; fm = f(m)
if (fa < 0) == (fm < 0): a, fa = m, fm
else: b = m
return (a + b) / 2

65
scripts/ridge-bake/tests/mutate.py Normal file
View file

@ -0,0 +1,65 @@
"""Mutation check: re-introduce a cause into a COPY of the baker, run the fixtures against it, and require that at least one fails.
venv/bin/python tests/mutate.py (prints one line per mutant; exit 1 if a mutant survives)"""
import os, shutil, subprocess, sys, tempfile
HERE = os.path.dirname(os.path.abspath(__file__)); SRC = os.path.join(HERE, '..')
MUTANTS = {
'round-1 bug: both chain sides take the END value (near_end always x1)': ("near_end = np.where(x - x1 <= W - x, x1, 0)", "near_end = np.where(x - x1 <= W - x, x1, x1)"),
'linking does not wrap (look-back columns clipped at 0)': ("xb = (x - back) % Wc", "xb = max(x - back, 0)"),
'a chain absent from a column still competes in the slot union (held valley value wins)': ("np.where(l['wt'] > 0, l['top'], -1e3 - col_gap(l['wt'] > 0))", "l['top']"),
'shadow march goes AWAY from the sun': ("lat + s * ca / 111320.0, lon + s * sa", "lat - s * ca / 111320.0, lon - s * sa"),
'shadow ignores the sun elevation (ray stays at the sample height)': ("(h + s * t)", "(h + 0 * t)"),
'taper length ignores the drop (always the 8-column minimum)': ("Tn = taper_len(drop)", "Tn = np.full_like(drop, 8.0)"),
'shadow march mirrors east/west (longitude sign)': ("lon + s * sa / (111320.0 * cosl)", "lon - s * sa / (111320.0 * cosl)"),
'shadow longitude scale ignores cos(lat)': ("lon + s * sa / (111320.0 * cosl)", "lon + s * sa / 111320.0"),
'shadow lightens instead of darkens the face': ("(gain - g_amb) * (1 - sp)[..., None]", "(gain - g_amb) * sp[..., None]"),
'taper cap lowered to 8 columns': ("TAPER_K, TAPER_MAX = 30.0, 72", "TAPER_K, TAPER_MAX = 30.0, 8"),
'taper per degree lowered to 2': ("TAPER_K, TAPER_MAX = 30.0, 72", "TAPER_K, TAPER_MAX = 2.0, 72"),
'taper floor ignores the chain end (top held, no fall)': ("res['top'] = res['top'] - drop * t", "res['top'] = res['top']"),
'painter shades one normal per column (the column mean of the per-pixel face)': ("sl, la = Ly['sl'], Ly['la']\n Lr", "sl, la = np.broadcast_to(Ly['sl'].mean(0), Ly['sl'].shape), np.broadcast_to(Ly['la'].mean(0), Ly['la'].shape)\n Lr"),
'history collapsed to ONE slope per peak (the strip mean)': ("Z, Dp, sl = lerp(hh), lerp(hd), lerp(sn)", "Z, Dp, sl = lerp(hh), lerp(hd), np.broadcast_to(mean_s, Z.shape)"),
'history not recorded: only valley and peak samples kept': ("return np.interp(t, va, vd), np.interp(t, va, vh)", "return np.interp(t, va[[0, -1]], vd[[0, -1]]), np.interp(t, va[[0, -1]], vh[[0, -1]])"),
'no lateral tilt from the neighbouring columns': ("np.clip((c_ - a_ * sl) / b_, -1.5, 1.5))", "0 * np.clip((c_ - a_ * sl) / b_, -1.5, 1.5))"),
'NEAR_CAP shoulder committed without the prominence test': ("if hh[j] - lowh >= prom: span(", "if True: span("),
'prominence measured as an angle at the CURRENT sample (round 2)': ("if cand - v >= hy[k] and ch - lowh >= prom:", "if cand - v >= hy[k] and cand - mn >= hy[k]:"),
'contour interval chosen from the shallowest quarter of the face (banding)': ("np.percentile(dzdpx[faceok], 75)", "np.percentile(dzdpx[faceok], 10)"),
'a tapering end keeps the 28 px minimum body (stub / curtain at every chain end)': ("28 * np.clip(wt, 0, 1)", "28"),
'absent columns of a slot take the first chain in list order (rotation-dependent)': ("np.where(l['wt'] > 0, l['top'], -1e3 - col_gap(l['wt'] > 0))", "np.where(l['wt'] > 0, l['top'], -1e9)"),
'history keeps samples the eye cannot see (no running maximum)': ("runmax = np.maximum.accumulate(seg)", "runmax = seg"),
'lateral tilt scale doubled': ("/ b_, -1.5, 1.5))", "/ b_ * 2, -1.5, 1.5))"),
'lateral tilt kept across occlusion edges': ("np.abs(a_) > 3 * b_", "np.abs(a_) > 3e9 * b_"),
'every foreground stub kept (round 2 arches)': ("if not cov[xs].any(): keep.append(c)", "if True: keep.append(c)"),
'no foreground stub kept (holes in the nearest layer)': ("if not cov[xs].any(): keep.append(c)", "if False: keep.append(c)"),
'prominence from the height of the lowest-angle sample (near, high valley)': ("if cand - v >= hy[k] and ch - lowh >= prom:", "if cand - v >= hy[k] and ch - mnh >= prom:"),
'layer body painted to the bottom of the frame (no foot fade: walls)': ("alpha = (m * (1 - fade)).astype(np.float32)", "alpha = m.astype(np.float32)"),
'relief: shading uses the real sun again (face flat with the sun behind the viewer)': ("L = relief_vec(st['az']); az", "L = st['L']; az"),
'relief: 45 deg off the viewer-ward axis (the previous reading) instead of the view axis': ("ctr = win_centre_az(); side", "ctr = (win_centre_az() + 180.0) % 360.0; side"),
'relief: always the + side, whatever side the real sun is on': ("side = 1.0 if math.sin(math.radians(sun_az - ctr)) >= 0 else -1.0", "side = 1.0"),
'relief: colour / strength / haze frozen at a 40 deg sun (night does not go dark)': ("st = sun_state(az, el); shs", "st = sun_state(az, 40.0); shs"),
'relief: cast shadows march toward the relief light instead of the real sun': ("cast_shadows(layers, az, el) if shadows", "cast_shadows(layers, math.degrees(math.atan2(relief_vec(az)[0], relief_vec(az)[1])) % 360.0, RELIEF_EL) if shadows"),
'compositor (relight) reads the layer mask m instead of its alpha a (walls return)': ("img = img * (1 - Ly['a'][..., None]) + face * Ly['a'][..., None]", "img = img * (1 - Ly['m'][..., None]) + face * Ly['m'][..., None]"),
'foot fade starts at 30% of the foot instead of 60%': ("(depthpx - 0.6) / 0.4", "(depthpx - 0.3) / 0.7"),
'foot fade starts at 90% of the foot instead of 60%': ("(depthpx - 0.6) / 0.4", "(depthpx - 0.9) / 0.1"),
'TAPER_MAX 72 -> 40': ("TAPER_K, TAPER_MAX = 30.0, 72", "TAPER_K, TAPER_MAX = 30.0, 40"),
'RELIEF_EL 32 -> 3 (grazing relief light)': ("RELIEF_OFF_DEG, RELIEF_EL = 45.0, 32.0", "RELIEF_OFF_DEG, RELIEF_EL = 45.0, 3.0"),
'spec half vector ignores the light': ("Hv = L[None, :] + np.stack(", "Hv = 0 * L[None, :] + np.stack("),
'lateral tilt not blurred across columns': ("blurw(la, 2.0).astype", "la.astype"),
'contour interval percentile includes the tapering columns (weight condition removed)': ("& (wt[None, :] > 0.5)", ""),
'Fresnel sign flipped (rim strongest on faces standing at the viewer)': ("np.clip(1 - sl / nrm, 0, 1) ** 3", "np.clip(1 + sl / nrm, 0, 1) ** 3"),
'rays cache key ignores the refraction constant': ("R_EARTH, KREF, W, bool(detail)", "R_EARTH, W, bool(detail)"),
'rays cache key ignores the eye elevation': ("round(float(P['eye_m']), 3), ", ""),
'mosaic cache key ignores the tile list': ("sorted(tiles_for(p)), 3600", "3600"),
}
EQUIVALENT = {} # mutant name -> one line on why it is behaviourally equivalent (reported, not run)
def run(item):
name, (a, b) = item; d = tempfile.mkdtemp(); shutil.copytree(SRC, os.path.join(d, 'src'), ignore=shutil.ignore_patterns('tests', '__pycache__')); p = os.path.join(d, 'src', 'ridge.py'); s = open(p).read(); assert a in s, name
open(p, 'w').write(s.replace(a, b))
r = subprocess.run([sys.executable, os.path.join(HERE, 'test_ridge.py')], env=dict(os.environ, RIDGE_SRC=os.path.join(d, 'src')), capture_output=True, text=True); shutil.rmtree(d)
fails = [l.split('|')[0].strip() for l in r.stdout.splitlines() if l.startswith('FAIL')]
return ('KILLED ' if r.returncode else 'SURVIVED ') + name + ' -> ' + (', '.join(fails) or 'no fixture failed'), r.returncode == 0
if __name__ == '__main__':
from concurrent.futures import ThreadPoolExecutor
with ThreadPoolExecutor(int(os.environ.get('MUTATE_JOBS', '6'))) as ex: res = list(ex.map(run, [(n, v) for n, v in MUTANTS.items() if n not in EQUIVALENT]))
for line, _ in res: print(line, flush=True)
for n, why in EQUIVALENT.items(): print('EQUIVALENT ' + n + ' -> ' + why)
sys.exit(1 if any(sv for _, sv in res) else 0)

7
scripts/ridge-bake/tests/run.sh Executable file
View file

@ -0,0 +1,7 @@
#!/bin/sh
# Usage: RIDGE_PY=<venv python> RIDGE_TILES=<tiles dir> tests/run.sh [test_ridge.py|mutate.py]
# RIDGE_WORK (cache dir) defaults to a fresh directory under $TMPDIR.
: "${RIDGE_PY:?RIDGE_PY must name the python of the baker venv}"
export RIDGE_WORK="${RIDGE_WORK:-${TMPDIR:-/tmp}/ridge-bake-work}"; mkdir -p "$RIDGE_WORK"
HERE=$(cd "$(dirname "$0")" && pwd)
exec "$RIDGE_PY" -I "$HERE/${1:-test_ridge.py}"

382
scripts/ridge-bake/tests/test_ridge.py Normal file
View file

@ -0,0 +1,382 @@
"""Fixtures for the A+D sweep/chain/taper. Plain asserts, no pytest needed: venv/bin/python tests/test_ridge.py (exit code 1 on failure).
RIDGE_SRC selects the source directory (default ..) (the mutation runner points it at a patched copy).
Synthetic terrain has a known right answer; the Vienna test pins the skyline of the real bake (its rays come from stage0 for the vienna preset: RIDGE_TILES needs the tiles, or a cached rays file must exist)."""
import os, sys, math, traceback
SRC = os.environ.get('RIDGE_SRC') or os.path.join(os.path.dirname(os.path.abspath(__file__)), '..'); sys.path.insert(0, SRC) # RIDGE_WORK and RIDGE_TILES must be set (tests/run.sh sets them)
import numpy as np, ridge
from ridge import W, TAPER, TAPER_MAX, taper_len
EYE = 100.0; GROUND = 50.0 # ground below the eye: every ray climbs, so only the placed blocks make ridges
ridge.P = dict(id='synthetic', rmax_km=60)
D = ridge.schedule(False)
def az_of(x): return (x + 0.5) / W * 360.0
def col_of(az): return int(round(az / 360.0 * W - 0.5)) % W
def rays(hfun):
az = (np.arange(W) + 0.5) / W * 360.0; A = np.zeros((len(D), W), np.float32); Hh = np.zeros_like(A)
for k, dk in enumerate(D):
h = hfun(dk, az); Hh[k] = h; A[k] = np.degrees(np.arctan2(h - EYE - dk * dk / (2 * ridge.R_EARTH) * (1 - ridge.KREF), dk))
return A, Hh
def theta(h, d): return math.degrees(math.atan2(h - EYE - d * d / (2 * ridge.R_EARTH) * (1 - ridge.KREF), d))
def in_arc(az, a0, a1): return ((az - a0) % 360.0) <= ((a1 - a0) % 360.0)
def block(a0, a1, d0, d1, h): return lambda dk, az: np.where(in_arc(az, a0, a1) & (dk >= d0) & (dk < d1), h, 0.0)
def terrain(*blocks): return lambda dk, az: np.maximum.reduce([np.full_like(az, GROUND)] + [np.asarray(b(dk, az), float) for b in blocks])
def slot_at(r, col, ang):
"""the slot whose top at `col` is closest to `ang` (the layer that owns the block)"""
return min(r['slots'], key=lambda S: abs(S['top'][col] - ang))
def jump_bound(S, x, w=26):
"""steepest step allowed inside a taper: smoothstep slope 1.5 x the end drop / taper_len(drop) (+ margin), the end drop taken as max((highest top within w columns) - lo, 0.4)"""
seg = S['top'][[(x + k) % W for k in range(-w, w + 1)]]; drop = max(seg.max() - S['lo'][x], 0.4); return 1.6 * drop / float(taper_len(drop)) + 0.02 # the taper always falls at least 0.4 deg (top - 0.4 when the valley is shallower)
def active_pairs(S): return [x for x in range(W) if S['wt'][x] > 0 and S['wt'][(x + 1) % W] > 0]
def test_plateau_has_no_spike_and_no_cliff_beyond_taper():
"""Flat-topped plateau, az 100-140, 20-24 km, 900 m. Right answer: top = theta(900, 20 km) inside, nothing above it anywhere, outside the plateau the top only falls
(monotonically, within the taper), and no single-column step exceeds the taper bound."""
A, Hh = rays(terrain(block(100, 140, 20000, 24000, 900))); r = ridge.sweep(D, A, Hh); th = theta(900, 20000)
S = slot_at(r, col_of(120), th); top = S['top']
inside = np.arange(col_of(101), col_of(139)); assert np.abs(top[inside] - th).max() < 0.03, ('plateau top not flat/right', np.abs(top[inside] - th).max())
assert top.max() <= th + 0.03, ('spike above the plateau', top.max() - th)
for side in (-1, 1):
x0 = col_of(100) if side < 0 else col_of(140)
n = int(taper_len(max(th - S['lo'][col_of(120)], 0.4))); seg = np.array([top[(x0 + side * j) % W] for j in range(1, n)])
assert (np.diff(seg) <= 1e-9).all(), ('taper is not monotone away from the plateau', side, np.diff(seg).max())
assert seg[-1] < seg[0] - 0.3, 'taper does not die away'
assert S['wt'][(x0 + side * (n + 3)) % W] == 0, 'the layer is still drawn far outside the plateau'
j = np.abs(np.diff(top, append=top[0])); bad = [(x, round(float(j[x]), 3)) for x in active_pairs(S) if j[x] > jump_bound(S, x)]; assert not bad, ('column jump above the bound', bad[:5])
drop = th - S['lo'][col_of(120)]; width = int(((S['wt'] > 0) & (S['wt'] < 1))[col_of(140):col_of(140) + 90].sum()); assert drop > 1.5 and width >= 0.9 * min(72, 30 * drop), ('a tall end must fall over ~30 columns per degree of drop up to 72 (absolute, not the constant) (its drop is %.2f deg)' % drop, width)
def test_ridge_whose_two_ends_differ_does_not_borrow_the_far_end():
"""Plateau az 100-140 whose height climbs 600 -> 1200 m with azimuth. Right answer: left of the left end nothing is higher than the left end's own top, right of the right end
nothing is higher than the right end's top (the round-1 bug held the END value on both sides: a 0.9 deg step on the low side)."""
hfun = lambda dk, az: np.where(in_arc(az, 100, 140) & (dk >= 20000) & (dk < 24000), 600.0 + 600.0 * np.clip((az - 100.0) / 40.0, 0, 1), 0.0)
A, Hh = rays(terrain(hfun)); r = ridge.sweep(D, A, Hh); tl, tr = theta(600, 20000), theta(1200, 20000)
S = slot_at(r, col_of(120), (tl + tr) / 2); top = S['top']; xl, xr = col_of(100), col_of(140)
left = top[[(xl - j) % W for j in range(0, TAPER_MAX + 3)]]; right = top[[(xr + j) % W for j in range(0, TAPER_MAX + 3)]]
assert left.max() <= tl + 0.05, ('left of the low end rises above it', left.max() - tl)
assert right.max() <= tr + 0.05 and right[0] >= tr - 0.1, ('right end value wrong', right[:3])
assert (np.diff(left) <= 1e-9).all() and (np.diff(right) <= 1e-9).all(), 'tapers are not monotone'
def test_slot_union_only_lets_present_chains_compete():
"""Two chains in one slot: L1 present on cols 0-99 (top 0.2), L2 present on cols 200-299 (top 2.0) but its HELD value elsewhere is 0.5 > 0.2. Right answer: on cols 0-99
the slot is L1 (top 0.2, wt 1); a chain that is absent there must not win with its held value."""
def lay(c0, c1, top, held):
wt = np.zeros(W); wt[c0:c1] = 1.0; tp = np.full(W, held); tp[c0:c1] = top
return dict(top=tp, lo=np.full(W, -1.0), ph=np.full(W, 100.0), vh=np.full(W, 0.0), slope=np.full(W, 5.0), wt=wt, pd=np.full(W, 20000.0), vd=np.full(W, 19000.0), hd=np.full((ridge.NH, W), 19500.0), hh=np.full((ridge.NH, W), 50.0), dist=20000.0)
S = ridge.to_slots([lay(0, 100, 0.2, -1.0), lay(200, 300, 2.0, 0.5)])
assert len(S) == 1 and (S[0]['wt'][0:100] == 1).all() and np.allclose(S[0]['top'][0:100], 0.2), ('absent chain won the slot union', S[0]['top'][:3], S[0]['wt'][:3])
assert (S[0]['wt'][200:300] == 1).all() and np.allclose(S[0]['top'][200:300], 2.0)
def S_ok(S, i): return S['wt'][i] > 0 and S['wt'][(i + 1) % W] > 0
def test_step_cliff_is_one_step_and_levels_are_right():
"""Two plateaus side by side, 700 m (az 60-180) then 1300 m (180-300). Right answer: level theta_A left, theta_B right, one step at az 180 (col 3200)."""
A, Hh = rays(terrain(block(60, 180, 20000, 24000, 700), block(180, 300, 20000, 24000, 1300))); r = ridge.sweep(D, A, Hh)
tA, tB = theta(700, 20000), theta(1300, 20000); x = col_of(180)
SA = slot_at(r, col_of(120), tA); top = SA['top']
assert np.abs(top[col_of(62):x - TAPER_MAX - 2] - tA).max() < 0.03, 'level A wrong'
assert np.abs(top[x + 2:col_of(298)] - tB).max() < 0.03, 'level B wrong'
assert top.max() <= tB + 0.03, 'spike above the higher plateau'
j = np.abs(np.diff(top)); ends = (col_of(60), col_of(300)); big = [i for i in np.where(j > 0.2)[0] if S_ok(SA, i) and not (x - TAPER_MAX - 2 <= i <= x + 2) and min(abs(i - e) for e in ends) > TAPER_MAX + 2]
assert not big, ('a jump away from the step and the two outer ends', big[:6])
assert (top[x - TAPER_MAX - 2:x + 3].max() - tA) <= (tB - tA) + 0.03
def test_single_peak_has_one_maximum():
"""Gaussian peak 1800 m at az 90, 30 km, sigma 1.5 km. Right answer: one maximum at az 90 +-0.2, monotone either side inside the chain, no secondary bump."""
def peak(dk, az):
lat = dk * np.sin(np.radians(az - 90.0)); return 1800.0 * np.exp(-(lat ** 2 + (dk - 30000.0) ** 2) / (2 * 1500.0 ** 2))
A, Hh = rays(terrain(peak)); r = ridge.sweep(D, A, Hh); top_peak = theta(1800, 30000)
S = max(r['slots'], key=lambda S: S['top'][col_of(90)]); top = S['top']
xm = int(np.argmax(top[col_of(60):col_of(120)])) + col_of(60); assert abs(az_of(xm) - 90.0) < 0.25, ('maximum at', az_of(xm))
assert abs(top[xm] - top_peak) < 0.1, (top[xm], top_peak)
seg = top[col_of(60):col_of(120)]; peaks = [i for i in range(1, len(seg) - 1) if seg[i] >= seg[i - 1] and seg[i] > seg[i + 1] and (seg[i] - min(seg[max(0, i - 30):i + 30])) > 0.01]
assert len(peaks) == 1, ('secondary local maxima', len(peaks))
def test_wraparound_ridge_across_north_is_seamless_and_rotation_invariant():
"""Plateau az 350 -> 10 (across 0/360). Right answer: flat top, no step across the seam, one chain crossing it, and the result does not depend on where column 0 is."""
A, Hh = rays(terrain(block(350, 10, 20000, 24000, 900))); r = ridge.sweep(D, A, Hh); th = theta(900, 20000)
S = slot_at(r, 0, th); top = S['top']
assert abs(top[W - 1] - top[0]) < 0.003, ('seam step', top[W - 1] - top[0])
assert np.abs(top[[col_of(351), col_of(355), W - 1, 0, col_of(5), col_of(9)]] - th).max() < 0.03
near20 = [l for l in r['keep'] if abs(math.log(l['dist'] / 20000.0)) < 0.3]
assert len(near20) == 1 and near20[0]['x0'] + near20[0]['span'] > W, ('the 20 km plateau must be ONE chain crossing the seam', [(l['x0'], l['span']) for l in near20])
R = 1600; r2 = ridge.sweep(D, np.roll(A, -R, axis=1), np.roll(Hh, -R, axis=1))
for S1 in r['slots']:
S2 = next(s for s in r2['slots'] if s['slot'] == S1['slot']); d = np.abs(np.roll(S2['top'], R) - S1['top']).max(); assert d < 1e-9, ('not rotation invariant', S1['slot'], d)
def test_vienna_skyline_has_no_invented_jump_and_seam_is_seamless():
"""Real bake: no per-column step above the stated bound. In a chain (wt == 1) a step above 0.12 deg is allowed only if both columns carry the raw span peak that wins there (a real step of the DEM, e.g. a roof edge); in a taper the
steepest allowed step is 1.6 x (top-lo)/TAPER + 0.02. The envelope across the seam must not jump."""
saved = ridge.P
ridge.configure(ridge.load_preset('vienna')); d, A, Hh, _ = ridge.stage0(False) # cached under the key of its own inputs (computed once if absent)
ridge.P = dict(id='vienna', rmax_km=95); r = ridge.sweep(d, A, Hh); ridge.P = saved
worst = []; cols = r['cols']
real = lambda x, top: any(abs(s[0] - top) < 1e-4 for s in cols[x % W]) # the top at x is an actual raw span peak: a step between two such columns is a step of the DEM
for S in r['slots']:
j = np.abs(np.diff(S['top'], append=S['top'][0]))
for x in active_pairs(S):
bound = jump_bound(S, x) if min(S['wt'][x], S['wt'][(x + 1) % W]) < 1 else 0.12
if j[x] > bound and not (S['wt'][x] == 1 and S['wt'][(x + 1) % W] == 1 and real(x, S['top'][x]) and real(x + 1, S['top'][(x + 1) % W])):
worst.append((S['slot'], x, round(float(j[x]), 3), round(float(bound), 3), round(float(S['wt'][x]), 2)))
assert not worst, ('Vienna jumps above the bound', worst[:6])
sky = np.max([S['top'] for S in r['slots']], axis=0); assert abs(sky[W - 1] - sky[0]) < 0.05, ('envelope jumps across the seam', sky[W - 1] - sky[0])
for S in r['slots']:
if S['wt'][W - 1] > 0 and S['wt'][0] > 0: assert abs(S['top'][W - 1] - S['top'][0]) < 0.05, ('layer jumps across the seam', S['slot'])
def _shadow_fixture():
"""flat 100 m terrain with one 1000 m wall at lat 48.45-48.46 (5.6-6.7 km north of the sample point), sample point lat 48.40 lon 16.5 at 103 m"""
m = np.full((3600, 3600), 100.0, np.float32); y0 = int((49 - 48.46) * 3600); y1 = int((49 - 48.45) * 3600); m[y0:y1, :] = 1000.0
ridge.P = dict(id='syn', lat=48.40, lon=16.5, rmax_km=60); ridge._mos[('syn', ridge.mosaic_key(ridge.P))] = (m, 49, 16)
return [dict(spts=dict(lat=np.full((3, W), 48.40), lon=np.full((3, W), 16.5), h=np.full((3, W), 103.0)))]
def test_cast_shadow_follows_the_sun_ray_over_the_dem():
"""Right answer: sun due north at 5 deg: the ray is 580 m at the wall, the wall is 1000 m => fully shadowed (1.0). Sun due north at 30 deg clears the wall (0.0). Sun due south: the wall is behind the sun (0.0).
Sun at 0.2 deg: no direct light, no shadow (0.0)."""
L = _shadow_fixture(); f = lambda az, el: float(ridge.cast_shadows(L, az, el)[0].mean())
assert f(0.0, 5.0) > 0.99, ('wall does not shadow a low sun behind it', f(0.0, 5.0))
assert f(0.0, 30.0) < 0.01, ('high sun clears the wall', f(0.0, 30.0))
assert f(180.0, 5.0) < 0.01, ('wall casts a shadow toward the sun', f(180.0, 5.0))
assert f(0.0, 0.2) == 0.0
def test_cast_shadow_east_west_and_longitude_scale():
"""Wall of 1000 m at lon 16.575-16.585 (5.5-6.2 km EAST of the sample point at lat 48.4, cos(lat) matters), flat 100 m elsewhere. Sun due east at 5 deg: shadowed (580 m ray vs 1000 m).
Sun due west: clear. Without the cos(lat) factor the wall would sit at 8.3 km and a 6 deg ray (980 m) would barely touch it, so the fully-shadowed assertion at 6 deg pins the scale."""
m = np.full((3600, 3600), 100.0, np.float32); x0 = int((16.575 - 16) * 3600); x1 = int((16.585 - 16) * 3600); m[:, x0:x1] = 1000.0
ridge.P = dict(id='syn2', lat=48.40, lon=16.5, rmax_km=60); ridge._mos[('syn2', ridge.mosaic_key(ridge.P))] = (m, 49, 16)
L = [dict(spts=dict(lat=np.full((3, W), 48.40), lon=np.full((3, W), 16.5), h=np.full((3, W), 103.0)))]; f = lambda az, el: float(ridge.cast_shadows(L, az, el)[0].mean())
assert f(90.0, 5.0) > 0.99, ('east wall does not shadow a low east sun', f(90.0, 5.0))
assert f(270.0, 5.0) < 0.01, ('east wall shadows a west sun', f(270.0, 5.0))
assert f(90.0, 6.0) > 0.99 and f(90.0, 12.0) < 0.01, ('east-wall shadow edge is at the wrong elevation (wall is 5.5 km east of the point: the sun ray at 6 deg is 660 m there, at 12 deg 1.2 km)', f(90.0, 6.0), f(90.0, 12.0))
def test_shadow_darkens_the_face_it_covers():
"""Whole pipeline on a synthetic ramp (slope 0.35, lit by the relief light): with the shadow march forced to 1 every face pixel is darker than with it forced to 0 (the shadow removes direct sun); forced 0 equals shadows off."""
A, Hh = rays(terrain(lambda dk, az: np.where(in_arc(az, 100, 140) & (dk >= 20000) & (dk < 24000), 50.0 + 0.35 * (dk - 20000), 0.0))); ridge.P = dict(id='synthetic', rmax_km=60, lat=48.0, lon=16.0, floor_agl_m=30.0, s_ref_deg=8.0, win0=(col_of(120) - 640) % W)
r = ridge.sweep(D, A, Hh); ridge.set_view(r['slots']); L = ridge.bake(r['slots']); pf = ridge.statics(); real = ridge.cast_shadows
try:
ridge.cast_shadows = lambda layers, az, el: [np.ones((3, W)) for _ in layers]; dark = ridge.relight(L, 300.0, 30.0, *pf)
ridge.cast_shadows = lambda layers, az, el: [np.zeros((3, W)) for _ in layers]; lit = ridge.relight(L, 300.0, 30.0, *pf)
finally: ridge.cast_shadows = real
off = ridge.relight(L, 300.0, 30.0, *pf, shadows=False); c = slice(col_of(105), col_of(135)); rows = slice(0, ridge.H); Ly = max(L, key=lambda l: ((l['m'][rows, c] > 0.99) & (l['w'][rows, c] < 0.05)).sum()); m = (Ly['m'][rows, c] > 0.99) & (Ly['w'][rows, c] < 0.05) # body pixels (not the mist foot) of the layer that owns the plateau
assert m.any(), 'no solid face pixels to compare'
assert dark[rows, c][m].mean() < lit[rows, c][m].mean() - 0.02, ('shadow does not darken the face', dark[rows, c][m].mean(), lit[rows, c][m].mean())
assert np.abs(lit - off).max() < 1e-4, 'a zero shadow differs from shadows off'
# ---------------------------------------------------------------------------------------- shading fixtures (round 3): the face history -> per-pixel normal
def _bake_view(hfun, az0):
A, Hh = rays(terrain(hfun)); ridge.P = dict(id='synthetic', rmax_km=60, lat=48.0, lon=16.0, floor_agl_m=30.0, s_ref_deg=8.0, win0=(col_of(az0) - 640) % W)
r = ridge.sweep(D, A, Hh); ridge.set_view(r['slots']); return r, ridge.bake(r['slots'])
def _owner(r, L, col, ang):
S = slot_at(r, col, ang); return S, next(Ly for Ly in L if Ly['dist'] == S['dist'])
CONE = dict(Dc=20000.0, Rb=6000.0, Hc=1800.0, az0=90.0)
def _cone(dk, az):
E = dk * np.sin(np.radians(az)) - CONE['Dc'] * math.sin(math.radians(CONE['az0'])); N = dk * np.cos(np.radians(az)) - CONE['Dc'] * math.cos(math.radians(CONE['az0']))
return np.where(np.hypot(E, N) < CONE['Rb'], CONE['Hc'] * (1 - np.hypot(E, N) / CONE['Rb']), 0.0)
def _cone_pixels(sun):
"""(analytic n.L, baked n.L) at pixels of the cone layer: the analytic one marches the ray finely to the first sample whose running-max angle reaches the pixel's elevation and takes the cone's own normal there."""
r, L = _bake_view(_cone, CONE['az0']); S, Ly = _owner(r, L, col_of(CONE['az0']), theta(CONE['Hc'], CONE['Dc'])); nd, _ = ridge.ndl_field(Ly, sun); dd = np.arange(500.0, 60000.0, 4.0); m = CONE['Hc'] / CONE['Rb']; out = []
for x in range(300, 1000, 7):
c = (ridge.P['win0'] + x) % W
if S['wt'][c] < 1: continue
az = (c + 0.5) / W * 360.0; h = _cone(dd, np.full_like(dd, az)); ang = np.degrees(np.arctan2(h - EYE - dd * dd / (2 * ridge.R_EARTH) * (1 - ridge.KREF), dd)); rm = np.maximum.accumulate(ang)
for f in (0.6, 0.7, 0.8, 0.9, 0.95):
row = int(round((ridge.TOP_DEG - (Ly['lo'][c] + f * (Ly['ang'][c] - Ly['lo'][c]))) * ridge.PXDEG_V - 0.5)); d = dd[min(np.searchsorted(rm, ridge.rowang[row]), len(dd) - 1)]
E = d * math.sin(math.radians(az)) - CONE['Dc'] * math.sin(math.radians(CONE['az0'])); N = d * math.cos(math.radians(az)) - CONE['Dc'] * math.cos(math.radians(CONE['az0'])); rr = math.hypot(E, N)
if rr >= 0.85 * CONE['Rb']: continue # on the cone proper, away from the foot kink (the +-2 node slope window rounds the kink)
n = np.array([m * E / rr, m * N / rr, 1.0]); n /= np.linalg.norm(n); out.append((float(n @ sun), float(nd[row, c])))
return np.array(out), Ly, nd
def test_cone_is_shaded_per_pixel_with_the_analytic_normal():
"""Cone (apex 1800 m, base radius 6 km, 20 km away), sun from the west at 20 deg. Right answer: the shading of each pixel is the cone's own n.L at the point that pixel sees, so it changes down every
column (the radial direction turns as the ray climbs) and across the face; the round-2 baker stored ONE slope per peak and gave every pixel of a column the same value (within-column std 0)."""
sun = ridge.sun_vec(250.0, 20.0); v, Ly, nd = _cone_pixels(sun); assert len(v) > 100, len(v)
err = np.abs(v[:, 0] - v[:, 1]); assert np.median(err) < 0.05 and err.mean() < 0.07, ('baked n.L is not the analytic one', np.median(err), err.mean())
assert np.ptp(v[:, 1]) > 0.8 * np.ptp(v[:, 0]), ('baked shading range is compressed', np.ptp(v[:, 1]), np.ptp(v[:, 0]))
cs = [nd[int((ridge.TOP_DEG - Ly['ang'][c]) * ridge.PXDEG_V) + 6:int((ridge.TOP_DEG - Ly['lo'][c]) * ridge.PXDEG_V) - 6, c].std() for c in (ridge.P['win0'] + np.arange(400, 900, 20)) % W]
assert np.mean(cs) > 0.03, ('shading is constant down a column', np.mean(cs))
def test_plane_face_has_the_planes_slope_and_no_lateral_tilt():
"""Planar ramp facing the viewer, az 60-120, rising 0.35 per metre of distance from 20 km to 24 km, then a drop. Right answer: every face pixel's front slope is 0.35 (19.3 deg) and its lateral tilt is 0,
whatever its row (one slope per peak could not tell a plane from a curved face, only the history can)."""
ramp = lambda dk, az: np.where(in_arc(az, 60, 120) & (dk >= 20000) & (dk < 24000), 50.0 + 0.35 * (dk - 20000), 0.0)
r, L = _bake_view(ramp, 90.0); S, Ly = _owner(r, L, col_of(90), theta(50.0 + 0.35 * 4000, 24000)); cols = (ridge.P['win0'] + np.arange(560, 720)) % W
rows = np.arange(int((ridge.TOP_DEG - Ly['ang'][col_of(90)]) * ridge.PXDEG_V) + 8, int((ridge.TOP_DEG - 1.0) * ridge.PXDEG_V)); sl = Ly['sl'][np.ix_(rows, cols)]; la = Ly['la'][np.ix_(rows, cols)]
assert len(rows) > 20 and abs(np.median(sl) - 0.35) < 0.03 and np.percentile(np.abs(sl - 0.35), 90) < 0.08, ('front slope is not the plane\'s', np.median(sl), np.percentile(np.abs(sl - 0.35), 90))
assert np.abs(la).max() < 0.08, ('lateral tilt on a plane', np.abs(la).max())
def test_ridge_step_tilts_the_face_only_at_the_step():
"""Two plateaus (700 m / 1300 m) meeting at az 180. Right answer: the face is level laterally (la ~ 0) away from the step, and strongly tilted at it (the height changes by 600 m between two neighbouring columns)."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 60, 180) & (dk >= 20000) & (dk < 24000), 700.0, np.where(in_arc(az, 180, 300) & (dk >= 20000) & (dk < 24000), 1300.0, 0.0)), 180.0)
S, Ly = _owner(r, L, col_of(120), theta(700, 20000)); cols = lambda a, b: (ridge.P['win0'] + np.arange(a, b)) % W; rows = slice(int((ridge.TOP_DEG - theta(1300, 20000)) * ridge.PXDEG_V), int((ridge.TOP_DEG - 0.3) * ridge.PXDEG_V))
far = np.abs(Ly['la'][rows, :][:, cols(200, 500)]); assert far.max() < 0.05, ('lateral tilt away from the step', far.max())
near = Ly['la'][rows, :][:, cols(620, 660)]; assert np.abs(near).max() > 0.2, ('no tilt at the step', np.abs(near).max())
def test_roof_height_bump_in_the_foreground_is_not_a_ridge():
"""Flat ground with a 19 m x 40 m block 700 m away (a roof / tree in a DSM). Right answer: no layer: a peak needs PROM_M (25 m) of height over its valley, and the 4 km NEAR_CAP shoulder rule obeys the same
rule (round 2 committed the shoulder with no prominence test: an arch 66 columns wide stood in the foreground). A 60 m block IS a ridge."""
bump = lambda h: (lambda dk, az: np.where(in_arc(az, 88, 92) & (dk >= 700) & (dk < 740), GROUND + h, 0.0))
A, Hh = rays(terrain(bump(19.0))); r = ridge.sweep(D, A, Hh); assert not [l for l in r['keep'] if l['dist'] < 4000], ('a 19 m bump became a ridge', [(l['x0'], l['span'], round(l['dist'])) for l in r['keep']])
A, Hh = rays(terrain(bump(60.0))); r = ridge.sweep(D, A, Hh); assert [l for l in r['keep'] if l['dist'] < 1000], 'a 60 m block must be a ridge'
def test_contour_lines_are_never_denser_than_9_px_over_most_of_the_face():
"""Right answer (cone): the iso-height lines drawn on a layer are at least CONTOUR_MIN_PX apart over three quarters of its face (round 2 picked the interval from the SHALLOWEST columns, so the steep centre
of the mountain carried a line every 2-3 px: the fine horizontal banding)."""
r, L = _bake_view(_cone, CONE['az0']); S, Ly = _owner(r, L, col_of(CONE['az0']), theta(CONE['Hc'], CONE['Dc'])); sp = []
for c in (ridge.P['win0'] + np.arange(450, 850, 10)) % W:
r0, r1 = int((ridge.TOP_DEG - Ly['ang'][c]) * ridge.PXDEG_V) + 4, int((ridge.TOP_DEG - Ly['lo'][c]) * ridge.PXDEG_V) - 4; line = np.where(Ly['dpx'][r0:r1, c] < 0.5)[0]
if len(line) > 1: g = np.diff(line); sp += list(g[g > 1])
assert len(sp) > 20 and np.percentile(sp, 25) >= 8, ('contour lines too dense', np.percentile(sp, 25), len(sp))
def test_no_layer_paints_a_wall_down_the_frame():
"""Right answer: a layer's alpha is 0 further below its own crest than its foot allows (<= 210 px): its column boundary leaves no vertical edge through the floor (round 2 painted the body to the bottom of the
frame wherever the chain existed: hard walls below every chain end)."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 100, 140) & (dk >= 20000) & (dk < 24000), 900.0, 0.0), 120.0)
for Ly in L: below = np.arange(ridge.H)[:, None] - (ridge.TOP_DEG - Ly['ang'])[None, :] * ridge.PXDEG_V; assert Ly['a'][below > 215].max() == 0, ('layer painted below its foot', Ly['rank'])
S = max(r['slots'], key=lambda S: S['top'][col_of(120)])
for Ly in L:
if Ly['dist'] != S['dist']: continue
end = [c for c in range(W) if 0 < S['wt'][c] < 0.12]; assert end, 'fixture has no taper end'
bad = [(c, int((Ly['a'][:, c] > 0.1).sum())) for c in end if (Ly['a'][:, c] > 0.1).sum() > 0.85 * max(Ly['ang'][c] - Ly['lo'][c], 0) * ridge.PXDEG_V + 4]; assert not bad, ('a chain end hangs a curtain thicker than the ridge is tall', bad[:4])
def test_slot_is_rotation_invariant_in_columns_where_no_chain_is_present():
"""Two plateaus in one distance slot (az 100-110 at 700 m, az 200-215 at 1100 m, both 20-24 km). Right answer: rotating the ring by 1600 columns changes nothing, in EVERY column of the slot including those where no chain is present
(round 3 first picked the first chain in list order there; list order depends on where column 0 is, and those values feed the blurs at every chain end: Half Dome's rotation test went from 0 to 97 px)."""
A, Hh = rays(terrain(block(100, 110, 20000, 24000, 700), block(200, 215, 20000, 24000, 1100))); r = ridge.sweep(D, A, Hh); R = 2800; r2 = ridge.sweep(D, np.roll(A, -R, axis=1), np.roll(Hh, -R, axis=1)) # 2800 columns carries the az 200-215 block through column 0: the chains swap places in list order
for S1 in r['slots']:
S2 = next(s for s in r2['slots'] if s['slot'] == S1['slot'])
for k in ('top', 'lo', 'wt', 'pd', 'ph', 'vh', 'slope'): assert np.abs(np.roll(S2[k], R) - S1[k]).max() < 1e-9, ('slot not rotation invariant', S1['slot'], k, np.abs(np.roll(S2[k], R) - S1[k]).max())
for k in ('hd', 'hh'): assert np.abs(np.roll(S2[k], R, axis=1) - S1[k]).max() < 1e-6, ('slot history not rotation invariant', S1['slot'], k)
def test_prominent_far_ridge_behind_a_basin_is_found():
"""Terrace (50 m, eye 100 m) ending at 10 km, a 0 m basin, a 60 m ridge at 30 km (az 80-100). Right answer: the ridge is a layer at 30 km: it stands 0.5 deg above its visible valley, which is 4 times the 25 m / 30 km bar of 0.048 deg.
(Round 3's first cut compared raw heights, ridge 60 m vs the valley's 0 m... and lost ridges whose valley sample is near and high.)"""
hf = lambda dk, az: np.where(dk < 10000, 50.0, np.where(in_arc(az, 80, 100) & (dk >= 30000) & (dk < 30500), 60.0, np.where(in_arc(az, 80, 100) & (dk >= 50000) & (dk < 55000), 300.0, 0.0))) # the 300 m wall at 50 km keeps the ridge from being committed as the ray's last peak
A, Hh = rays(hf); r = ridge.sweep(D, A, Hh); far = [l for l in r['keep'] if abs(l['dist'] - 30000) < 1500]; assert far and far[0]['span'] > 200, ('the 30 km ridge is missing', [(l['x0'], l['span'], round(l['dist'])) for l in r['keep']])
def test_lone_foreground_stub_is_kept_but_a_stub_inside_a_near_ridge_is_not():
"""(a) a 0.5 deg wide, 60 m block at 1.7 km with nothing else near: the nearest layer must cover it (a hole would show the far ridge's body). (b) the same stub in front of a long near ridge (az 100-160, 150 m, 1.5 km):
its columns are covered by that ridge, so it is dropped (round 2 kept every foreground stub: lone arches)."""
A, Hh = rays(terrain(block(90, 90.5, 1700, 1740, 110))); r = ridge.sweep(D, A, Hh); assert [l for l in r['keep'] if l['dist'] < 2500], 'lone foreground stub dropped: hole in the nearest layer'
A, Hh = rays(terrain(block(100, 160, 1500, 1600, 150), block(129, 129.5, 600, 640, 160))); r = ridge.sweep(D, A, Hh); assert not [l for l in r['keep'] if l['dist'] < 1000], ('a stub inside a near ridge became an arch', [(l['x0'], l['span'], round(l['dist'])) for l in r['keep']])
def test_tilted_plane_has_the_planes_front_slope_and_lateral_tilt():
"""Plane h = 50 + 0.3 (E - 20000) - 0.2 N over az 85-95, 20-24 km, then a drop. At az 90 the radial is +E and the tangent is -N, so n ~ (-0.3 r - la t, 1) with la = -(dh/dN) = +0.2. Right answer: front slope 0.3, lateral tilt 0.2, each within 0.04
(a lateral scale of 2x, or no tilt, fails)."""
def plane(dk, az):
E = dk * np.sin(np.radians(az)); N = dk * np.cos(np.radians(az)); return np.where(in_arc(az, 85, 95) & (dk >= 20000) & (dk < 24000), 50.0 + 0.3 * (E - 20000.0) - 0.2 * N, 0.0)
r, L = _bake_view(plane, 90.0); S, Ly = _owner(r, L, col_of(90), theta(50.0 + 0.3 * 4000, 24000)); cols = (ridge.P['win0'] + np.arange(600, 680)) % W
rows = np.arange(int((ridge.TOP_DEG - Ly['ang'][col_of(90)]) * ridge.PXDEG_V) + 8, int((ridge.TOP_DEG - 1.0) * ridge.PXDEG_V)); sl = Ly['sl'][np.ix_(rows, cols)]; la = Ly['la'][np.ix_(rows, cols)]
assert len(rows) > 15 and abs(np.median(sl) - 0.3) < 0.04 and abs(np.median(la) - 0.2) < 0.04, ('plane slope / tilt wrong', np.median(sl), np.median(la))
def test_history_keeps_only_what_the_eye_sees():
"""Samples (angle, distance, height): (0, 100, 0) (2, 200, 10) (1, 300, 5) (1.5, 350, 8) (3, 400, 30). The third and fourth are hidden behind the second (the fourth RISES from the third, so only a running maximum drops it). Right answer: at elevation 2.5 the face is half way between the second and the fourth sample: distance 300, height 20
(without the running maximum the non-monotone angles make the interpolation garbage)."""
d, h = ridge.strip_history([0.0, 2.0, 1.0, 1.5, 3.0], [100.0, 200.0, 300.0, 350.0, 400.0], [0.0, 10.0, 5.0, 8.0, 30.0], 0, 4); t = np.linspace(0, 3.0, ridge.NH); k = int(np.argmin(np.abs(t - 2.5)))
assert abs(d[k] - (200 + 200 * (t[k] - 2) / 1.0)) < 1e-6 and abs(h[k] - (10 + 20 * (t[k] - 2) / 1.0)) < 1e-6 and d[0] == 100 and d[-1] == 400, (d[k], h[k])
def test_no_lateral_tilt_across_an_occlusion_edge():
"""Block A (az 60-130, 15 km out, 700 m) and block B (az 110-180, 21 km out, 940 m): the same top angle, the same distance slot, so the slot owner switches where their tapers cross. Right answer: at every column where the owner's
distance jumps by more than 2 km the two sides are not one surface, so the lateral tilt there is 0 (it saturated at the 1.5 clamp over ~8 columns on Matterhorn)."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 60, 130) & (dk >= 15000) & (dk < 19000), 700.0, np.where(in_arc(az, 110, 180) & (dk >= 21000) & (dk < 25000), 940.0, 0.0)), 120.0)
S, Ly = _owner(r, L, col_of(80), theta(700, 15000)); assert S['dist'] == ridge.KNOTS[3]
sw = [x for x in range(W - 1) if S['wt'][x] > 0 and S['wt'][x + 1] > 0 and abs(S['pd'][x + 1] - S['pd'][x]) > 2000]; assert sw, 'fixture has no owner switch'
for x in sw:
rows = np.arange(max(int((ridge.TOP_DEG - Ly['ang'][x]) * ridge.PXDEG_V) + 4, 0), min(int((ridge.TOP_DEG - Ly['lo'][x]) * ridge.PXDEG_V) - 4, ridge.H)); la = Ly['la'][np.ix_(rows, [x - 1, x, x + 1, x + 2])]
assert np.abs(la).max() < 0.5, ('tilt across an occlusion edge at column', x, np.abs(la).max())
# ---------------------------------------------------------------------------------------- relief light (round 3): n.L uses an authored light, colour / haze / shadows keep the real sun
def _gully(dk, az):
"""A ramp facing the viewer (0.35 slope, az 50-130, 20-24 km) cut by gullies that run down the slope: height +300 m sin(2 pi N / 3 km), so the lateral tilt swings between -0.63 and +0.63 across the face."""
N = dk * np.cos(np.radians(az)); return np.where(in_arc(az, 50, 130) & (dk >= 20000) & (dk < 24000), 50.0 + 0.35 * (dk - 20000) + 300.0 * np.sin(2 * np.pi * N / 3000.0), 0.0)
def _relief_probe(az_s, el_s):
"""Relight the gullied ramp seen from az 90 with spies on ndl_field (the L it shades with and the n.L it returns for the ramp's layer) and on cast_shadows (the sun it marches toward). Returns (image, Ls, n.L of the ramp layer, shadow suns, layer)."""
r, L = _bake_view(_gully, 90.0); pf = ridge.statics(); S = max(r['slots'], key=lambda S: S['top'][col_of(90)]); Ly = next(l for l in L if l['dist'] == S['dist']); real_nd, real_cs = ridge.ndl_field, ridge.cast_shadows; Ls, nds, suns = [], [], []
def nd_spy(layer, Lv): out = real_nd(layer, Lv); Ls.append(np.array(Lv)); nds.append((layer is Ly, out[0])); return out
def cs_spy(layers, az, el): suns.append((az, el)); return [np.zeros((3, W)) for _ in layers]
ridge.ndl_field, ridge.cast_shadows = nd_spy, cs_spy
try: img = ridge.relight(L, az_s, el_s, *pf)
finally: ridge.ndl_field, ridge.cast_shadows = real_nd, real_cs
return img, Ls, next(n for own, n in nds if own), suns, Ly
def _face_std(nd, Ly):
"""n.L std over the ramp's face pixels (6 px under the crest down to 1 deg elevation), columns 500-780 of the window"""
return float(np.concatenate([nd[int((ridge.TOP_DEG - Ly['ang'][c]) * ridge.PXDEG_V) + 6:int((ridge.TOP_DEG - 1.0) * ridge.PXDEG_V), c] for c in (ridge.P['win0'] + np.arange(500, 780, 5)) % W]).std())
def test_sun_behind_the_viewer_still_models_the_face():
"""Gullied ramp seen from az 90, sun directly behind the viewer (az 270, 20 deg up). Lit head-on, a gully's two flanks shade alike (n.L std over the face 0.034 with the real sun as the shading light);
the relief light 45 deg off the axis lights one flank and leaves the other dim (0.23). Right answer: std over the face pixels > 0.12."""
img, Ls, nd, suns, Ly = _relief_probe(270.0, 20.0); b = _face_std(nd, Ly); assert b > 0.12, ('face is flat with the sun behind the viewer', b)
def test_relief_light_is_45_deg_off_the_view_axis_on_the_suns_side():
"""View centre az 90 (the viewer looks east). The relief light is the view centre +/- 45 deg on the side of the view axis the real sun is on. Sun at az 300 (north of the axis, the viewer's left) -> relief az 90-45 = 45;
sun at az 240 (south, the viewer's right) -> relief az 90+45 = 135. Elevation is 25..40 deg (an absolute range, not the constant) both times, whatever the sun's elevation."""
for sun_az, want in ((300.0, 45.0), (240.0, 135.0)):
_, Ls, _, _, _ = _relief_probe(sun_az, 25.0); assert Ls, 'no shading call'; L = Ls[0]; az = math.degrees(math.atan2(L[0], L[1])) % 360.0; el = math.degrees(math.asin(L[2]))
assert abs(az - want) < 0.5 and 25.0 <= el <= 40.0, ('relief direction wrong', sun_az, az, el)
def test_sun_elevation_changes_colour_but_not_the_relief_direction():
"""Same sun azimuth, elevation 10 deg vs 40 deg: the relief direction is identical, the picture is not (sun colour / haze / strength change), and at -10 deg (night) the face is darker than at 40."""
i1, L1, nd1, _, Ly = _relief_probe(300.0, 10.0); i2, L2, nd2, _, _ = _relief_probe(300.0, 40.0); i3, _, _, _, _ = _relief_probe(300.0, -10.0)
assert np.abs(L1[0] - L2[0]).max() < 1e-12 and np.abs(nd1 - nd2).max() < 1e-12, 'the relief direction follows the sun elevation'
c = (ridge.P['win0'] + 640) % W; rows = slice(int((ridge.TOP_DEG - Ly['ang'][c]) * ridge.PXDEG_V) + 8, int((ridge.TOP_DEG - Ly['ang'][c]) * ridge.PXDEG_V) + 40)
assert np.abs(i1[rows, 640] - i2[rows, 640]).mean() > 0.01, 'colour does not follow the real sun'; assert i3[rows, 640].mean() < 0.6 * i2[rows, 640].mean(), ('night is not dark', i3[rows, 640].mean(), i2[rows, 640].mean())
def test_cast_shadows_march_toward_the_real_sun_not_the_relief_light():
"""relight must hand cast_shadows the real sun's azimuth and elevation, not the relief light's."""
_, _, _, suns, _ = _relief_probe(300.0, 17.0); assert suns == [(300.0, 17.0)], suns
# ---------------------------------------------------------------------------------------- Soul round-3 survivors
def test_floor_band_has_no_column_wise_luminance_step():
"""Plateau az 100-140 (20-24 km) seen from az 120. Right answer: in the floor band (the frame's bottom 80 rows) no column differs from its neighbour by more than the grain (baseline 0.009, a painted wall 0.25): the compositor must use the layer's
alpha (body dissolved at the foot), not its mask (body to the bottom of the frame, a wall at every chain end)."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 100, 140) & (dk >= 20000) & (dk < 24000), 900.0, 0.0), 120.0); img = ridge.relight(L, 200.0, 20.0, *ridge.statics(), shadows=False)
lum = img @ np.array([0.2126, 0.7152, 0.0722]); st = np.abs(np.diff(lum[ridge.H - 80:], axis=1)).mean(0); assert st.max() < 0.03, ('wall: column-wise luminance step in the floor band', st.max(), int(np.argmax(st)))
def test_foot_fade_starts_at_60_percent_of_the_foot_and_the_body_is_opaque_above_it():
"""A layer's alpha is 1 down to 60% of its foot length, then smoothsteps to 0 at the foot. Right answer on a plateau column: rows with alpha > 0.99 / rows with alpha > 0.001 = 0.6 (+-0.04), and it has the full 210 px foot at most
(a fade that starts early makes the face translucent, one that starts late leaves a wall)."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 100, 140) & (dk >= 20000) & (dk < 24000), 900.0, 0.0), 120.0); S, Ly = _owner(r, L, col_of(120), theta(900, 20000)); a = Ly['a'][:, col_of(120)]
full, nz = int((a > 0.99).sum()), int((a > 0.001).sum()); assert nz > 100, nz; assert abs(full / nz - 0.6) < 0.04, ('fade does not start at 60% of the foot', full, nz)
def test_lateral_tilt_is_smoothed_across_columns():
"""Two plateaus (700 / 1300 m) meeting at az 180: the raw lateral tilt is a spike at the step. Right answer: the baked tilt is blurred over a few columns: no column-to-column change above 0.5 (unblurred: 0.9), and the tilt is still > 0.2 at the step."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 60, 180) & (dk >= 20000) & (dk < 24000), 700.0, np.where(in_arc(az, 180, 300) & (dk >= 20000) & (dk < 24000), 1300.0, 0.0)), 180.0)
S, Ly = _owner(r, L, col_of(120), theta(700, 20000)); rows = slice(int((ridge.TOP_DEG - theta(1300, 20000)) * ridge.PXDEG_V), int((ridge.TOP_DEG - 0.3) * ridge.PXDEG_V)); la = np.abs(Ly['la'][rows, :][:, (ridge.P['win0'] + np.arange(600, 680)) % W]).max(0)
assert la.max() > 0.2 and np.abs(np.diff(la)).max() < 0.5, ('tilt not smoothed', la.max(), np.abs(np.diff(la)).max())
def test_contour_interval_is_chosen_from_full_weight_columns_only():
"""Ramp (slope 1.0, az 110-130, 20-22 km) whose tapering ends carry coarser metres-per-pixel. Right answer: the contour spacing in the core is 9-10 px (interval from the chain's own face), not 13-14 (the taper columns pulled the percentile up). median <= 12 and p25 >= 8."""
r, L = _bake_view(lambda dk, az: np.where(in_arc(az, 110, 130) & (dk >= 20000) & (dk < 22000), 50.0 + 1.0 * (dk - 20000), 0.0), 120.0); S, Ly = _owner(r, L, col_of(120), theta(2050, 22000)); sp = []
for c in range(col_of(120) - 60, col_of(120) + 61, 6):
r0, r1 = int((ridge.TOP_DEG - Ly['ang'][c]) * ridge.PXDEG_V) + 4, int((ridge.TOP_DEG - Ly['lo'][c]) * ridge.PXDEG_V) - 4; line = np.where(Ly['dpx'][r0:r1, c] < 0.5)[0]
if len(line) > 1: g = np.diff(line); sp += list(g[g > 1])
assert len(sp) > 100 and np.percentile(sp, 25) >= 8 and np.median(sp) <= 12, ('contour spacing', len(sp), np.percentile(sp, 25), np.median(sp))
def test_spec_highlight_follows_the_half_vector_of_light_and_viewer():
"""spec_field against an independent evaluation: half vector H = normalise(L + V), V the unit vector back to the viewer (-r) in every column, n = (-sl r - la t + z)/|.|, GGX alpha 0.3. Two different lights, random faces: the field must match to 1e-9
(a half vector that ignores the light, or any change of the lobe, shows)."""
rng = np.random.default_rng(5); sl = rng.uniform(-0.3, 1.2, (7, W)); la = rng.uniform(-0.6, 0.6, (7, W)); nrm = np.sqrt(1 + sl * sl + la * la); az = ridge.azrad(); r = np.stack([np.sin(az), np.cos(az), 0 * az], 1); t = np.stack([np.cos(az), -np.sin(az), 0 * az], 1)
for laz, lel in ((45.0, 32.0), (200.0, 15.0)):
Lv = ridge.sun_vec(laz, lel); ndl = (Lv[2] - sl * (r @ Lv)[None] - la * (t @ Lv)[None]) / nrm; got = ridge.spec_field(Lv, sl, la, nrm, ndl)
H = Lv[None, :] - r; H /= np.linalg.norm(H, axis=1, keepdims=True); nh = np.clip((H[:, 2][None] - sl * (r * H).sum(1)[None] - la * (t * H).sum(1)[None]) / nrm, 0, 1)
want = 0.30 * 0.0081 / (1 + nh * nh * (0.09 - 1)) ** 2 * np.clip(ndl, 0, 1)
assert np.abs(got - want).max() < 1e-9, ('spec does not match the half-vector model', laz, np.abs(got - want).max())
def test_fresnel_rim_is_strongest_on_faces_that_lie_away_from_the_viewer():
"""fresnel(sl, nrm) = (1 - n.V)^3 with n.V = sl / |n|: 1 for a flat or back-leaning face (sl <= 0), falling monotonically to ~0 for a face that stands up at the viewer (sl = 5: 0.00001)."""
s = np.array([-1.0, 0.0, 0.2, 0.5, 1.0, 2.0, 5.0]); f = ridge.fresnel(s, np.sqrt(1 + s * s))
assert f[0] == 1.0 and f[1] == 1.0 and (np.diff(f[1:]) < 0).all() and f[-1] < 1e-3, ('fresnel shape', f)
def test_cache_keys_follow_every_input_of_the_cached_array():
"""mosaic-/rays- cache files are keyed by a digest of their inputs. Right answer: identical inputs -> identical key (the preset id does not enter), a changed preset field (eye height, range, position), tile list or baker constant (refraction, detail noise) -> a different key."""
saved = ridge.P
try:
base = dict(id='k1', lat=47.0, lon=11.0, rmax_km=40, eye_m=2000.0); ridge.P = dict(base); k0 = (ridge.rays_key(False), ridge.rays_key(True), ridge.mosaic_key(ridge.P))
ridge.P = dict(base, id='other'); assert (ridge.rays_key(False), ridge.rays_key(True), ridge.mosaic_key(ridge.P)) == k0, 'the preset id changed the key'
for fld, v in (('eye_m', 2010.0), ('rmax_km', 41), ('lon', 11.001), ('lat', 47.001)):
ridge.P = dict(base, **{fld: v}); assert ridge.rays_key(False) != k0[0], ('rays key ignores', fld)
ridge.P = dict(base, lat=47.0, lon=13.5); assert ridge.mosaic_key(ridge.P) != k0[2], 'mosaic key ignores the tile list'
ridge.P = dict(base); old = ridge.KREF; ridge.KREF = 0.2
try: assert ridge.rays_key(False) != k0[0], 'rays key ignores the refraction constant'
finally: ridge.KREF = old
old = ridge.DETAIL; ridge.DETAIL = dict(old, amp_m=30.0)
try: assert ridge.rays_key(True) != k0[1] and ridge.rays_key(False) == k0[0], 'detail key ignores the noise settings (or the plain key follows them)'
finally: ridge.DETAIL = old
finally: ridge.P = saved
def test_paths_come_from_the_environment_with_no_default():
"""Importing the baker without RIDGE_WORK / RIDGE_TILES fails with a message naming the variable (no location is assumed)."""
import subprocess; env = {k: v for k, v in os.environ.items() if not k.startswith('RIDGE_')}
for var in ('RIDGE_WORK', 'RIDGE_TILES'):
e = dict(env); e.update({'RIDGE_WORK': '.', 'RIDGE_TILES': '.'}); del e[var]
r = subprocess.run([sys.executable, '-I', '-c', 'import sys; sys.path.insert(0, %r); import ridge' % SRC], env=e, capture_output=True, text=True); assert r.returncode and var in r.stderr and 'not set' in r.stderr, (var, r.stderr[-200:])
TESTS = [v for k, v in sorted(globals().items()) if k.startswith('test_')]
if __name__ == '__main__':
fails = 0
for t in TESTS:
try: t(); print('PASS', t.__name__)
except Exception as e: fails += 1; print('FAIL', t.__name__, '|', repr(e)[:300])
print('failures:', fails); sys.exit(1 if fails else 0)