veloren_common/terrain/
map.rs

1use super::{
2    NEIGHBOR_DELTA, TERRAIN_CHUNK_BLOCKS_LG, TerrainChunkSize, neighbors, quadratic_nearest_point,
3    river_spline_coeffs, uniform_idx_as_vec2, vec2_as_uniform_idx,
4};
5use crate::vol::RectVolSize;
6use common_base::prof_span;
7use core::{f32, f64, iter, ops::RangeInclusive};
8use vek::*;
9
10/// Base two logarithm of the maximum size of the precomputed world, in meters,
11/// along the x (E/W) and y (N/S) dimensions.
12///
13/// NOTE: Each dimension is guaranteed to be a power of 2, so the logarithm is
14/// exact. This is so that it is possible (at least in theory) for compiler or
15/// runtime optimizations exploiting this are possible.  For example, division
16/// by the chunk size can turn into a bit shift.
17///
18/// NOTE: As an invariant, this value is at least [TERRAIN_CHUNK_BLOCKS_LG].
19///
20/// NOTE: As an invariant, `(1 << [MAX_WORLD_BLOCKS_LG])` fits in an i32.
21///
22/// TODO: Add static assertions for the above invariants.
23///
24/// Currently, we define the maximum to be 19 (corresponding to 2^19 m) for both
25/// directions. This value was derived by backwards reasoning from the following
26/// conservative estimate of the maximum landmass area (using an approximation
27/// of 1024 blocks / km instead of 1000 blocks / km, which will result in an
28/// estimate that is strictly lower than the real landmass):
29///
30/// Max area (km²)
31///     ≌ (2^19 blocks * 1 km / 1024 blocks)^2
32///     = 2^((19 - 10) * 2) km²
33///     = 2^18 km²
34///     = 262,144 km²
35///
36/// which is roughly the same area as the entire United Kingdom, and twice the
37/// horizontal extent of Dwarf Fortress's largest map.  Besides the comparison
38/// to other games without infinite or near-infinite maps (like Dwarf Fortress),
39/// there are other reasons to choose this as a good maximum size:
40///
41/// * It is large enough to include geological features of fairly realistic
42///   scale.  It may be hard to do justice to truly enormous features like the
43///   Amazon River, and natural temperature variation not related to altitude
44///   would probably not produce climate extremes on an Earth-like planet, but
45///   it can comfortably fit enormous river basins, Everest-scale mountains,
46///   large islands and inland lakes, vast forests and deserts, and so on.
47///
48/// * It is large enough that making it from one side of the map to another will
49///   take a *very* long time.  We show this with two examples.  In each
50///   example, travel is either purely horizontal or purely vertical (to
51///   minimize distance traveled) across the whole map, and we assume there are
52///   no obstacles or slopes.
53///
54/// In example 1, a human is walking at the (real-time) speed of the fastest
55/// marathon runners   (around 6 blocks / real-time s).  We assume the human can
56/// maintain this pace indefinitely   without stopping.  Then crossing the map
57/// will take about:
58///
59///   2^19 blocks * 1 real-time s / 6 blocks * 1 real-time min / 60 real-time s
60/// * 1 real-time hr / 60 real-time min * 1 real-time days / 24 hr = 2^19 / 6 /
61///   60 / 60 / 24 real-time days ≌ 1 real-time day.
62///
63/// That's right--it will take a full day of *real* time to cross the map at
64/// an apparent speed of   6 m / s.  Moreover, since in-game time passes at a
65/// rate of 1 in-game min / 1 in-game s, this   would also take *60 days* of
66/// in-game time.
67///
68/// Still though, this is the rate of an ordinary human.  And besides that, if
69/// we instead had a   marathon runner traveling at 6 m / in-game s, it would
70/// take just 1 day of in-game time for   the runner to cross the map, or a mere
71/// 0.4 hr of real time.   To show that this rate of travel is unrealistic (and
72/// perhaps make an eventual argument for   a slower real-time to in-game time
73/// conversion rate), our second example will consist of a   high-speed train
74/// running at 300 km / real-time h (the fastest real-world high speed train
75///  averages under 270 k m / h, with 300 km / h as the designed top speed).
76/// For a train traveling at this apparent speed (in real time), crossing the
77/// map would take:
78///
79///   2^19 blocks * 1 km / 1000 blocks * 1 real-time hr / 300 km
80///   = 2^19 / 1000 / 300 real-time hr
81///   ≌ 1.75 real-time hr
82///
83///   = 2^19 / 1000 / 300 real-time hr * 60 in-game hr / real-time hr
84///     * 1 in-game days / 24 in-game hr
85///   = 2^19 / 1000 / 300 * 60 / 24 in-game days
86///   ≌ 4.37 in-game days
87///
88/// In other words, something faster in real-time than any existing high-speed
89/// train would be   over 4 times slower (in real-time) than our hypothetical
90/// top marathon runner running at 6 m /   s in in-game speed.  This suggests
91/// that the gap between in-game time and real-time is   probably much too large
92/// for most purposes; however, what it definitely shows is that even
93/// extremely fast in-game transport across the world will not trivialize its
94/// size.
95///
96/// It follows that cities or towns of realistic scale, player housing,
97/// fields, and so on, will   all fit comfortably on a map of this size, while
98/// at the same time still being reachable by   non-warping, in-game mechanisms
99/// (such as high-speed transit).  It also provides plenty of   room for mounts
100/// of varying speeds, which can help ensure that players don't feel cramped or
101///   deliberately slowed down by their own speed.
102///
103/// * It is small enough that it is (barely) plausible that we could still
104///   generate maps for a world of this size using detailed and realistic
105///   erosion algorithms.  At 1/4 of this map size along each dimension,
106///   generation currently takes around 5 hours on a good computer, and one
107///   could imagine (since the bottleneck step appears to be roughly O(n)) that
108///   with a smart implementation generation times of under a week could be
109///   achievable.
110///
111/// * The map extends further than the resolution of human eyesight under
112///   Earthlike conditions, even from tall mountains across clear landscapes.
113///   According to one calculation, even from Mt. Everest in the absence of
114///   cloud cover, you could only see for about 339 km before the Earth's
115///   horizon prevented you from seeing further, and other sources suggest that
116///   in practice the limit is closer to 160 km under realistic conditions. This
117///   implies that making the map much larger in a realistic way would require
118///   incorporating curvature, and also implies that any features that cannot
119///   fit on the map would not (under realistic atmospheric conditions) be fully
120///   visible from any point on Earth.  Therefore, even if we cannot represent
121///   features larger than this accurately, nothing should be amiss from a
122///   visual perspective, so this should not significantly impact the player
123///   experience.
124pub const MAX_WORLD_BLOCKS_LG: Vec2<u32> = Vec2 { x: 19, y: 19 };
125
126/// Base two logarithm of a world size, in chunks, per dimension
127/// (each dimension must be a power of 2, so the logarithm is exact).
128///
129/// NOTE: As an invariant, each dimension must be between 0 and
130/// `[MAX_WORLD_BLOCKS_LG] - [TERRAIN_CHUNK_BLOCKS_LG]`.
131///
132/// NOTE: As an invariant, `(1 << ([DEFAULT_WORLD_CHUNKS_LG] +
133/// [TERRAIN_CHUNK_BLOCKS_LG]))` fits in an i32 (derived from the invariant
134/// on [MAX_WORLD_BLOCKS_LG]).
135///
136/// NOTE: As an invariant, each dimension (in chunks) must fit in a i16.
137///
138/// NOTE: As an invariant, the product of dimensions (in chunks) must fit in a
139/// usize.
140///
141/// These invariants are all checked on construction of a `MapSizeLg`.
142#[derive(Clone, Copy, Debug)]
143pub struct MapSizeLg(Vec2<u32>);
144
145impl MapSizeLg {
146    // FIXME: We cannot use is_some() here because it is not currently marked as a
147    // `const fn`.  Since being able to use conditionals in constant expressions has
148    // not technically been stabilized yet, Clippy probably doesn't check for this
149    // case yet.  When it can, or when is_some() is stabilized as a `const fn`,
150    // we should deal with this.
151    /// Construct a new `MapSizeLg`, returning an error if the needed invariants
152    /// do not hold and the vector otherwise.
153    ///
154    /// TODO: In the future, we may use unsafe code to assert to the compiler
155    /// that these invariants indeed hold, safely opening up optimizations
156    /// that might not otherwise be available at runtime.
157    #[inline(always)]
158    #[expect(clippy::result_unit_err)]
159    pub const fn new(map_size_lg: Vec2<u32>) -> Result<Self, ()> {
160        // Assertion on dimensions: must be between
161        // 0 and ([MAX_WORLD_BLOCKS_LG] - [TERRAIN_CHUNK_BLOCKS_LG])
162        let is_le_max = map_size_lg.x <= MAX_WORLD_BLOCKS_LG.x - TERRAIN_CHUNK_BLOCKS_LG
163            && map_size_lg.y <= MAX_WORLD_BLOCKS_LG.y - TERRAIN_CHUNK_BLOCKS_LG;
164        // Assertion on dimensions: chunks must fit in a i16.
165        let chunks_in_range =
166            /* 1u15.checked_shl(map_size_lg.x).is_some() &&
167            1u15.checked_shl(map_size_lg.y).is_some(); */
168            map_size_lg.x <= 15 &&
169            map_size_lg.y <= 15;
170        if is_le_max && chunks_in_range {
171            // Assertion on dimensions: blocks must fit in a i32.
172            let blocks_in_range =
173                /* 1i32.checked_shl(map_size_lg.x + TERRAIN_CHUNK_BLOCKS_LG).is_some() &&
174                1i32.checked_shl(map_size_lg.y + TERRAIN_CHUNK_BLOCKS_LG).is_some(); */
175                map_size_lg.x + TERRAIN_CHUNK_BLOCKS_LG < 32 &&
176                map_size_lg.y + TERRAIN_CHUNK_BLOCKS_LG < 32;
177            // Assertion on dimensions: product of dimensions must fit in a usize.
178            let chunks_product_in_range =
179                1usize.checked_shl(map_size_lg.x + map_size_lg.y).is_some();
180            if blocks_in_range && chunks_product_in_range {
181                // Cleared all invariants.
182                Ok(MapSizeLg(map_size_lg))
183            } else {
184                Err(())
185            }
186        } else {
187            Err(())
188        }
189    }
190
191    #[inline(always)]
192    /// Acquire the `MapSizeLg`'s inner vector.
193    pub const fn vec(self) -> Vec2<u32> { self.0 }
194
195    #[inline(always)]
196    /// Get the size of this map in chunks.
197    pub const fn chunks(self) -> Vec2<u16> { Vec2::new(1 << self.0.x, 1 << self.0.y) }
198
199    /// Get the size of an array of the correct size to hold all chunks.
200    pub const fn chunks_len(self) -> usize { 1 << (self.0.x + self.0.y) }
201
202    #[inline(always)]
203    /// Determine whether a chunk position is in bounds.
204    pub const fn contains_chunk(&self, chunk_key: Vec2<i32>) -> bool {
205        let map_size = self.chunks();
206        chunk_key.x >= 0
207            && chunk_key.y >= 0
208            && chunk_key.x == chunk_key.x & ((map_size.x as i32) - 1)
209            && chunk_key.y == chunk_key.y & ((map_size.y as i32) - 1)
210    }
211}
212
213impl From<MapSizeLg> for Vec2<u32> {
214    #[inline(always)]
215    fn from(size: MapSizeLg) -> Self { size.vec() }
216}
217
218pub struct MapConfig<'a> {
219    /// Base two logarithm of the chunk dimensions of the base map.
220    /// Has no default; set explicitly during initial orthographic projection.
221    pub map_size_lg: MapSizeLg,
222    /// Dimensions of the window being written to.
223    ///
224    /// Defaults to `1 << [MapConfig::map_size_lg]`.
225    pub dimensions: Vec2<usize>,
226    /// x, y, and z of top left of map.
227    ///
228    /// Default x and y are 0.0; no reasonable default for z, so set during
229    /// initial orthographic projection.
230    pub focus: Vec3<f64>,
231    /// Altitude is divided by gain and clamped to [0, 1]; thus, decreasing gain
232    /// makes smaller differences in altitude appear larger.
233    ///
234    /// No reasonable default for z; set during initial orthographic projection.
235    pub gain: f32,
236    /// `fov` is used for shading purposes and refers to how much impact a
237    /// change in the z direction has on the perceived slope relative to the
238    /// same change in x and y.
239    ///
240    /// It is stored as cos θ in the range (0, 1\] where θ is the FOV
241    /// "half-angle" used for perspective projection.  At 1.0, we treat it
242    /// as the limit value for θ = 90 degrees, and use an orthographic
243    /// projection.
244    ///
245    /// Defaults to 1.0.
246    ///
247    /// FIXME: This is a hack that tries to incorrectly implement a variant of
248    /// perspective projection (which generates ∂P/∂x and ∂P/∂y for screen
249    /// coordinate P by using the hyperbolic function \[assuming frustum of
250    /// \[l, r, b, t, n, f\], rh coordinates, and output from -1 to 1 in
251    /// s/t, 0 to 1 in r, and NDC is left-handed \[so visible z ranges from
252    /// -n to -f\]\]):
253    ///
254    /// P.s(x, y, z) = -1 +  2(-n/z x -  l) / ( r -  l)
255    /// P.t(x, y, z) = -1 +  2(-n/z y -  b) / ( t -  b)
256    /// P.r(x, y, z) =  0 + -f(-n/z   -  1) / ( f -  n)
257    ///
258    /// Then arbitrarily using W_e_x = (r - l) as the width of the projected
259    /// image, we have W_e_x = 2 n_e tan θ ⇒ tan Θ = (r - l) / (2n_e), for a
260    /// perspective projection
261    ///
262    /// (where θ is the half-angle of the FOV).
263    ///
264    /// Taking the limit as θ → 90, we show that this degenerates to an
265    /// orthogonal projection:
266    ///
267    /// lim{n → ∞}(-f(-n / z - 1) / (f - n)) = -(z - -n) / (f - n).
268    ///
269    /// (Proof not currently included, but has been formalized for the P.r case
270    /// in Coq-tactic notation; the proof can be added on request, but is
271    /// large and probably not well-suited to Rust documentation).
272    ///
273    /// For this reason, we feel free to store `fov` as cos θ in the range (0,
274    /// 1\].
275    ///
276    /// However, `fov` does not actually work properly yet, so for now we just
277    /// treat it as a visual gimmick.
278    pub fov: f64,
279    /// Scale is like gain, but for x and y rather than z.
280    ///
281    /// Defaults to (1 << world_size_lg).x / dimensions.x (NOTE: fractional, not
282    /// integer, division!).
283    pub scale: f64,
284    /// Vector that indicates which direction light is coming from, if shading
285    /// is turned on.
286    ///
287    /// Right-handed coordinate system: light is going left, down, and
288    /// "backwards" (i.e. on the map, where we translate the y coordinate on
289    /// the world map to z in the coordinate system, the light comes from -y
290    /// on the map and points towards +y on the map).  In a right
291    /// handed coordinate system, the "camera" points towards -z, so positive z
292    /// is backwards "into" the camera.
293    ///
294    /// "In world space the x-axis will be pointing east, the y-axis up and the
295    /// z-axis will be pointing south"
296    ///
297    /// Defaults to (-0.8, -1.0, 0.3).
298    pub light_direction: Vec3<f64>,
299    /// If Some, uses the provided horizon map.
300    ///
301    /// Defaults to None.
302    pub horizons: Option<&'a [(Vec<f32>, Vec<f32>); 2]>,
303    /// If true, only the basement (bedrock) is used for altitude; otherwise,
304    /// the surface is used.
305    ///
306    /// Defaults to false.
307    pub is_basement: bool,
308    /// If true, water is rendered; otherwise, the surface without water is
309    /// rendered, even if it is underwater.
310    ///
311    /// Defaults to true.
312    pub is_water: bool,
313    /// When `is_water` is true, controls whether an ice layer should appear on
314    /// that water.
315    ///
316    /// Defaults to true.
317    pub is_ice: bool,
318    /// If true, 3D lighting and shading are turned on.  Otherwise, a plain
319    /// altitude map is used.
320    ///
321    /// Defaults to true.
322    pub is_shaded: bool,
323    /// If true, the red component of the image is also used for temperature
324    /// (redder is hotter). Defaults to false.
325    pub is_temperature: bool,
326    /// If true, the blue component of the image is also used for humidity
327    /// (bluer is wetter).
328    ///
329    /// Defaults to false.
330    pub is_humidity: bool,
331    /// Record debug information.
332    ///
333    /// Defaults to false.
334    pub is_debug: bool,
335    /// If true, contour lines are drawn on top of the base rbg
336    ///
337    /// Defaults to false.
338    pub is_contours: bool,
339    /// If true, a yellow/terracotta heightmap shading is applied to the
340    /// terrain and water is a faded blue.
341    ///
342    /// Defaults to false
343    pub is_height_map: bool,
344    /// Applies contour lines as well as color modifications
345    ///
346    /// Defaults to false
347    pub is_stylized_topo: bool,
348}
349
350pub const QUADRANTS: usize = 4;
351
352pub struct MapDebug {
353    pub quads: [[u32; QUADRANTS]; QUADRANTS],
354    pub rivers: u32,
355    pub lakes: u32,
356    pub oceans: u32,
357}
358
359/// Connection kind (per edge).  Currently just supports rivers, but may be
360/// extended to support paths or at least one other kind of connection.
361#[derive(Clone, Copy, Debug)]
362pub enum ConnectionKind {
363    /// Connection forms a visible river.
364    River,
365}
366
367/// Map connection (per edge).
368#[derive(Clone, Copy, Debug)]
369pub struct Connection {
370    /// The kind of connection this is (e.g. river or path).
371    pub kind: ConnectionKind,
372    /// Assumed to be the "b" part of a 2d quadratic function.
373    pub spline_derivative: Vec2<f32>,
374    /// Width of the connection.
375    pub width: f32,
376}
377
378/// Per-chunk data the map needs to be able to sample in order to correctly
379/// render.
380#[derive(Clone, Debug)]
381pub struct MapSample {
382    /// the base RGB color for a particular map pixel using the current settings
383    /// (i.e. the color *without* lighting).
384    pub rgb: Rgb<u8>,
385    /// Surface altitude information
386    /// (correctly reflecting settings like is_basement and is_water)
387    pub alt: f64,
388    /// Downhill chunk (may not be meaningful on ocean tiles, or at least edge
389    /// tiles)
390    pub downhill_wpos: Vec2<i32>,
391    /// Connection information about any connections to/from this chunk (e.g.
392    /// rivers).
393    ///
394    /// Connections at each index correspond to the same index in
395    /// NEIGHBOR_DELTA.
396    pub connections: Option<[Option<Connection>; 8]>,
397}
398
399impl MapConfig<'_> {
400    /// Constructs the configuration settings for an orthographic projection of
401    /// a map from the top down, rendering (by default) the complete map to
402    /// an image such that the chunk:pixel ratio is 1:1.
403    ///
404    /// Takes two arguments: the base two logarithm of the horizontal map extent
405    /// (in chunks), and the z bounds of the projection.
406    pub fn orthographic(map_size_lg: MapSizeLg, z_bounds: RangeInclusive<f32>) -> Self {
407        assert!(z_bounds.start() <= z_bounds.end());
408        // NOTE: Safe cast since map_size_lg is restricted by the prior assert.
409        let dimensions = map_size_lg.chunks().map(usize::from);
410        Self {
411            map_size_lg,
412            dimensions,
413            focus: Vec3::new(0.0, 0.0, f64::from(*z_bounds.start())),
414            gain: z_bounds.end() - z_bounds.start(),
415            fov: 1.0,
416            scale: 1.0,
417            light_direction: Vec3::new(-1.2, -1.0, 0.8),
418            horizons: None,
419
420            is_basement: false,
421            is_water: true,
422            is_ice: true,
423            is_shaded: true,
424            is_temperature: false,
425            is_humidity: false,
426            is_debug: false,
427            is_contours: false,
428            is_height_map: false,
429            is_stylized_topo: false,
430        }
431    }
432
433    /// Get the base 2 logarithm of the underlying map size.
434    pub fn map_size_lg(&self) -> MapSizeLg { self.map_size_lg }
435
436    /// Generates a map image using the specified settings.  Note that it will
437    /// write from left to write from (0, 0) to dimensions - 1, inclusive,
438    /// with 4 1-byte color components provided as (r, g, b, a).  It is up
439    /// to the caller to provide a function that translates this information
440    /// into the correct format for a buffer and writes to it.
441    ///
442    /// sample_pos is a function that, given a chunk position, returns enough
443    /// information about the chunk to attempt to render it on the map.
444    /// When in doubt, try using `MapConfig::sample_pos` for this.
445    ///
446    /// sample_wpos is a simple function that, given a *column* position,
447    /// returns the approximate altitude at that column.  When in doubt, try
448    /// using `MapConfig::sample_wpos` for this.
449    pub fn generate(
450        &self,
451        sample_pos: impl Fn(Vec2<i32>) -> MapSample,
452        sample_wpos: impl Fn(Vec2<i32>) -> f32,
453        mut write_pixel: impl FnMut(Vec2<usize>, (u8, u8, u8, u8)),
454    ) -> MapDebug {
455        prof_span!("MapConfig::generate");
456        let MapConfig {
457            map_size_lg,
458            dimensions,
459            focus,
460            gain,
461            fov,
462            scale,
463            light_direction,
464            horizons,
465            is_shaded,
466            is_stylized_topo,
467            // is_debug,
468            ..
469        } = *self;
470
471        let light_direction = Vec3::new(
472            light_direction.x,
473            light_direction.y,
474            0.0, // we currently ignore light_direction.z.
475        );
476        let light_shadow_dir = usize::from(light_direction.x < 0.0);
477        let horizon_map = horizons.map(|horizons| &horizons[light_shadow_dir]);
478        let light = light_direction.normalized();
479        let /*mut */quads = [[0u32; QUADRANTS]; QUADRANTS];
480        let /*mut */rivers = 0u32;
481        let /*mut */lakes = 0u32;
482        let /*mut */oceans = 0u32;
483
484        let focus_rect = Vec2::from(focus);
485
486        let chunk_size = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
487
488        /* // NOTE: Asserting this to enable LLVM optimizations.  Ideally we should come up
489        // with a principled way to do this (especially one with no runtime
490        // cost).
491        assert!(
492            map_size_lg
493                .vec()
494                .cmple(&(MAX_WORLD_BLOCKS_LG - TERRAIN_CHUNK_BLOCKS_LG))
495                .reduce_and()
496        ); */
497        let world_size = map_size_lg.chunks();
498
499        (0..dimensions.y * dimensions.x).for_each(|chunk_idx| {
500            let i = chunk_idx % dimensions.x;
501            let j = chunk_idx / dimensions.x;
502
503            let wposf = focus_rect + Vec2::new(i as f64, j as f64) * scale;
504            let pos = wposf.map(|e: f64| e as i32);
505            let wposf = wposf * chunk_size;
506
507            let chunk_idx = if pos.reduce_partial_min() >= 0
508                && pos.x < world_size.x as i32
509                && pos.y < world_size.y as i32
510            {
511                Some(vec2_as_uniform_idx(map_size_lg, pos))
512            } else {
513                None
514            };
515
516            let MapSample {
517                rgb,
518                alt,
519                downhill_wpos,
520                ..
521            } = sample_pos(pos);
522
523            let alt = alt as f32;
524            let wposi = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32);
525            let mut rgb = rgb.map(|e| e as f64 / 255.0);
526
527            // Material properties:
528            //
529            // For each material in the scene,
530            //  k_s = (RGB) specular reflection constant
531            let mut k_s = Rgb::new(1.0, 1.0, 1.0);
532            //  k_d = (RGB) diffuse reflection constant
533            let mut k_d = rgb;
534            //  k_a = (RGB) ambient reflection constant
535            let mut k_a = rgb;
536            //  α = (per-material) shininess constant
537            let mut alpha = 4.0; // 4.0;
538
539            // Compute connections
540            let mut has_river = false;
541            // NOTE: consider replacing neighbors with local_cells, since it is more
542            // accurate (though I'm not sure if it can matter for these
543            // purposes).
544            chunk_idx
545                .into_iter()
546                .flat_map(|chunk_idx| {
547                    neighbors(map_size_lg, chunk_idx).chain(iter::once(chunk_idx))
548                })
549                .for_each(|neighbor_posi| {
550                    let neighbor_pos = uniform_idx_as_vec2(map_size_lg, neighbor_posi);
551                    let neighbor_wpos = neighbor_pos.map(|e| e as f64) * chunk_size;
552                    let MapSample { connections, .. } = sample_pos(neighbor_pos);
553                    NEIGHBOR_DELTA
554                        .iter()
555                        .zip(connections.iter().flatten())
556                        .for_each(|(&delta, connection)| {
557                            let connection = if let Some(connection) = connection {
558                                connection
559                            } else {
560                                return;
561                            };
562                            let downhill_wpos = neighbor_wpos
563                                + Vec2::from(delta).map(|e: i32| e as f64) * chunk_size;
564                            let coeffs = river_spline_coeffs(
565                                neighbor_wpos,
566                                connection.spline_derivative,
567                                downhill_wpos,
568                            );
569                            let (_t, _pt, dist) = if let Some((t, pt, dist)) =
570                                quadratic_nearest_point(
571                                    &coeffs,
572                                    wposf,
573                                    Vec2::new(neighbor_wpos, downhill_wpos),
574                                ) {
575                                (t, pt, dist)
576                            } else {
577                                let ndist = wposf.distance_squared(neighbor_wpos);
578                                let ddist = wposf.distance_squared(downhill_wpos);
579                                if ndist <= ddist {
580                                    (0.0, neighbor_wpos, ndist)
581                                } else {
582                                    (1.0, downhill_wpos, ddist)
583                                }
584                            };
585                            let connection_dist =
586                                (dist.sqrt() - (connection.width as f64 * 0.5).max(1.0)).max(0.0);
587                            if connection_dist == 0.0 {
588                                match connection.kind {
589                                    ConnectionKind::River => {
590                                        has_river = true;
591                                    },
592                                }
593                            }
594                        });
595                });
596
597            // Color in connections.
598            let water_color_factor = 2.0;
599            let g_water = 32.0 * water_color_factor;
600            let b_water = 64.0 * water_color_factor;
601            if has_river {
602                // Rudimentary ice check
603                if !rgb.map(|e| e > 0.35).reduce_and() {
604                    let water_rgb = Rgb::new(0, ((g_water) * 1.0) as u8, ((b_water) * 1.0) as u8)
605                        .map(|e| e as f64 / 255.0);
606                    rgb = water_rgb;
607                    k_s = Rgb::new(1.0, 1.0, 1.0);
608                    k_d = water_rgb;
609                    k_a = water_rgb;
610                    alpha = 0.255;
611                }
612            }
613
614            let downhill_alt = sample_wpos(downhill_wpos);
615            let cross_pos = wposi
616                + ((downhill_wpos - wposi)
617                    .map(|e| e as f32)
618                    .rotated_z(f32::consts::FRAC_PI_2)
619                    .map(|e| e as i32));
620            let cross_alt = sample_wpos(cross_pos);
621            // TODO: Fix use of fov to match proper perspective projection, as described in
622            // the doc comment.
623            // Pointing downhill, forward
624            // (index--note that (0,0,1) is backward right-handed)
625            let forward_vec = Vec3::new(
626                (downhill_wpos.x - wposi.x) as f64,
627                ((downhill_alt - alt) * gain) as f64 * fov,
628                (downhill_wpos.y - wposi.y) as f64,
629            );
630            // Pointing 90 degrees left (in horizontal xy) of downhill, up
631            // (middle--note that (1,0,0), 90 degrees CCW backward, is right right-handed)
632            let up_vec = Vec3::new(
633                (cross_pos.x - wposi.x) as f64,
634                ((cross_alt - alt) * gain) as f64 * fov,
635                (cross_pos.y - wposi.y) as f64,
636            );
637            // let surface_normal = Vec3::new(fov* (f.y * u.z - f.z * u.y), -(f.x * u.z -
638            // f.z * u.x), fov* (f.x * u.y - f.y * u.x)).normalized();
639            // Then cross points "to the right" (upwards) on a right-handed coordinate
640            // system. (right-handed coordinate system means (0, 0, 1.0) is
641            // "forward" into the screen).
642            let surface_normal = forward_vec.cross(up_vec).normalized();
643
644            // TODO: Figure out if we can reimplement debugging.
645            /* if is_debug {
646                let quad =
647                    |x: f32| ((x as f64 * QUADRANTS as f64).floor() as usize).min(QUADRANTS - 1);
648                if river_kind.is_none() || humidity != 0.0 {
649                    quads[quad(humidity)][quad(temperature)] += 1;
650                }
651                match river_kind {
652                    Some(RiverKind::River { .. }) => {
653                        rivers += 1;
654                    },
655                    Some(RiverKind::Lake { .. }) => {
656                        lakes += 1;
657                    },
658                    Some(RiverKind::Ocean { .. }) => {
659                        oceans += 1;
660                    },
661                    None => {},
662                }
663            } */
664
665            let rgb = if is_shaded {
666                let shade_frac = horizon_map
667                    .and_then(|(angles, heights)| {
668                        chunk_idx
669                            .and_then(|chunk_idx| angles.get(chunk_idx))
670                            .map(|&e| (e as f64, heights))
671                    })
672                    .and_then(|(e, heights)| {
673                        chunk_idx
674                            .and_then(|chunk_idx| heights.get(chunk_idx))
675                            .map(|&f| (e, f as f64))
676                    })
677                    .map(|(angle, height)| {
678                        let w = 0.1;
679                        let height = (height - f64::from(alt * gain)).max(0.0);
680                        if angle != 0.0 && light_direction.x != 0.0 && height != 0.0 {
681                            let deltax = height / angle;
682                            let lighty = (light_direction.y / light_direction.x * deltax).abs();
683                            let deltay = lighty - height;
684                            let s = (deltay / deltax / w).clamp(0.0, 1.0);
685                            // Smoothstep
686                            s * s * (3.0 - 2.0 * s)
687                        } else {
688                            1.0
689                        }
690                    })
691                    .unwrap_or(1.0);
692
693                // Phong reflection model with shadows:
694                //
695                // I_p = k_a i_a + shadow * Σ {m ∈ lights} (k_d (L_m ⋅ N) i_m,d + k_s (R_m ⋅
696                // V)^α i_m,s)
697                //
698                // where for the whole scene,
699                //  i_a = (RGB) intensity of ambient lighting component
700                let i_a = if is_stylized_topo {
701                    Rgb::new(0.4, 0.4, 0.4)
702                } else {
703                    Rgb::new(0.1, 0.1, 0.1)
704                };
705                //  V = direction pointing towards the viewer (e.g. virtual camera).
706                let v = Vec3::new(0.0, 0.0, -1.0).normalized();
707
708                // for each light m,
709                //  i_m,d = (RGB) intensity of diffuse component of light source m
710                let i_m_d = Rgb::new(1.0, 1.0, 1.0);
711                //  i_m,s = (RGB) intensity of specular component of light source m
712                let i_m_s = Rgb::new(0.45, 0.45, 0.45);
713
714                // for each light m and point p,
715                //  L_m = (normalized) direction vector from point on surface to light source m
716                let l_m = light;
717                //  N = (normalized) normal at this point on the surface,
718                let n = surface_normal;
719                //  R_m = (normalized) direction a perfectly reflected ray of light from m would
720                // take from point p      = 2(L_m ⋅ N)N - L_m
721                let r_m = (-l_m).reflected(n); // 2 * (l_m.dot(n)) * n - l_m;
722                //
723                // and for each point p in the scene,
724                //  shadow = computed shadow factor at point p
725                // FIXME: Should really just be shade_frac, but with only ambient light we lose
726                // all local lighting detail... some sort of global illumination (e.g.
727                // radiosity) is of course the "right" solution, but maybe we can find
728                // something cheaper?
729                let shadow = 0.2 + 0.8 * shade_frac;
730
731                let lambertian = l_m.dot(n).max(0.0);
732                let spec_angle = r_m.dot(v).max(0.0);
733
734                let ambient = k_a * i_a;
735                let diffuse = k_d * lambertian * i_m_d;
736                let specular = k_s * spec_angle.powf(alpha) * i_m_s;
737                (ambient + shadow * (diffuse + specular)).map(|e| e.min(1.0))
738            } else {
739                rgb
740            }
741            .map(|e| (e * 255.0) as u8);
742
743            let rgba = (rgb.r, rgb.g, rgb.b, 255);
744            write_pixel(Vec2::new(i, j), rgba);
745        });
746
747        MapDebug {
748            quads,
749            rivers,
750            lakes,
751            oceans,
752        }
753    }
754}