README.md

@@ -2,8 +2,6 @@
 
 Make beautiful maps.
 
-Very rough proof-of-concept for real-time raytracing or digital terrain models.
-
 ![Image of Spectacled Flying Fox](Spectacled_flying_fox.jpg)
 
 Spectacled Flying Fox image By CSIRO, [Creative Commons Attribution

src/app/dem_renderer/dem_renderer.wgsl

@@ -33,78 +33,60 @@ fn vs_main(
     return result;
 }
 
-const RGSS_WIGGLE = array(vec2(0.125, 0.375),
-                          vec2(-0.125, -0.375),
-                          vec2(0.375, 0.125),
-                          vec2(-0.375, -0.125));
-
 @fragment
 fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
     var ray: Ray;
     ray.origin = vertex.world_space_position.xyz;
     ray.direction = normalize(vertex.world_space_position.xyz - uniforms.camera_position.xyz);
-    // Spread rays into RGSS antialiasing pattern.
-    let ray_dx = dpdy(ray.direction);
-    let ray_dy = dpdy(ray.direction);
-    //var ray_directions: array<vec3<f32>,4>;
-    //ray_directions[0] = ray.direction + 0.125*ray_dx + 0.375*ray_dy;
-    //ray_directions[1] = ray.direction - 0.125*ray_dx - 0.375*ray_dy;
-    //ray_directions[2] = ray.direction + 0.375*ray_dx + 0.125*ray_dy;
-    //ray_directions[3] = ray.direction - 0.375*ray_dx - 0.125*ray_dy;
 
-    // Possibly these ray directions could be stored in a matrix and we could
-    // evaluate them all at once instead of looping.
+    var root_node: BoundingNode;
+    root_node.index = vec2<u32>(0);
+    root_node.level = textureNumLevels(dembvh_texture) - 1;
 
-    let sun_direction = vec3<f32>(0.761904762, 0.380952381, 0.19047619);
-    
-    var color_accumulator = vec4<f32>(0);
-    let rgss_value = vec2(0.125, 0.375);
-    var wiggled_ray = ray;
-    for(var i=0; i<4; i++) {
-      let rgss_wiggle = select(rgss_value, rgss_value.yx, vec2(bool(i&2))) * select(1.0, -1.0, bool(i&1));
-      wiggled_ray.direction = ray.direction + rgss_wiggle.x * ray_dx + rgss_wiggle.y * ray_dy;
-
-       var root_node: BoundingNode;
-       root_node.index = vec2<u32>(0);
-       root_node.level = textureNumLevels(dembvh_texture) - 1;
-
-       var hit_location: vec3<f32>;
-       var hit_normal: vec3<f32>;
-       if intersect_ray_with_node(wiggled_ray, root_node, &hit_location, &hit_normal) {
-           var shadow_value = 0.0f;
-           var shadow_ray :Ray;
-           shadow_ray.origin = hit_location + hit_normal * 0.1;
-           for(var j=0; i<8; i++) {
-             let shadow_sample = 0.02 * (vec3(1.0, 0.0, 0.0) * f32(bool((i*4+j)&1))
-                                         + vec3(0.0, 1.0, 0.0) * f32(bool((i*4+j)&2))
-                                         + vec3(-1.0, 0.0, 0.0) * f32(bool((i*4+j)&4))
-                                         + vec3(0.0, -1.0, 0.0) * f32(bool((i*4+j)&8))
-                                         + vec3(0.0, 0.0, 1.0) * f32(bool((i*4+j)&16)));
-               // Calculate light
-               let sun_direction =
-                 (uniforms.camera_to_world_matrix
-                  * vec4<f32>(sun_direction + shadow_sample, 0.0)).xyz;
-               shadow_ray.direction = sun_direction;
-               var dummy0: vec3<f32>;
-               var dummy1: vec3<f32>;
-               shadow_value +=
-                 select(0.25, 0.0, intersect_ray_with_node(shadow_ray,
-                                                              root_node,
-                                                              &dummy0,
-                                                              &dummy1));
-           }
-           let sun_direction =
-             (uniforms.camera_to_world_matrix
-              * vec4<f32>(sun_direction, 0.0)).xyz;
-           let sun_lambertian_value = dot(hit_normal, sun_direction);
-           let fill_lambertian_value = dot(hit_normal, sun_direction * vec3(1.0, -1.0, 1.0));
-           let fill_strength = 0.25;
-           let l = fill_strength * fill_lambertian_value + (1.0 - fill_strength) * sun_lambertian_value * shadow_value;
-           color_accumulator += 0.25 * vec4<f32>(vec3<f32>(l), 1.0);
-        }
-    }
-    if color_accumulator.a == 0.0 {
+    var hit_index = vec2<u32>(0);
+    if intersect_ray_with_node(dembvh_texture,
+                               uniforms.dem_min_corner,
+                               uniforms.dem_cell_size,
+                               uniforms.dem_z_range,
+                               ray,
+                               root_node,
+                               &hit_index) {
+        //Calculate x-values of cell corners
+        let v00 = textureLoad(dem_texture, hit_index, 0).r;
+        let v01 = textureLoad(dem_texture, hit_index + vec2<u32>(0,1), 0).r;
+        let v10 = textureLoad(dem_texture, hit_index + vec2<u32>(1,0), 0).r;
+        let v11 = textureLoad(dem_texture, hit_index + vec2<u32>(1,1), 0).r;
+        let z00 = 12.0 * f32(v00) / 65535.0;
+        let z01 = 12.0 * f32(v01) / 65535.0;
+        let z10 = 12.0 * f32(v10) / 65535.0;
+        let z11 = 12.0 * f32(v11) / 65535.0;
+        // Calculate xyz of cell corners
+        let xy00 = uniforms.dem_min_corner.xy
+          + uniforms.dem_cell_size * vec2<f32>(hit_index);
+        let p00 = vec3<f32>(xy00, z00);
+        let p01 = vec3<f32>(xy00 + vec2<f32>(0, 1), z01);
+        let p10 = vec3<f32>(xy00 + vec2<f32>(1, 0), z10);
+        let p11 = vec3<f32>(xy00 + vec2<f32>(1, 1), z11);
+        // Intersect ray with the plane of each triangle and then take the
+        // point that's inside it's triangle.
+        // p is the point and p0, p1 and p2 are the triangle corners.
+        let p_t0 = intersect_ray_with_triangle_plane(ray, p00, p01, p11);
+        let p_t1 = intersect_ray_with_triangle_plane(ray, p00, p10, p11);
+        let point_is_in_first_triangle = point_is_in_triangle(p_t0, p00, p01, p11);
+        let p = select(p_t1, p_t0, point_is_in_first_triangle);
+        let p0 = p00;
+        let p1 = select(p10, p01, point_is_in_first_triangle);
+        let p2 = p11;
+        // Estimate normal as cross product of triangle for now
+        var normal = normalize(cross(p1-p0, p2-p0));
+        normal *= sign(normal.z); //Ensure normal points up
+        // Calculate light
+        let sun_direction = (uniforms.camera_to_world_matrix
+                             * vec4<f32>(normalize(vec3<f32>(1.0, 1.0, 1.0)), 0.0)).xyz;
+        let l = dot(normal, sun_direction);
+        return vec4<f32>(l, l, l, 1.0);
+        //return vec4<f32>(normal.xy, 0.0, 1.0);
+    } else {
         discard;
     }
-    return color_accumulator;
 }

src/app/dem_renderer/mod.rs

@@ -494,7 +494,7 @@ fn create_dembvh_texture(
 }
 
 fn get_animated_camera_position(animation_time: std::time::Duration, dem_size: f32) -> glam::Vec3 {
-    let animation_phase = 2.0 * std::f32::consts::PI * (animation_time.as_secs_f32() % 100.0) / 100.0;
+    let animation_phase = 2.0 * std::f32::consts::PI * (animation_time.as_secs_f32() % 25.0) / 25.0;
     glam::Vec3::new(
         dem_size * f32::sin(animation_phase),
         dem_size * f32::cos(animation_phase),

src/app/dem_renderer/ray_intersection.wgsl

@@ -13,7 +13,7 @@ struct AABB {
 struct Ray { origin: vec3<f32>, direction: vec3<f32> };
 
 fn invert_ray_direction(v: vec3<f32> ) -> vec3<f32> {
-  return select(vec3<f32>(1.0e+30),
+    return select(vec3<f32>(1.0e+30),
                   vec3<f32>(1.0) / v,
                   vec3<bool>(v));
 }
@@ -31,44 +31,43 @@ fn intersect_ray_with_aabb_optimized(ray_origin: vec3<f32>, inv_ray_direction: v
     let t_min = max(t_mins.x, max(t_mins.y, t_mins.z));
     let t_maxs = max(t1, t2);
     let t_max = min(t_maxs.x, min(t_maxs.y, t_maxs.z));
-    return select(-1.0, max(t_min, 0.0f),  t_min <= t_max && t_max > 0.0);
+    return select(-1.0, max(t_min, 0.0f),  t_min <= t_max);
 }
 
-fn get_xy_min_for_node(node: BoundingNode) -> vec2<f32> {
-    return uniforms.dem_min_corner.xy + uniforms.dem_cell_size * vec2<f32>(node.index) * exp2(f32(node.level));
+fn get_xy_min_for_node(dem_min_corner: vec2<f32>,
+    dem_cell_size: vec2<f32>,
+    node: BoundingNode) -> vec2<f32> {
+    return dem_min_corner.xy + dem_cell_size * vec2<f32>(node.index) * exp2(f32(node.level));
 }
 
-fn get_xy_max_for_node(node: BoundingNode) -> vec2<f32> {
-    return uniforms.dem_min_corner.xy + uniforms.dem_cell_size * vec2<f32>(node.index + 1) * exp2(f32(node.level));
+fn get_xy_max_for_node(dem_min_corner: vec2<f32>,
+    dem_cell_size: vec2<f32>,
+    node: BoundingNode) -> vec2<f32> {
+    return dem_min_corner.xy + dem_cell_size * vec2<f32>(node.index + 1) * exp2(f32(node.level));
 }
 
 struct NodeStack {
-  stack: array<u32,64>,
+    stack: array<BoundingNode,64>,
     count: u32
 }
 
 fn push_node_stack(stack: ptr<function,NodeStack>, node: BoundingNode) {
-    var packed_node = (node.index.x & 0x1fff) << 19;
-    packed_node += (node.index.y & 0x1fff) << 6;
-    packed_node += node.level & 0x3f;
-    (*stack).stack[(*stack).count] = packed_node;
+    (*stack).stack[(*stack).count] = node;
     (*stack).count++;
 }
 
 fn pop_node_stack(stack: ptr<function,NodeStack>) -> BoundingNode {
     (*stack).count--;
-    let packed_node = (*stack).stack[(*stack).count];
-    var node: BoundingNode;
-    node.index.x = (packed_node >> 19) & 0x1fff;
-    node.index.y = (packed_node >> 6) & 0x1fff;
-    node.level = packed_node & 0x3f;
-    return node;
+    return (*stack).stack[(*stack).count];
 }
 
-fn intersect_ray_with_node(ray: Ray,
+fn intersect_ray_with_node(tree_texture: texture_2d<u32>,
+                           dem_min_corner: vec2<f32>,
+                           dem_cell_size: vec2<f32>,
+                           dem_z_range: vec2<f32>,
+                           ray: Ray,
                            root_node: BoundingNode,
-                           hit_location: ptr<function, vec3<f32>>,
-                           hit_normal: ptr<function, vec3<f32>>) -> bool {
+                           hit_cell: ptr<function, vec2<u32>>) -> bool {
     let inv_ray_direction = invert_ray_direction(ray.direction);
     var node_stack: NodeStack;
     node_stack.count = 0u;
@@ -76,29 +75,27 @@ fn intersect_ray_with_node(ray: Ray,
     var closest_hit_distance = 1.0e+30f;
     while node_stack.count > 0 {
         let node = pop_node_stack(&node_stack);
-        let minmax_z = textureLoad(dembvh_texture, node.index, i32(node.level)).rg;
+        let minmax_z = textureLoad(tree_texture, node.index, i32(node.level)).rg;
         if minmax_z.r == 0 {
             return false;
         }
-        let min_z = scale_u16(minmax_z.r, uniforms.dem_z_range);
-        let max_z = scale_u16(minmax_z.g, uniforms.dem_z_range);
+        let min_z = scale_u16(minmax_z.r, dem_z_range);
+        let max_z = scale_u16(minmax_z.g, dem_z_range);
         var aabb: AABB ;
-        aabb.min_corner = vec3<f32>(get_xy_min_for_node(node), min_z);
-        aabb.max_corner = vec3<f32>(get_xy_max_for_node(node), max_z);
+        aabb.min_corner = vec3<f32>(get_xy_min_for_node(dem_min_corner,
+            dem_cell_size,
+            node),
+            min_z);
+        aabb.max_corner = vec3<f32>(get_xy_max_for_node(dem_min_corner,
+            dem_cell_size,
+            node),
+            max_z);
         let hit_distance = intersect_ray_with_aabb_optimized(ray.origin, inv_ray_direction, aabb);
         if hit_distance >= 0.0 {
             if node.level == 0 {
-                var location: vec3<f32>;
-                var normal: vec3<f32>;
                 if hit_distance < closest_hit_distance {
-                    if intersect_ray_with_grid_cell(ray, node.index, &location, &normal)
-                    {
-                        // Node bounding boxes are non-overlapping, so we don't need
-                        // to calculate a more precise hit_distance
-                        closest_hit_distance = hit_distance;
-                        *hit_location = location;
-                        *hit_normal = normal;
-                    }
+                    closest_hit_distance = hit_distance;
+                    *hit_cell = node.index;
                 }
             } else {
                 let next_index = node.index * 2;
@@ -118,66 +115,19 @@ fn intersect_ray_with_node(ray: Ray,
     return closest_hit_distance < 1.0e+20;
 }
 
-fn intersect_ray_with_grid_cell(ray: Ray,
-                                cell_index: vec2<u32>,
-                                location: ptr<function, vec3<f32>>,
-                                normal: ptr<function, vec3<f32>>) -> bool {
-    //Get z-values of cell corners
-    // TODO: Handle missing data
-    let v00 = textureLoad(dem_texture, cell_index, 0).r;
-    let v01 = textureLoad(dem_texture, cell_index + vec2<u32>(0,1), 0).r;
-    let v10 = textureLoad(dem_texture, cell_index + vec2<u32>(1,0), 0).r;
-    let v11 = textureLoad(dem_texture, cell_index + vec2<u32>(1,1), 0).r;
-    let z00 = scale_u16(v00, uniforms.dem_z_range);
-    let z01 = scale_u16(v01, uniforms.dem_z_range);
-    let z10 = scale_u16(v10, uniforms.dem_z_range);
-    let z11 = scale_u16(v11, uniforms.dem_z_range);
-    // Calculate xyz of cell corners
-    let xy00 = uniforms.dem_min_corner.xy + uniforms.dem_cell_size * vec2<f32>(cell_index);
-    let p00 = vec3<f32>(xy00, z00);
-    let p01 = vec3<f32>(xy00 + vec2<f32>(0, 1) * uniforms.dem_cell_size, z01);
-    let p10 = vec3<f32>(xy00 + vec2<f32>(1, 0) * uniforms.dem_cell_size, z10);
-    let p11 = vec3<f32>(xy00 + vec2<f32>(1, 1) * uniforms.dem_cell_size, z11);
-
-    // Intersect ray with the plane of each triangle and then take the
-    // point that's inside it's triangle.
-    // p is the point and p0, p1 and p2 are the triangle corners.
-    var hit_t0: bool;
-    let p_t0 = intersect_ray_with_triangle_plane(ray, p00, p01, p11, &hit_t0);
-    var hit_t1: bool;
-    let p_t1 = intersect_ray_with_triangle_plane(ray, p00, p10, p11, &hit_t1);
-    let point_is_in_triangle_0 = hit_t0 && point_is_in_triangle(p_t0, p00, p01, p11);
-    let point_is_in_triangle_1 = hit_t1 && point_is_in_triangle(p_t1, p00, p10, p11);
-
-    // Always calculate location and normal, even if ray doesn't intersect
-    // either triangle. It might be better to do an early return if there's
-    // no intersection but it's probably not worth the thread divergence.
-    *location = select(p_t1, p_t0, point_is_in_triangle_0);
-    let p0 = p00;
-    let p1 = select(p10, p01, point_is_in_triangle_0);
-    let p2 = p11;
-    // Estimate normal as cross product of triangle for now
-    *normal = normalize(cross(p1-p0, p2-p0));
-    *normal *= sign((*normal).z); //Ensure normal points up
-
-    return point_is_in_triangle_0 || point_is_in_triangle_1;
-}
 
 fn intersect_ray_with_triangle_plane(ray: Ray,
                                      p0: vec3<f32>,
                                      p1: vec3<f32>,
-                                     p2: vec3<f32>,
-                                     hit: ptr<function,bool>) -> vec3<f32> {
+                                     p2: vec3<f32>) -> vec3<f32> {
     let n = cross(p1-p0, p2-p0);
-    let t = dot(p0 - ray.origin, n) / dot(ray.direction, n);
-    *hit = t >= 0.0;
-    return ray.origin + ray.direction * t;
+    return ray.origin + ray.direction * dot(p0 - ray.origin, n) / dot(ray.direction, n);
 }
 
 
 fn point_is_inside_triangle_side(p: vec3<f32>,
                                  p0: vec3<f32>, p1: vec3<f32>, p2: vec3<f32>) -> bool {
-  return dot(cross(p2-p1, p-p1), cross(p2-p1, p0-p1)) >= -0.0;
+  return dot(cross(p2-p1, p-p1), cross(p2-p1, p0-p1)) >= 0;
 }
 
 

src/app/dem_renderer/tests.rs

@@ -140,14 +140,14 @@ fn test_shaders() {
             aabb_max_corner: vec3(1.0, 1.0, 1.0),
         },
     ];
-    let output_buffer: Vec<f32> = bytemuck::cast_slice(&block_on(run_compute_shader_test(
+    let output_buffer: Vec<u32> = bytemuck::cast_slice(&block_on(run_compute_shader_test(
         wgsl_module!("tests.wgsl"),
         "test_intersect_ray_with_aabb",
         bytemuck::cast_slice(&input_buffer),
         input_buffer.len() * 4,
     )))
-        .to_vec();
-    assert!(dbg!(output_buffer[0]) >= 0.0);
-    assert!(dbg!(output_buffer[1]) < 0.0);
-    assert!(dbg!(output_buffer[2]) >= 0.0);
+    .to_vec();
+    assert_eq!(output_buffer[0], 1);
+    assert_eq!(output_buffer[1], 0);
+    assert_eq!(output_buffer[2], 1);
 }

src/app/dem_renderer/tests.wgsl

@@ -3,7 +3,7 @@
 var<storage, read> input_data: array<Input>;
 @group(0)
 @binding(1)
-var<storage, read_write> output_data: array<f32>;
+var<storage, read_write> output_data: array<u32>;
 
 //#include ray_intersection.wgsl
 
@@ -24,11 +24,10 @@ fn test_intersect_ray_with_aabb() {
       var aabb: AABB;
       aabb.min_corner = input_data[i].aabb_min_corner;
       aabb.max_corner = input_data[i].aabb_max_corner;
-      output_data[i] = intersect_ray_with_aabb(ray, aabb);
-      /*if intersect_ray_with_aabb(ray, aabb) >= 0.0 {
+      if intersect_ray_with_aabb(ray, aabb) >= 0.0 {
           output_data[i] = 1u;
       }else {
           output_data[i] = 0u;
-          }*/
+      }
     }
 }

src/app/raster/dem.rs

@@ -174,20 +174,18 @@ impl Dem {
         let (num_cells_x, num_cells_y) = tiff.dimensions().unwrap();
         match tiff.read_image().unwrap() {
             tiff::decoder::DecodingResult::F32(f32_values) => {
-                let x_min = -500.0;
-                let x_max = 500.0;
-                let y_min = -500.0;
-                let y_max = 500.0;
+                let x_min = -5000.0;
+                let x_max = 5000.0;
+                let y_min = -5000.0;
+                let y_max = 5000.0;
                 let (z_min, z_max) = f32_values[1..]
                     .iter()
-                    .map(|z| z * 0.1)
                     .fold((f32_values[0], f32_values[0]), |(min, max), elem| {
-                        (min.min(elem), max.max(elem))
+                        (min.min(*elem), max.max(*elem))
                     });
                 let grid = f32_values
                     .iter()
-                    .map(|z| z * 0.1)
-                    .map(|v| normalized_f32_to_u16(normalize_f32(v, z_min, z_max)))
+                    .map(|v| normalized_f32_to_u16(normalize_f32(*v, z_min, z_max)))
                     .collect();
                 Dem {
                     num_cells_x,