Calculate normals and apply lambertian lighting

src/app/dem_renderer/dem_renderer.wgsl

@@ -4,7 +4,8 @@ struct VertexOutput {
 }
 
 struct Uniforms {
-    transform: mat4x4<f32>,
+    camera_to_world_matrix: mat4x4<f32>,
+    world_to_ndc_matrix: mat4x4<f32>,
     camera_position: vec4<f32>,
     dem_min_corner: vec2<f32>,
     dem_cell_size: vec2<f32>,
@@ -27,7 +28,7 @@ fn vs_main(
     @location(0) position: vec4<f32>,
 ) -> VertexOutput {
     var result: VertexOutput;
-    result.position = uniforms.transform * position;
+    result.position = uniforms.world_to_ndc_matrix * position;
     result.world_space_position = position;
     return result;
 }
@@ -50,8 +51,41 @@ fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
                                ray,
                                root_node,
                                &hit_index) {
-        let v_i = textureLoad(dem_texture, hit_index, 0).r;
-        return vec4<f32>(vec3<f32>(v_i) / 65535.0, 1.0);
+        //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;
     }

src/app/dem_renderer/mod.rs

@@ -62,15 +62,20 @@ impl Camera {
             glam::Mat4::perspective_rh(std::f32::consts::FRAC_PI_4, ratio, 1.0, 10.0);
     }
 
-    fn get_matrix(&self) -> glam::Mat4 {
+    fn get_world_to_ndc_matrix(&self) -> glam::Mat4 {
         self.projection_matrix * self.view_matrix
     }
+
+    fn get_camera_to_world_matrix(&self) -> glam::Mat4 {
+        self.view_matrix.inverse()
+    }
 }
 
 #[repr(C)]
 #[derive(Clone, Copy, Pod, Zeroable)]
 struct UniformBufferContents {
-    camera_matrix: [f32; 16],
+    camera_to_world_matrix: [f32; 16],
+    world_to_ndc_matrix: [f32; 16],
     camera_position: [f32; 4],
     dem_min_corner: [f32; 2],
     dem_cell_size: [f32; 2],
@@ -80,7 +85,7 @@ struct UniformBufferContents {
 // member with the largest alignment. (The mat4x4<f32> type has a 16-byte
 // alignment.) So this value may be rounded up from the actual size of
 // UniformBufferContents.
-const UNIFORM_BUFFER_SIZE: u64 = 112;
+const UNIFORM_BUFFER_SIZE: u64 = 176;
 
 struct UniformBufferManager {
     cpu_buffer: UniformBufferContents,
@@ -101,10 +106,16 @@ impl UniformBufferManager {
         }
     }
 
-    fn set_camera_matrix(&mut self, camera_matrix: glam::Mat4) {
+    fn set_camera_to_world_matrix(&mut self, camera_to_world_matrix: glam::Mat4) {
         self.cpu_buffer
-            .camera_matrix
-            .clone_from_slice(&camera_matrix.to_cols_array())
+            .camera_to_world_matrix
+            .clone_from_slice(&camera_to_world_matrix.to_cols_array())
+    }
+
+    fn set_world_to_ndc_matrix(&mut self, world_to_ndc_matrix: glam::Mat4) {
+        self.cpu_buffer
+            .world_to_ndc_matrix
+            .clone_from_slice(&world_to_ndc_matrix.to_cols_array())
     }
 
     fn set_camera_position(&mut self, camera_position: glam::Vec4) {
@@ -301,7 +312,8 @@ impl DemRenderer {
         );
         self.camera
             .set_position(camera_position, self.get_dem_centre());
-        self.uniforms.set_camera_matrix(self.camera.get_matrix());
+        self.uniforms.set_camera_to_world_matrix(self.camera.get_camera_to_world_matrix());
+        self.uniforms.set_world_to_ndc_matrix(self.camera.get_world_to_ndc_matrix());
         self.uniforms
             .set_camera_position(camera_position.extend(1.0));
         self.uniforms.update_buffer(queue);
@@ -482,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() % 10.0) / 10.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

@@ -114,3 +114,26 @@ fn intersect_ray_with_node(tree_texture: texture_2d<u32>,
     }
     return closest_hit_distance < 1.0e+20;
 }
+
+
+fn intersect_ray_with_triangle_plane(ray: Ray,
+                                     p0: vec3<f32>,
+                                     p1: vec3<f32>,
+                                     p2: vec3<f32>) -> vec3<f32> {
+    let n = cross(p1-p0, p2-p0);
+    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;
+}
+
+
+fn point_is_in_triangle(p: vec3<f32>,
+                        p0: vec3<f32>, p1: vec3<f32>, p2: vec3<f32>) -> bool {
+  return point_is_inside_triangle_side(p, p0, p1, p2)
+    && point_is_inside_triangle_side(p, p1, p2, p0)
+    && point_is_inside_triangle_side(p, p2, p0, p1);
+}