It's a mess but it works.

src/app/dem_renderer/dem_renderer.wgsl

@@ -51,8 +51,8 @@ fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
     // evaluate them all at once instead of looping.
 
     var color_accumulator = vec4<f32>(0);
-    for(var i=0; i<4; i++) {
+    for(var i=0; i<1; i++) {
-       ray.direction = ray_directions[i];
+      //ray.direction = ray_directions[i];
 
        var root_node: BoundingNode;
        root_node.index = vec2<u32>(0);
@@ -61,7 +61,8 @@ fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
        var hit_index = vec2<u32>(0);
        var hit_location: vec3<f32>;
        var hit_normal: vec3<f32>;
-       if intersect_ray_with_node(dembvh_texture,
+       if intersect_ray_with_node(dem_texture,
+                                  dembvh_texture,
                                   uniforms.dem_min_corner,
                                   uniforms.dem_cell_size,
                                   uniforms.dem_z_range,
@@ -70,31 +71,48 @@ fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
                                   &hit_index,
                                   &hit_location,
                                   &hit_normal) {
+           let sun_direction_samples = array<vec3<f32>,7>(vec3<f32>(1.0, 0.5, 0.25),
+                                                          vec3<f32>(1.0, 0.33, 0.25),
+                                                          vec3<f32>(1.1, 0.3, 0.25),
+                                                          vec3<f32>(1.1, 0.33, 0.25),
+                                                          vec3<f32>(1.0, 0.27, 0.25),
+                                                          vec3<f32>(0.9, 0.3, 0.25),
+                                                          vec3<f32>(0.9, 0.27, 0.25));
+           var shadow_value = 0.0f;
+           var shadow_ray :Ray;
+           shadow_ray.origin = hit_location + hit_normal * 0.1;
+           for(var i=0; i<1; i++) {
                // Calculate light
                let sun_direction =
                  (uniforms.camera_to_world_matrix
-              * vec4<f32>(normalize(vec3<f32>(1.0, 1.0, 1.0)), 0.0)).xyz;
+                  * vec4<f32>(normalize(sun_direction_samples[i]), 0.0)).xyz;
-           let l = dot(hit_normal, sun_direction);
-           var shadow_ray :Ray;
-           shadow_ray.origin = hit_location + hit_normal * 0.1;
                shadow_ray.direction = sun_direction;
-           var dummy0: vec3<f32>;
+               var dummy0: vec2<u32>;
                var dummy1: vec3<f32>;
-           let shadow_value =
+               var dummy2: vec3<f32>;
-             select(1.0, 0.0, intersect_ray_with_node(dembvh_texture,
+               shadow_value +=
+                 select(1.0, 0.0, intersect_ray_with_node(dem_texture,
+                                                          dembvh_texture,
                                                           uniforms.dem_min_corner,
                                                           uniforms.dem_cell_size,
                                                           uniforms.dem_z_range,
                                                           shadow_ray,
                                                           root_node,
-                                                      &hit_index,
                                                           &dummy0,
-                                                      &dummy1));
+                                                          &dummy1,
-           color_accumulator +=  vec4<f32>(vec3<f32>(l)*shadow_value, 1.0);
+                                                          &dummy2));
+           }
+           let sun_direction =
+             (uniforms.camera_to_world_matrix
+              * vec4<f32>(normalize(sun_direction_samples[0]), 0.0)).xyz;
+           let lambertian_value = dot(hit_normal, sun_direction);
+           let ambient_strength = 0.25;
+           let l = ambient_strength + (1.0 - ambient_strength) * lambertian_value * shadow_value;
+           color_accumulator += vec4<f32>(vec3<f32>(l), 1.0);
         }
     }
-    if color_accumulator.w == 0.0 {
+    if color_accumulator.a == 0.0 {
         discard;
     }
-    return color_accumulator * 0.25;
+    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() % 25.0) / 25.0;
+    let animation_phase = 2.0 * std::f32::consts::PI * (animation_time.as_secs_f32() % 100.0) / 100.0;
     glam::Vec3::new(
         dem_size * f32::sin(animation_phase),
         dem_size * f32::cos(animation_phase),

src/app/dem_renderer/ray_intersection.wgsl

@@ -61,7 +61,8 @@ fn pop_node_stack(stack: ptr<function,NodeStack>) -> BoundingNode {
     return (*stack).stack[(*stack).count];
 }
 
-fn intersect_ray_with_node(tree_texture: texture_2d<u32>,
+fn intersect_ray_with_node(dem_texture: texture_2d<u32>,
+                           tree_texture: texture_2d<u32>,
                            dem_min_corner: vec2<f32>,
                            dem_cell_size: vec2<f32>,
                            dem_z_range: vec2<f32>,
@@ -154,17 +155,19 @@ fn intersect_ray_with_grid_cell(ray: Ray,
     // Calculate xyz of cell corners
     let xy00 = dem_min_corner.xy + dem_cell_size * vec2<f32>(cell_index);
     let p00 = vec3<f32>(xy00, z00);
-    let p01 = vec3<f32>(xy00 + vec2<f32>(0, 1), z01);
+    let p01 = vec3<f32>(xy00 + vec2<f32>(0, 1) * dem_cell_size, z01);
-    let p10 = vec3<f32>(xy00 + vec2<f32>(1, 0), z10);
+    let p10 = vec3<f32>(xy00 + vec2<f32>(1, 0) * dem_cell_size, z10);
-    let p11 = vec3<f32>(xy00 + vec2<f32>(1, 1), z11);
+    let p11 = vec3<f32>(xy00 + vec2<f32>(1, 1) * 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.
-    let p_t0 = intersect_ray_with_triangle_plane(ray, p00, p01, p11);
+    var hit_t0: bool;
-    let p_t1 = intersect_ray_with_triangle_plane(ray, p00, p10, p11);
+    let p_t0 = intersect_ray_with_triangle_plane(ray, p00, p01, p11, &hit_t0);
-    let point_is_in_triangle_0 = point_is_in_triangle(p_t0, p00, p01, p11);
+    var hit_t1: bool;
-    let point_is_in_triangle_1 = point_is_in_triangle(p_t0, p00, p10, p11);
+    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
@@ -183,9 +186,12 @@ fn intersect_ray_with_grid_cell(ray: Ray,
 fn intersect_ray_with_triangle_plane(ray: Ray,
                                      p0: vec3<f32>,
                                      p1: vec3<f32>,
-                                     p2: vec3<f32>) -> vec3<f32> {
+                                     p2: vec3<f32>,
+                                     hit: ptr<function,bool>) -> vec3<f32> {
     let n = cross(p1-p0, p2-p0);
-    return ray.origin + ray.direction * dot(p0 - ray.origin, n) / dot(ray.direction, n);
+    let t = dot(p0 - ray.origin, n) / dot(ray.direction, n);
+    *hit = t >= 0.0;
+    return ray.origin + ray.direction * t;
 }
 
 

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<u32> = bytemuck::cast_slice(&block_on(run_compute_shader_test(
+    let output_buffer: Vec<f32> = 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_eq!(output_buffer[0], 1);
+    assert!(dbg!(output_buffer[0]) >= 0.0);
-    assert_eq!(output_buffer[1], 0);
+    assert!(dbg!(output_buffer[1]) < 0.0);
-    assert_eq!(output_buffer[2], 1);
+    assert!(dbg!(output_buffer[2]) >= 0.0);
 }

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<u32>;
+var<storage, read_write> output_data: array<f32>;
 
 //#include ray_intersection.wgsl
 
@@ -24,10 +24,11 @@ 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;
-      if intersect_ray_with_aabb(ray, aabb) >= 0.0 {
+      output_data[i] = intersect_ray_with_aabb(ray, aabb);
+      /*if intersect_ray_with_aabb(ray, aabb) >= 0.0 {
           output_data[i] = 1u;
       }else {
           output_data[i] = 0u;
-      }
+          }*/
     }
 }