Use bilinear patches in height map rendering

Based on "Maximum Mipmaps for Fast, Accurate, and Scalable Dynamic
Height Field Rendering"

README.md

@@ -4,7 +4,7 @@ Peaks is a (work in progress) tool for rendering 3d maps.
 
 ## TODO
 
-* [ ] Ray bilinear patch intersection for terrain
+* [x] Ray bilinear patch intersection for terrain
 * [ ] Better sub sampling method (Poisson disc...)
 * [ ] Orthographic camera
 * [ ] Use simulated film values for setting field of view

examples/ben_nevis/main.rs

@@ -82,15 +82,10 @@ pub fn main() -> Result<()> {
         Texture::blank(raw_height_data.width, raw_height_data.height);
     peaks::ops::scale(&raw_height_data, &mut height_map, vertical_exageration);
 
-    let w = height_map.width;
-    let h = height_map.height;
-
-    let mut normal_map = Texture::blank(w, h);
+    let mut normal_map = Texture::blank(height_map.width, height_map.height);
     peaks::ops::normals(&height_map, &mut normal_map);
-    let quad_tree = peaks::spatial::height_map_quad_tree(&height_map);
 
-    let height_map =
-        HeightMap::new(transform, height_map, normal_map, quad_tree);
+    let height_map = HeightMap::new(transform, height_map, normal_map);
 
     let scene = Scene {
         background: Vec3::new(254.0, 254.0, 200.0) / 255.0,

examples/height_map/main.rs

@@ -37,8 +37,6 @@ pub fn main() -> Result<()> {
     peaks::ops::rgb_height_map(&mono_map, &mut height_map, 20.0);
     peaks::ops::normals(&height_map, &mut normal_map);
 
-    let quad_tree = peaks::spatial::height_map_quad_tree(&height_map);
-
     let height_map = HeightMap::new(
         [
             -(mono_map.width as f64 / 2.0),
@@ -48,7 +46,6 @@ pub fn main() -> Result<()> {
         ],
         height_map,
         normal_map,
-        quad_tree,
     );
 
     let width = 960;

examples/mljet/main.rs

@@ -73,9 +73,7 @@ pub fn main() -> Result<()> {
     let mut normal_map = Texture::blank(w, h);
     peaks::ops::normals(&height_map, &mut normal_map);
 
-    let quad_tree = peaks::spatial::height_map_quad_tree(&height_map);
-    let height_map =
-        HeightMap::new(transform, height_map, normal_map, quad_tree);
+    let height_map = HeightMap::new(transform, height_map, normal_map);
 
     let scene = Scene {
         background: Vec3::new(254.0, 254.0, 200.0) / 255.0,

src/lib.rs

@@ -25,7 +25,6 @@ pub mod ops;
 mod primitives;
 mod render;
 mod samplers;
-pub mod spatial;
 mod textures;
 
 pub use cameras::PinholeCamera;
@@ -36,5 +35,4 @@ pub use render::{
     Object, Renderer, Scene,
 };
 pub use samplers::RegularGridSampler;
-pub use spatial::{QuadTree, QuadTreeNode};
 pub use textures::Texture;

src/math.rs

@@ -19,6 +19,13 @@ use std::ops::{
     Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign,
 };
 
+/// Round to the nearest power of two
+pub fn ceil_pow2(num: usize) -> usize {
+    let num = num as f64;
+    let exp = (num.log2() / 2.0_f64.log2()).ceil();
+    2.0_f64.powf(exp) as usize
+}
+
 #[derive(Debug, Default, Copy, Clone, PartialEq)]
 pub struct Color {
     pub r: u8,

src/ops.rs

@@ -82,6 +82,43 @@ where
     }
 }
 
+/// Create map of bilinear patches and its first mipmap level from a height map
+pub fn height_map_to_bilinear_patch(
+    input: &Texture<f64>,
+    level0: &mut Texture<[f64; 4]>,
+    level1: &mut Texture<f64>,
+) {
+    assert_eq!(level0.width, level1.width);
+    assert_eq!(level0.height, level1.height);
+    assert_eq!(input.width - 1, level0.width);
+    assert_eq!(input.height - 1, level0.height);
+
+    for y in 0..level0.width {
+        for x in 0..level0.height {
+            // Read data in `z` order but write out in a clockwise order
+            let [nw, ne, sw, se] = input.lookup2x2(x, y);
+            level0.write1x1(x, y, [nw, ne, se, sw]);
+            level1.write1x1(x, y, nw.max(ne).max(se).max(sw));
+        }
+    }
+}
+
+/// Create the next maximum mipmap level for a floating point texture
+pub fn maximum_mipmap_bilinear_patch(
+    input: &Texture<f64>,
+    output: &mut Texture<f64>,
+) {
+    assert_eq!(input.width / 2, output.width);
+    assert_eq!(input.height / 2, output.height);
+
+    for y in 0..output.height {
+        for x in 0..output.width {
+            let [p1, p2, p3, p4] = input.lookup2x2(x * 2, y * 2);
+            output.write1x1(x, y, p1.max(p2).max(p3).max(p4));
+        }
+    }
+}
+
 /// Scale a surface by `n`
 pub fn scale<T>(input: &Texture<T>, output: &mut Texture<T>, n: f64)
 where
@@ -207,4 +244,70 @@ mod tests {
         assert_eq!(dest.lookup1x1(6, 6), 250.0);
         assert_eq!(dest.lookup1x1(7, 7), 255.0);
     }
+
+    #[test]
+    fn test_maximum_mipmaps_bilinear_patches() {
+        #[cfg_attr(rustfmt, rustfmt_skip)]
+        let mut bilinear_patches = Texture::blank(4, 4);
+        let mut bilinear_patches_mipmap0 = Texture::blank(4, 4);
+        let mut bilinear_patches_mipmap1 = Texture::blank(2, 2);
+        let mut bilinear_patches_mipmap2 = Texture::blank(1, 1);
+        let height_map = Texture::new(
+            5,
+            5,
+            vec![
+                1.0, 3.5, 6.0, 8.5, 11.0, 1.5, 4.0, 6.5, 9.0, 11.5, 2.0, 4.5,
+                7.0, 9.5, 12.0, 2.5, 5.0, 7.5, 10.0, 12.5, 3.0, 5.5, 8.0, 10.5,
+                13.0,
+            ],
+        );
+
+        height_map_to_bilinear_patch(
+            &height_map,
+            &mut bilinear_patches,
+            &mut bilinear_patches_mipmap0,
+        );
+        maximum_mipmap_bilinear_patch(
+            &bilinear_patches_mipmap0,
+            &mut bilinear_patches_mipmap1,
+        );
+        maximum_mipmap_bilinear_patch(
+            &bilinear_patches_mipmap1,
+            &mut bilinear_patches_mipmap2,
+        );
+
+        #[cfg_attr(rustfmt, rustfmt_skip)]
+        assert_eq!(bilinear_patches.buffer, [
+            // Row 1
+            [1.0, 3.5, 4.0, 1.5],
+            [3.5, 6.0, 6.5, 4.0],
+            [6.0, 8.5, 9.0, 6.5],
+            [8.5, 11.0, 11.5, 9.0],
+            // Row 2
+            [1.5, 4.0, 4.5, 2.0],
+            [4.0, 6.5, 7.0, 4.5],
+            [6.5, 9.0, 9.5, 7.0],
+            [9.0, 11.5, 12.0, 9.5],
+            // Row 3
+            [2.0, 4.5, 5.0, 2.5],
+            [4.5, 7.0, 7.5, 5.0],
+            [7.0, 9.5, 10.0, 7.5],
+            [9.5, 12.0, 12.5, 10.0],
+            // Row 4
+            [2.5, 5.0, 5.5, 3.0],
+            [5.0, 7.5, 8.0, 5.5],
+            [7.5, 10.0, 10.5, 8.0],
+            [10.0, 12.5, 13.0, 10.5],
+        ]);
+
+        #[cfg_attr(rustfmt, rustfmt_skip)]
+        assert_eq!(bilinear_patches_mipmap0.buffer, [
+            4.0, 6.5,  9.0, 11.5,
+            4.5, 7.0,  9.5, 12.0,
+            5.0, 7.5, 10.0, 12.5,
+            5.5, 8.0, 10.5, 13.0,
+        ]);
+        assert_eq!(bilinear_patches_mipmap1.buffer, [7.0, 12.0, 8.0, 13.0]);
+        assert_eq!(bilinear_patches_mipmap2.buffer, [13.0]);
+    }
 }