diff --git a/scripts/ridge-bake/__pycache__/ridge.cpython-313.pyc b/scripts/ridge-bake/__pycache__/ridge.cpython-313.pyc new file mode 100644 index 00000000..0a31c5c6 Binary files /dev/null and b/scripts/ridge-bake/__pycache__/ridge.cpython-313.pyc differ diff --git a/scripts/ridge-bake/places.json b/scripts/ridge-bake/places.json new file mode 100644 index 00000000..60da1825 --- /dev/null +++ b/scripts/ridge-bake/places.json @@ -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 + ] + } + } + ] +} \ No newline at end of file diff --git a/scripts/ridge-bake/ridge.py b/scripts/ridge-bake/ridge.py new file mode 100644 index 00000000..07b5c94e --- /dev/null +++ b/scripts/ridge-bake/ridge.py @@ -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 diff --git a/scripts/ridge-bake/sun.py b/scripts/ridge-bake/sun.py new file mode 100644 index 00000000..f1eaa182 --- /dev/null +++ b/scripts/ridge-bake/sun.py @@ -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 diff --git a/scripts/ridge-bake/tests/mutate.py b/scripts/ridge-bake/tests/mutate.py new file mode 100644 index 00000000..19c0fc66 --- /dev/null +++ b/scripts/ridge-bake/tests/mutate.py @@ -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) diff --git a/scripts/ridge-bake/tests/run.sh b/scripts/ridge-bake/tests/run.sh new file mode 100755 index 00000000..f97a126c --- /dev/null +++ b/scripts/ridge-bake/tests/run.sh @@ -0,0 +1,7 @@ +#!/bin/sh +# Usage: RIDGE_PY= RIDGE_TILES= 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}" diff --git a/scripts/ridge-bake/tests/test_ridge.py b/scripts/ridge-bake/tests/test_ridge.py new file mode 100644 index 00000000..49e6d34f --- /dev/null +++ b/scripts/ridge-bake/tests/test_ridge.py @@ -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)