Adding shadows.

Cargo.lock

@@ -3321,6 +3321,7 @@ dependencies = [
 name = "veloren-world"
 version = "0.5.0"
 dependencies = [
+ "approx 0.1.1 (registry+https://github.com/rust-lang/crates.io-index)",
  "arr_macro 0.1.3 (registry+https://github.com/rust-lang/crates.io-index)",
  "bincode 1.2.0 (registry+https://github.com/rust-lang/crates.io-index)",
  "bitvec 0.15.2 (registry+https://github.com/rust-lang/crates.io-index)",

world/Cargo.toml

@@ -5,6 +5,7 @@ authors = ["Joshua Barretto <joshua.s.barretto@gmail.com>"]
 edition = "2018"
 
 [dependencies]
+approx = "0.1.1"
 bincode = "1.2.0"
 common = { package = "veloren-common", path = "../common" }
 bitvec = "0.15.2"

world/examples/water.rs

@@ -1,8 +1,12 @@
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
+use rayon::prelude::*;
 use std::{f64, io::Write, path::PathBuf, time::SystemTime};
 use vek::*;
 use veloren_world::{
-    sim::{self, MapConfig, MapDebug, WorldOpts, WORLD_SIZE},
+    sim::{
+        self, get_shadows, uniform_idx_as_vec2, Alt, MapConfig, MapDebug, WorldOpts, WORLD_SIZE,
+    },
+    util::Sampler,
     World, CONFIG,
 };
 
@@ -18,26 +22,64 @@ fn main() {
     let map_file =
         // "map_1575990726223.bin";
         // "map_1575987666972.bin";
-        "map_1576046079066.bin";
+        // "map_1576046079066.bin";
+        "map_1579539133272.bin";
     let mut _map_file = PathBuf::from("./maps");
     _map_file.push(map_file);
 
     let world = World::generate(5284, WorldOpts {
         seed_elements: false,
-        // world_file: sim::FileOpts::Load(_map_file),
+        // world_file: sim::FileOpts::LoadAsset(veloren_world::sim::DEFAULT_WORLD_MAP.into()),
-        world_file: sim::FileOpts::Save,
+        world_file: sim::FileOpts::Load(_map_file),
+        // world_file: sim::FileOpts::Save,
         ..WorldOpts::default()
     });
+    log::info!("Sampling data...");
+    let sampler = world.sim();
+
+    let samples_data = {
+        let column_sample = world.sample_columns();
+        (0..WORLD_SIZE.product())
+            .into_par_iter()
+            .map(|posi| {
+                column_sample
+                    .get(uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+            })
+            .collect::<Vec<_>>()
+            .into_boxed_slice()
+    };
+
+    let refresh_shadows = |light_direction: Vec3<f64>, lgain, scale, is_basement, is_water| {
+        get_shadows(
+        //Vec3::new(-0.8, 0.3, /*-1.0*/-(1.0 / TerrainChunkSize::RECT_SIZE.x as Alt)),
+        Vec3::new(light_direction.x, light_direction.z, light_direction.y/* / lgain*/),
+        lgain,
+        TerrainChunkSize::RECT_SIZE.x as f64/* * scale*/,
+        TerrainChunkSize::RECT_SIZE.y as f64/* * scale*/,
+        Aabr {
+            min: Vec2::new(0.0, 0.0), // focus.into(),
+            max: WORLD_SIZE.map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64)/* * scale*//* + focus.into() */,
+        },
+        CONFIG.sea_level as f64, // focus.z,
+        (CONFIG.sea_level + sampler.max_height) as f64, // (focus.z + self.max_height) as Alt,
+        |posi| {
+            let sample = sampler.get(uniform_idx_as_vec2(posi)).unwrap();
+            if is_basement {
+                sample.alt as f64
+            } else {
+                sample.basement as f64
+            }.max(if is_water { sample.water_alt as f64 } else { -f64::INFINITY })
+        },
+    ).ok()
+    };
 
     let mut win =
         minifb::Window::new("World Viewer", W, H, minifb::WindowOptions::default()).unwrap();
 
-    let sampler = world.sim();
-
     let mut focus = Vec3::new(0.0, 0.0, CONFIG.sea_level as f64);
     // Altitude is divided by gain and clamped to [0, 1]; thus, decreasing gain
     // makes smaller differences in altitude appear larger.
-    let mut gain = CONFIG.mountain_scale;
+    let mut gain = /*CONFIG.mountain_scale*/sampler.max_height;
     // The Z component during normal calculations is multiplied by gain; thus,
     let mut lgain = 1.0;
     let mut scale = WORLD_SIZE.x as f64 / W as f64;
@@ -50,7 +92,7 @@ fn main() {
     //
     // "In world space the x-axis will be pointing east, the y-axis up and the
     // z-axis will be pointing south"
-    let mut light_direction = Vec3::new(-0.8, -1.0, 0.3);
+    let mut light_direction = Vec3::new(-/*0.8*/1.3, -1.0, 0.3);
 
     let mut is_basement = false;
     let mut is_water = true;
@@ -58,6 +100,9 @@ fn main() {
     let mut is_temperature = true;
     let mut is_humidity = true;
 
+    let mut shadows = None; //refresh_shadows(light_direction, lgain, scale, is_basement, is_water);
+    let mut samples = None;
+
     while win.is_open() {
         let config = MapConfig {
             dimensions: Vec2::new(W, H),
@@ -66,12 +111,15 @@ fn main() {
             lgain,
             scale,
             light_direction,
+            shadows: shadows.as_deref(),
+            samples,
 
             is_basement,
             is_water,
             is_shaded,
             is_temperature,
             is_humidity,
+            // is_sampled,
             is_debug: true,
         };
 
@@ -161,6 +209,9 @@ fn main() {
         let is_camera = win.is_key_down(minifb::Key::C);
         if win.is_key_down(minifb::Key::B) {
             is_basement ^= true;
+            shadows = shadows.and_then(|_| {
+                refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+            });
         }
         if win.is_key_down(minifb::Key::H) {
             is_humidity ^= true;
@@ -172,11 +223,25 @@ fn main() {
             is_water ^= true;
         }
         if win.is_key_down(minifb::Key::L) {
+            if is_camera {
+                shadows = match shadows {
+                    Some(_) => None,
+                    None => refresh_shadows(light_direction, lgain, scale, is_basement, is_water),
+                };
+            } else {
                 is_shaded ^= true;
             }
+        }
+        if win.is_key_down(minifb::Key::M) {
+            samples = samples.xor(Some(&*samples_data));
+            // is_sampled ^= true;
+        }
         if win.is_key_down(minifb::Key::W) {
             if is_camera {
                 light_direction.z -= lspd;
+                shadows = shadows.and_then(|_| {
+                    refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+                });
             } else {
                 focus.y -= spd * scale;
             }
@@ -184,6 +249,9 @@ fn main() {
         if win.is_key_down(minifb::Key::A) {
             if is_camera {
                 light_direction.x -= lspd;
+                shadows = shadows.and_then(|_| {
+                    refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+                });
             } else {
                 focus.x -= spd * scale;
             }
@@ -191,6 +259,9 @@ fn main() {
         if win.is_key_down(minifb::Key::S) {
             if is_camera {
                 light_direction.z += lspd;
+                shadows = shadows.and_then(|_| {
+                    refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+                });
             } else {
                 focus.y += spd * scale;
             }
@@ -198,6 +269,9 @@ fn main() {
         if win.is_key_down(minifb::Key::D) {
             if is_camera {
                 light_direction.x += lspd;
+                shadows = shadows.and_then(|_| {
+                    refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+                });
             } else {
                 focus.x += spd * scale;
             }
@@ -206,6 +280,9 @@ fn main() {
             if is_camera {
                 if (lgain * 2.0).is_normal() {
                     lgain *= 2.0;
+                    shadows = shadows.and_then(|_| {
+                        refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+                    });
                 }
             } else {
                 gain += 64.0;
@@ -215,6 +292,9 @@ fn main() {
             if is_camera {
                 if (lgain / 2.0).is_normal() {
                     lgain /= 2.0;
+                    shadows = shadows.and_then(|_| {
+                        refresh_shadows(light_direction, lgain, scale, is_basement, is_water)
+                    });
                 }
             } else {
                 gain = (gain - 64.0).max(64.0);
@@ -226,6 +306,8 @@ fn main() {
             } else {
                 if (scale * 2.0).is_normal() {
                     scale *= 2.0;
+                    // shadows = refresh_shadows(light_direction, lgain, scale,
+                    // is_basement);
                 }
             }
         }
@@ -235,6 +317,8 @@ fn main() {
             } else {
                 if (scale / 2.0).is_normal() {
                     scale /= 2.0;
+                    // shadows = refresh_shadows(light_direction, lgain, scale,
+                    // is_basement);
                 }
             }
         }

world/src/config.rs

@@ -64,5 +64,5 @@ pub const CONFIG: Config = Config {
     river_roughness: 0.06125,
     river_max_width: 2.0,
     river_min_height: 0.25,
-    river_width_to_depth: 1.0,
+    river_width_to_depth: 8.0,
 };

world/src/sim/map.rs

@@ -1,12 +1,13 @@
 use crate::{
-    sim::{RiverKind, WorldSim, WORLD_SIZE},
+    column::ColumnSample,
+    sim::{vec2_as_uniform_idx, Alt, RiverKind, WorldSim, WORLD_SIZE},
     CONFIG,
 };
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
 use std::{f32, f64};
 use vek::*;
 
-pub struct MapConfig {
+pub struct MapConfig<'a> {
     /// Dimensions of the window being written to.  Defaults to WORLD_SIZE.
     pub dimensions: Vec2<usize>,
     /// x, y, and z of top left of map (defaults to (0.0, 0.0,
@@ -43,6 +44,14 @@ pub struct MapConfig {
     ///
     /// Defaults to (-0.8, -1.0, 0.3).
     pub light_direction: Vec3<f64>,
+    /// If Some, uses the provided shadow map.
+    ///
+    /// Defaults to None.
+    pub shadows: Option<&'a [Alt]>,
+    /// If Some, uses the provided column samples to determine surface color.
+    ///
+    /// Defaults to None.
+    pub samples: Option<&'a [Option<ColumnSample<'a>>]>,
     /// If true, only the basement (bedrock) is used for altitude; otherwise,
     /// the surface is used.
     ///
@@ -81,7 +90,7 @@ pub struct MapDebug {
     pub oceans: u32,
 }
 
-impl Default for MapConfig {
+impl<'a> Default for MapConfig<'a> {
     fn default() -> Self {
         let dimensions = WORLD_SIZE;
         Self {
@@ -90,7 +99,9 @@ impl Default for MapConfig {
             gain: CONFIG.mountain_scale,
             lgain: TerrainChunkSize::RECT_SIZE.x as f64,
             scale: WORLD_SIZE.x as f64 / dimensions.x as f64,
-            light_direction: Vec3::new(-0.8, -1.0, 0.3),
+            light_direction: Vec3::new(-1.2, -1.0, 0.8),
+            shadows: None,
+            samples: None,
 
             is_basement: false,
             is_water: true,
@@ -102,7 +113,7 @@ impl Default for MapConfig {
     }
 }
 
-impl MapConfig {
+impl<'a> MapConfig<'a> {
     /// Generates a map image using the specified settings.  Note that it will
     /// write from left to write from (0, 0) to dimensions - 1, inclusive,
     /// with 4 1-byte color components provided as (r, g, b, a).  It is up
@@ -120,6 +131,8 @@ impl MapConfig {
             lgain,
             scale,
             light_direction,
+            shadows,
+            samples,
 
             is_basement,
             is_water,
@@ -129,11 +142,15 @@ impl MapConfig {
             is_debug,
         } = *self;
 
+        let light_direction = Vec3::new(light_direction.x, light_direction.y, light_direction.z);
+        // let light_direction = Vec3::new(light_direction.x * lgain, light_direction.y,
+        // light_direction.z * lgain);
         let light = light_direction.normalized();
         let mut quads = [[0u32; QUADRANTS]; QUADRANTS];
         let mut rivers = 0u32;
         let mut lakes = 0u32;
         let mut oceans = 0u32;
+        // let column_sample = ColumnGen::new(sampler);
 
         let focus_rect = Vec2::from(focus);
         let true_sea_level = (CONFIG.sea_level as f64 - focus.z) / gain as f64;
@@ -147,11 +164,20 @@ impl MapConfig {
                 let pos =
                     (focus_rect + Vec2::new(i as f64, j as f64) * scale).map(|e: f64| e as i32);
 
-                let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) =
+                let (
-                    sampler
+                    chunk_idx,
+                    alt,
+                    basement,
+                    water_alt,
+                    humidity,
+                    temperature,
+                    downhill,
+                    river_kind,
+                ) = sampler
                     .get(pos)
                     .map(|sample| {
                         (
+                            Some(vec2_as_uniform_idx(pos)),
                             sample.alt,
                             sample.basement,
                             sample.water_alt,
@@ -162,6 +188,7 @@ impl MapConfig {
                         )
                     })
                     .unwrap_or((
+                        None,
                         CONFIG.sea_level,
                         CONFIG.sea_level,
                         CONFIG.sea_level,
@@ -173,6 +200,30 @@ impl MapConfig {
                 let humidity = humidity.min(1.0).max(0.0);
                 let temperature = temperature.min(1.0).max(-1.0) * 0.5 + 0.5;
                 let pos = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32);
+                let column_rgb = samples
+                    .and_then(|samples| {
+                        chunk_idx
+                            .and_then(|chunk_idx| samples.get(chunk_idx))
+                            .map(Option::as_ref)
+                            .flatten()
+                    })
+                    .map(|sample| {
+                        if is_basement {
+                            sample.stone_col.map(|e| e as f64 / 255.0)
+                        } else {
+                            sample.surface_color.map(|e| e as f64)
+                        }
+                    });
+                /*let column_rgb = if is_sampled {
+                    column_sample.get(pos)
+                        .map(|sample| if is_basement {
+                            sample.stone_col.map(|e| e as f64 / 255.0)
+                        } else {
+                            sample.surface_color.map(|e| e as f64)
+                        })
+                } else {
+                    None
+                };*/
                 let downhill_pos = (downhill
                     .map(|downhill_pos| downhill_pos)
                     .unwrap_or(pos + TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
@@ -206,6 +257,8 @@ impl MapConfig {
                     (cross_alt - alt) as f64 * lgain,
                     (cross_pos.y - pos.y) as f64,
                 );
+                // let surface_normal = Vec3::new(lgain * (f.y * u.z - f.z * u.y), -(f.x * u.z -
+                // f.z * u.x), lgain * (f.x * u.y - f.y * u.x)).normalized();
                 // Then cross points "to the right" (upwards) on a right-handed coordinate
                 // system. (right-handed coordinate system means (0, 0, 1.0) is
                 // "forward" into the screen).
@@ -239,29 +292,84 @@ impl MapConfig {
                     }
                 }
 
+                let shade_frac = shadows
+                    .and_then(|shadows| chunk_idx.and_then(|chunk_idx| shadows.get(chunk_idx)))
+                    .copied()
+                    .map(|e| e as f64)
+                    .unwrap_or(alt);
                 let water_color_factor = 2.0;
-                let g_water = 32.0 * water_color_factor;
+                let g_water = 32.0
-                let b_water = 64.0 * water_color_factor;
+                    * water_color_factor
+                    * if is_shaded {
+                        0.2 + shade_frac * 0.8
+                    } else {
+                        1.0
+                    };
+                let b_water = 64.0
+                    * water_color_factor
+                    * if is_shaded {
+                        0.2 + shade_frac * 0.8
+                    } else {
+                        1.0
+                    };
+                let column_rgb = column_rgb.unwrap_or(Rgb::new(
+                    if is_shaded { shade_frac * 0.6 } else { alt },
+                    if is_shaded {
+                        0.4 + (shade_frac * 0.6)
+                    } else {
+                        alt
+                    },
+                    if is_shaded { shade_frac * 0.6 } else { alt },
+                ));
                 let rgba = match (river_kind, (is_water, true_alt >= true_sea_level)) {
                     (_, (false, _)) | (None, (_, true)) => {
                         let (r, g, b) = (
-                            (if is_shaded { alt } else { alt }
+                            (column_rgb.r/*if is_shaded { shade_frac * 0.6  } else { alt }*/
                                 * if is_temperature {
                                     temperature as f64
                                 } else if is_shaded {
+                                    if samples.is_some() {
+                                        // column_rgb.r
+                                        0.2 + shade_frac * 0.8
+                                    } else {
+                                        shade_frac * 0.6
+                                    }
+                                } else {
+                                    if samples.is_some() {
                                         alt
                                     } else {
                                         0.0
+                                    }
                                 })
                             .sqrt(),
-                            if is_shaded { 0.4 + (alt * 0.6) } else { alt },
+                            (column_rgb.g
-                            (if is_shaded { alt } else { alt }
+                                * if is_shaded {
+                                    if samples.is_some() {
+                                        // column_rgb.g
+                                        0.2 + shade_frac * 0.8
+                                    } else {
+                                        0.4 + shade_frac * 0.6
+                                    }
+                                } else {
+                                    alt
+                                })
+                            .sqrt(),
+                            (column_rgb.b/*if is_shaded { shade_frac * 0.6 } else { alt }*/
                                 * if is_humidity {
                                     humidity as f64
                                 } else if is_shaded {
+                                    if samples.is_some() {
+                                        // column_rgb.b
+                                        0.2 + shade_frac * 0.8
+                                    } else {
+                                        shade_frac * 0.6
+                                    }
+                                } else {
+                                    if samples.is_some() {
                                         alt
                                     } else {
                                         0.0
+                                    }
                                 })
                             .sqrt(),
                         );
@@ -281,15 +389,39 @@ impl MapConfig {
                     ),
                     (Some(RiverKind::River { .. }), _) => (
                         0,
-                        g_water as u8 + (alt * (127.0 - g_water)) as u8,
+                        g_water as u8
-                        b_water as u8 + (alt * (255.0 - b_water)) as u8,
+                            + (if is_shaded {
+                                0.2 + shade_frac * 0.8
+                            } else {
+                                alt
+                            } * (127.0 - g_water)) as u8,
+                        b_water as u8
+                            + (if is_shaded {
+                                0.2 + shade_frac * 0.8
+                            } else {
+                                alt
+                            } * (255.0 - b_water)) as u8,
                         255,
                     ),
                     (None, _) | (Some(RiverKind::Lake { .. }), _) => (
                         0,
-                        (((g_water + water_alt * (127.0 - 32.0)) + (-water_depth * g_water)) * 1.0)
+                        (((g_water
-                            as u8,
+                            + if is_shaded {
-                        (((b_water + water_alt * (255.0 - b_water)) + (-water_depth * b_water))
+                                0.2 + shade_frac * 0.8
+                            } else {
+                                1.0
+                            } * water_alt
+                                * (127.0 - g_water))
+                            + (-water_depth * g_water))
+                            * 1.0) as u8,
+                        (((b_water
+                            + if is_shaded {
+                                0.2 + shade_frac * 0.8
+                            } else {
+                                1.0
+                            } * water_alt
+                                * (255.0 - b_water))
+                            + (-water_depth * b_water))
                             * 1.0) as u8,
                         255,
                     ),

world/src/sim/mod.rs

@@ -17,7 +17,7 @@ pub use self::{
     map::{MapConfig, MapDebug},
     settlement::Settlement,
     util::{
-        cdf_irwin_hall, downhill, get_oceans, local_cells, map_edge_factor, neighbors,
+        cdf_irwin_hall, downhill, get_oceans, get_shadows, local_cells, map_edge_factor, neighbors,
         uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, InverseCdf, ScaleBias,
         NEIGHBOR_DELTA,
     },
@@ -1101,9 +1101,9 @@ impl WorldSim {
             )
         };
         let flux_old = get_multi_drainage(&mstack, &mrec, &*mwrec, boundary_len);
-        let flux_rivers = get_drainage(&water_alt_pos, &dh, boundary_len);
+        // let flux_rivers = get_drainage(&water_alt_pos, &dh, boundary_len);
         // TODO: Make rivers work with multi-direction flux as well.
-        // let flux_rivers = flux_old.clone();
+        let flux_rivers = flux_old.clone();
 
         let water_height_initial = |chunk_idx| {
             let indirection_idx = indirection[chunk_idx];
@@ -1196,6 +1196,21 @@ impl WorldSim {
             None => false,
         };
 
+        /* // Build a shadow map.
+        let shadows = get_shadows(
+            Vec3::new(-0.8, -1.0, 0.3),
+            TerrainChunkSize::RECT_SIZE.x,
+            TerrainChunkSize::RECT_SIZE.x as Alt,
+            TerrainChunkSize::RECT_SIZE.y as Alt,
+            Aabr {
+                min: Vec2::new(0.0, 0.0),
+                max: WORLD_SIZE.map(|e| e as Alt) * TerrainChunkSize::RECT_SIZE.map(|e| e as Alt),
+            },
+            0.0,
+            maxh,
+            |posi| alt[posi].max(water_alt[posi]),
+        ); */
+
         // Check whether any tiles around this tile are not water (since Lerp will
         // ensure that they are included).
         let pure_water = |posi: usize| {
@@ -1308,11 +1323,53 @@ impl WorldSim {
     /// Draw a map of the world based on chunk information.  Returns a buffer of
     /// u32s.
     pub fn get_map(&self) -> Vec<u32> {
+        let mut map_config = MapConfig::default();
+        map_config.lgain = 1.0;
+        // Build a shadow map.
+        let shadows = get_shadows(
+            Vec3::new(
+                map_config.light_direction.x,
+                map_config.light_direction.z,
+                map_config.light_direction.y,
+            ),
+            map_config.lgain,
+            TerrainChunkSize::RECT_SIZE.x as Alt,
+            TerrainChunkSize::RECT_SIZE.y as Alt,
+            Aabr {
+                min: Vec2::new(0.0, 0.0),
+                max: WORLD_SIZE.map(|e| e as Alt) * TerrainChunkSize::RECT_SIZE.map(|e| e as Alt),
+            },
+            CONFIG.sea_level as Alt,
+            (CONFIG.sea_level + self.max_height) as Alt,
+            |posi| {
+                let chunk = &self.chunks[posi];
+                chunk.alt.max(chunk.water_alt) as Alt
+            },
+        )
+        .unwrap();
+
+        let samples_data = {
+            let column_sample = ColumnGen::new(self);
+            (0..WORLD_SIZE.product())
+                .into_par_iter()
+                .map(|posi| {
+                    column_sample.get(
+                        uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32),
+                    )
+                })
+                .collect::<Vec<_>>()
+                .into_boxed_slice()
+        };
+
         let mut v = vec![0u32; WORLD_SIZE.x * WORLD_SIZE.y];
         // TODO: Parallelize again.
         MapConfig {
             gain: self.max_height,
-            ..MapConfig::default()
+            // lgain: 1.0,
+            shadows: Some(&shadows),
+            samples: Some(&samples_data),
+            // is_sampled: true,
+            ..map_config
         }
         .generate(&self, |pos, (r, g, b, a)| {
             v[pos.y * WORLD_SIZE.x + pos.x] = u32::from_le_bytes([r, g, b, a]);

world/src/sim/util.rs

@@ -1,4 +1,5 @@
 use super::WORLD_SIZE;
+use approx::ApproxEq;
 use bitvec::prelude::{bitbox, bitvec, BitBox};
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
 use noise::{MultiFractal, NoiseFn, Perlin, Point2, Point3, Point4, Seedable};
@@ -294,6 +295,52 @@ pub fn downhill<F: Float>(
         .into_boxed_slice()
 }
 
+/* /// Bilinear interpolation.
+///
+/// Linear interpolation in both directions (i.e. quadratic interpolation).
+fn get_interpolated_bilinear<T, F>(&self, pos: Vec2<i32>, mut f: F) -> Option<T>
+    where
+        T: Copy + Default + Signed + Float + Add<Output = T> + Mul<f32, Output = T>,
+        F: FnMut(Vec2<i32>) -> Option<T>,
+{
+    // (i) Find downhill for all four points.
+    // (ii) Compute distance from each downhill point and do linear interpolation on
+    // their heights. (iii) Compute distance between each neighboring point
+    // and do linear interpolation on       their distance-interpolated
+    // heights.
+
+    // See http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1990A%26A...239..443S&defaultprint=YES&page_ind=0&filetype=.pdf
+    //
+    // Note that these are only guaranteed monotone in one dimension; fortunately,
+    // that is sufficient for our purposes.
+    let pos = pos.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| {
+        e as f64 / sz as f64
+    });
+
+    // Orient the chunk in the direction of the most downhill point of the four.  If
+    // there is no "most downhill" point, then we don't care.
+    let x0 = pos.map2(Vec2::new(0, 0), |e, q| e.max(0.0) as i32 + q);
+    let y0 = f(x0)?;
+
+    let x1 = pos.map2(Vec2::new(1, 0), |e, q| e.max(0.0) as i32 + q);
+    let y1 = f(x1)?;
+
+    let x2 = pos.map2(Vec2::new(0, 1), |e, q| e.max(0.0) as i32 + q);
+    let y2 = f(x2)?;
+
+    let x3 = pos.map2(Vec2::new(1, 1), |e, q| e.max(0.0) as i32 + q);
+    let y3 = f(x3)?;
+
+    let z0 = y0
+        .mul(1.0 - pos.x.fract() as f32)
+        .mul(1.0 - pos.y.fract() as f32);
+    let z1 = y1.mul(pos.x.fract() as f32).mul(1.0 - pos.y.fract() as f32);
+    let z2 = y2.mul(1.0 - pos.x.fract() as f32).mul(pos.y.fract() as f32);
+    let z3 = y3.mul(pos.x.fract() as f32).mul(pos.y.fract() as f32);
+
+    Some(z0 + z1 + z2 + z3)
+} */
+
 /// Find all ocean tiles from a height map, using an inductive definition of
 /// ocean as one of:
 /// - posi is at the side of the world (map_edge_factor(posi) == 0.0)
@@ -335,6 +382,301 @@ pub fn get_oceans<F: Float>(oldh: impl Fn(usize) -> F + Sync) -> BitBox {
     is_ocean
 }
 
+/// Finds all shadowed chunks.
+///
+/// ray should be a nonzero vector (if it's not, we return an error).
+/// dx and dy should both be positive.
+pub fn get_shadows<F: std::iter::Sum + ApproxEq + Float + Send + Sync + std::fmt::Debug>(
+    ray: Vec3<F>,
+    lgain: F,
+    dx: F,
+    dy: F,
+    bounds: Aabr<F>,
+    minh: F,
+    maxh: F,
+    h: impl Fn(usize) -> F + Sync,
+) -> Result<Box<[F]>, ()> {
+    // First, make sure the ray and delta aren't zero.
+    let ray = -ray;
+    let ray_squared = ray.magnitude_squared();
+    if ray_squared == F::zero()
+        || ray_squared == F::neg_zero()
+        || !(dx > F::zero())
+        || !(dy > F::zero())
+    {
+        return Err(());
+    }
+    let hsize = Vec2::new(dx, dy);
+    /* if hstep.is_approx_zero() {
+        return Err(());
+    } */
+
+    // Find map sizes.
+    println!("Here?");
+    let wmap_size = Vec2::<F>::from(bounds.size()).map2(hsize, |e, f| e / f);
+    println!("Here?");
+    let map_size = if let Vec2 {
+        x: Some(x),
+        y: Some(y),
+    } = wmap_size.map(|e| F::to_usize(&e))
+    {
+        Vec2::new(x, y)
+    } else {
+        return Err(());
+    };
+    let map_len = map_size.product();
+    /* let distance_map = |wposf: Vec3<F>| {
+        // For all nodes in the height map, the minimum distance to a vertex is either
+        // the distance to the top of this chunk, or the distance to one of the neighbors.
+        let wpos = Vec2::new(i32::from(pos), i32::from(pos));
+        let posi = vec2_as_uniform_idx(wpos);
+        let neighbor_min = neighbors(posi)
+            .map(|posj| Vec3::new(vec2_as_uniform_idx(wpos), h(posj))
+            .min_by_key(|(pos, _)| pos.distance_squared(wpos))
+            .unwrap_or(F::infinity());
+        (wposf.z - ).min(neighobr_min)
+    }
+
+            .min_by
+        wposf.z.min(
+    }; */
+
+    // Make sure that the ray has a horizontal component at all; if not, the ray is
+    // vertical, so it can't cast a shadow, and we set the brightness for each
+    // chunk to 1.0.
+    if Vec2::<F>::from(ray).is_approx_zero() {
+        return Ok(vec![F::one(); map_len].into_boxed_slice());
+    }
+
+    // Conversely if the ray has no vertical component, each chunk must be entirely
+    // in shadow (at least for our purposes).
+    if ray.z == F::zero() || ray.z == F::neg_zero() {
+        return Ok(vec![F::zero(); map_len].into_boxed_slice());
+    }
+
+    // Otherwise, we can use step as the minimum length we need to raymarch in order
+    // to guarantee an accurate bounds check for the given dx and dy (i.e., all
+    // boundaries at spacings smaller than a "pixel" of size (dx,dy) should be
+    // taken into account).
+    let step_dot = if ray.x == F::zero() || ray.x == F::neg_zero() {
+        // Ray has no x component, so we must use y
+        dy * ray.y
+    } else if ray.y == F::zero() || ray.y == F::neg_zero() {
+        // Ray has no y component, so we must use x
+        dx * ray.x
+    } else {
+        // Ray has both an x and y component, so we must use the minimum.
+        if dy < dx { dy * ray.y } else { dx * ray.x }
+    }
+    .abs();
+    let step = ray * (step_dot / ray_squared);
+    let hstep = Vec2::from(step);
+    let zstep = step.z;
+
+    // Now, do the raymarching.
+    println!("Here?");
+    let max_steps = if let Some(max_steps) = ((maxh - minh) / zstep)
+        .abs()
+        .min(
+            wmap_size.reduce_partial_max(), /* / hstep.reduce_partial_min() */
+        )
+        .to_usize()
+    {
+        max_steps
+    } else {
+        return Err(());
+    };
+    println!("Here?");
+    let step_mag = step.magnitude();
+    let two = F::one() + F::one();
+    let three = two + F::one();
+    let w = F::from(0.5).unwrap();
+    let wstep_mag = lgain / (w * step_mag);
+    let wmax_steps = if let Some(wmax_steps) = F::from(max_steps) {
+        wmax_steps
+    } else {
+        return Err(());
+    };
+    println!("Here?");
+    let chunk_size = if let Vec2 {
+        x: Some(x),
+        y: Some(y),
+    } = TerrainChunkSize::RECT_SIZE.map(F::from)
+    {
+        Vec2::new(x, y)
+    } else {
+        return Err(());
+    };
+    println!(
+        "Here? map_len={:?}, map_size={:?}, hstep={:?}, zstep={:?} max_steps={:?}",
+        map_len, map_size, hstep, zstep, max_steps
+    );
+    Ok((0..map_len)
+        .into_par_iter()
+        .map(|posi| {
+            // Simple raymarch to determine whether we're in shadow.
+            // From https://www.iquilezles.org/www/articles/rmshadows/rmshadows.htm
+            // NOTE: How do we know these will succeed?
+            let wposf_orig = bounds.min
+                + Vec2::new(
+                    F::from(posi % map_size.x).unwrap(),
+                    F::from(posi / map_size.x).unwrap(),
+                ) * hsize;
+            let wpos_orig = wposf_orig.map2(chunk_size, |e, f| e / f);
+            // NOTE: How do we know these will succeed?
+            let wposl_orig = wpos_orig.map(|e| e.floor().to_i32().unwrap());
+            // NOTE: How do we know these will succeed?
+            // let wposr_orig = wpos_orig.map(|e| e.ceil().to_i32().unwrap());
+
+            let h_orig = if wposl_orig.reduce_partial_min() < 0 ||
+                // Casts are fine since we're a positive i32
+                /*wposr_orig*/wposl_orig.x as usize >= WORLD_SIZE.x ||
+                /*wposr_orig*/wposl_orig.y as usize >= WORLD_SIZE.y
+            {
+                // Out of bounds, assign minimum
+                minh
+            } else {
+                let hl_orig = h(vec2_as_uniform_idx(wposl_orig));
+                // let hr_orig = h(vec2_as_uniform_idx(wposr_orig));
+                // let wpos_frac = wposl_orig.map(|e| F::from(e).unwrap()).distance(wpos_orig);
+                hl_orig // hl_orig * wpos_frac + hr_orig * (F::one() - wpos_frac)
+            };
+
+            // let h_orig = h(vec2_as_uniform_idx(posi_orig.to_usize().unwrap()));
+            if h_orig < minh {
+                // Below the minimum height, always in shadow.
+                return F::zero();
+            }
+            let wmax_steps = wmax_steps.min(((maxh - h_orig) / zstep).abs());
+            // NOTE: How do we know these will succeed?
+            /* let wposf = bounds.min
+            + Vec2::new(
+                F::from(posi % map_size.x).unwrap(),
+                F::from(posi / map_size.x).unwrap(),
+            ) * hsize; */
+            let mut s = F::one();
+            // NOTE: How do we know these will succeed?
+            let mut wstep = F::zero();
+            let mut h_i = h_orig;
+            let mut wposf_i = wposf_orig;
+            while wstep < wmax_steps && s >= F::zero() {
+                wstep = wstep + F::one();
+                h_i = h_i + zstep;
+                wposf_i = wposf_i + hstep;
+                // Find height at this point.
+                // let h_i = h_orig + zstep * wstep;
+                // Find locations before and after h_i and use them to interpolate a height.
+                // let wposf_i = wposf + hstep * wstep;
+                let wpos_i = wposf_i.map2(chunk_size, |e, f| e / f);
+                // println!("h_orig={:?} h_i={:?}; wposf_orig={:?} wposf={:?}", h_orig, h_i,
+                // wposf, wposf_i);
+                // NOTE: How do we know these will succeed?
+                let wposl_i = wpos_i.map(|e| e.floor().to_i32().unwrap());
+                // NOTE: How do we know these will succeed?
+                // let wposr_i = wpos_i.map(|e| e.ceil().to_i32().unwrap());
+                if wposl_i.reduce_partial_min() < 0 ||
+                    // Casts are fine since we're a positive i32
+                    /*wposr_i*/wposl_i.x as usize >= WORLD_SIZE.x ||
+                    /*wposr_i*/wposl_i.y as usize >= WORLD_SIZE.y
+                {
+                    // Out of bounds, we're done!
+                    break;
+                }
+                let hl_i = h(vec2_as_uniform_idx(wposl_i));
+                // let hr_i = h(vec2_as_uniform_idx(wposr_i));
+                // let wpos_frac = wposl_i.map(|e| F::from(e).unwrap()).distance(wpos_i);
+                let h_j = hl_i; //hl_i * wpos_frac + hr_i * (F::one() - wpos_frac);
+                //
+                /* let posj = wpos_i.map(|e| e.into_i32());
+                let h_j = h(posj); */
+                // If we landed in shadow, we're done.
+                // println!("h_orig={:?} h_i={:?} h_j={:?}; wposf_orig={:?} wposf={:?}
+                // wpos_frac={:?}", h_orig, h_i, h_j, wposf, wposf_i, wpos_frac);
+                // NOTE: Approximation to real distance metric, which is hard to compute for an
+                // arbitrary point in a heightmap.
+                // s = s.min((h_i - h_j) * lgain / (wstep * w * step_mag));
+                s = s.min((h_i - h_j) * wstep_mag / wstep);
+                /* if s < F::zero()
+                /* h_i < h_j */
+                {
+                    // s = F::zero();
+                    // println!("A shadow!");
+                    break;
+                } */
+            }
+            /*while h_i + zstep * posj < maxh {
+                wposf += hstep;
+                let posj : Vec2<i32> = (wposf / chunk_size).into();
+                // Simple interpolation.
+                let h_j = h(posj);
+                let h_j = h_i + (h_j - h_i) * (hstep / ).magnitude();
+                h_i += zstep;
+                // If we landed in shadow, we're done.
+                if h_i > h_j {
+                    return 0.0;
+                }
+            } */
+            s = s.max(F::zero());
+            // Smoothstep
+            s * s * (three - two * s)
+            // // Above the maximum height, definitely not in shadow.
+            // return F::one();
+            /* let nh = h(posi);
+            if is_ocean(posi) {
+                -2
+            } else {
+                let mut best = -1;
+                let mut besth = nh;
+                for nposi in neighbors(posi) {
+                    let nbh = h(nposi);
+                    if nbh < besth {
+                        besth = nbh;
+                        best = nposi as isize;
+                    }
+                }
+                best
+            } */
+        })
+        .collect::<Vec<_>>()
+        .into_boxed_slice())
+
+    /* // Raymarching technique to quickly identify approxmate distances
+    // Use a shadow map to raymarch all the height entries.
+    // We can mark tiles as ocean candidates by scanning row by row, since the top
+    // edge is ocean, the sides are connected to it, and any subsequent ocean
+    // tiles must be connected to it.
+    let mut is_ocean = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
+    let mut stack = Vec::new();
+    let mut do_push = |pos| {
+        let posi = vec2_as_uniform_idx(pos);
+        if oldh(posi) <= F::zero() {
+            stack.push(posi);
+        }
+    };
+    for x in 0..WORLD_SIZE.x as i32 {
+        do_push(Vec2::new(x, 0));
+        do_push(Vec2::new(x, WORLD_SIZE.y as i32 - 1));
+    }
+    for y in 1..WORLD_SIZE.y as i32 - 1 {
+        do_push(Vec2::new(0, y));
+        do_push(Vec2::new(WORLD_SIZE.x as i32 - 1, y));
+    }
+    while let Some(chunk_idx) = stack.pop() {
+        // println!("Ocean chunk {:?}: {:?}", uniform_idx_as_vec2(chunk_idx),
+        // oldh(chunk_idx));
+        if *is_ocean.at(chunk_idx) {
+            continue;
+        }
+        *is_ocean.at(chunk_idx) = true;
+        stack.extend(neighbors(chunk_idx).filter(|&neighbor_idx| {
+            // println!("Ocean neighbor: {:?}: {:?}", uniform_idx_as_vec2(neighbor_idx),
+            // oldh(neighbor_idx));
+            oldh(neighbor_idx) <= F::zero()
+        }));
+    }
+    is_ocean */
+}
+
 /// A 2-dimensional vector, for internal use.
 type Vector2<T> = [T; 2];
 /// A 3-dimensional vector, for internal use.