Add hard shadows, although there are bad artifacts

src/app/dem_renderer/dem_renderer.wgsl

@@ -44,48 +44,37 @@ fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
     root_node.level = textureNumLevels(dembvh_texture) - 1;
 
     var hit_index = vec2<u32>(0);
+    var hit_location: vec3<f32>;
+    var hit_normal: vec3<f32>;
     if intersect_ray_with_node(dembvh_texture,
                                uniforms.dem_min_corner,
                                uniforms.dem_cell_size,
                                uniforms.dem_z_range,
                                ray,
                                root_node,
-                               &hit_index) {
+                               &hit_index,
-        //Calculate x-values of cell corners
+                               &hit_location,
-        let v00 = textureLoad(dem_texture, hit_index, 0).r;
+                               &hit_normal) {
-        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);
+        let l = dot(hit_normal, sun_direction);
-        return vec4<f32>(l, l, l, 1.0);
+        var shadow_ray :Ray;
-        //return vec4<f32>(normal.xy, 0.0, 1.0);
+        shadow_ray.origin = hit_location + hit_normal * 0.1;
+        shadow_ray.direction = sun_direction;
+        var dummy0: vec3<f32>;
+        var dummy1: vec3<f32>;
+        let shadow_value =
+          select(1.0, 0.0, intersect_ray_with_node(dembvh_texture,
+                                                   uniforms.dem_min_corner,
+                                                   uniforms.dem_cell_size,
+                                                   uniforms.dem_z_range,
+                                                   shadow_ray,
+                                                   root_node,
+                                                   &hit_index,
+                                                   &dummy0,
+                                                   &dummy1));
+        return vec4<f32>(vec3<f32>(l)*shadow_value, 1.0);
     } else {
         discard;
     }

src/app/dem_renderer/ray_intersection.wgsl

@@ -31,7 +31,7 @@ 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);
+    return select(-1.0, max(t_min, 0.0f),  t_min <= t_max && t_max > 0.0);
 }
 
 fn get_xy_min_for_node(dem_min_corner: vec2<f32>,
@@ -67,7 +67,9 @@ fn intersect_ray_with_node(tree_texture: texture_2d<u32>,
                            dem_z_range: vec2<f32>,
                            ray: Ray,
                            root_node: BoundingNode,
-                           hit_cell: ptr<function, vec2<u32>>) -> bool {
+                           hit_cell: ptr<function, vec2<u32>>,
+                           hit_location: ptr<function, vec3<f32>>,
+                           hit_normal: ptr<function, vec3<f32>>) -> bool {
     let inv_ray_direction = invert_ray_direction(ray.direction);
     var node_stack: NodeStack;
     node_stack.count = 0u;
@@ -93,9 +95,25 @@ fn intersect_ray_with_node(tree_texture: texture_2d<u32>,
         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,
+                                                    dem_texture,
+                                                    dem_min_corner,
+                                                    dem_cell_size,
+                                                    dem_z_range,
+                                                    &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_cell = node.index;
+                    *hit_location = location;
+                    *hit_normal = normal;
+                    }
                 }
             } else {
                 let next_index = node.index * 2;
@@ -115,6 +133,52 @@ fn intersect_ray_with_node(tree_texture: texture_2d<u32>,
     return closest_hit_distance < 1.0e+20;
 }
 
+fn intersect_ray_with_grid_cell(ray: Ray,
+                                cell_index: vec2<u32>,
+                                dem_texture: texture_2d<u32>,
+                                dem_min_corner: vec2<f32>,
+                                dem_cell_size: vec2<f32>,
+                                dem_z_range: vec2<f32>,
+                                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, dem_z_range);
+    let z01 = scale_u16(v01, dem_z_range);
+    let z10 = scale_u16(v10, dem_z_range);
+    let z11 = scale_u16(v11, dem_z_range);
+    // 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 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_triangle_0 = point_is_in_triangle(p_t0, p00, p01, p11);
+    let point_is_in_triangle_1 = point_is_in_triangle(p_t0, 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>,
@@ -127,7 +191,7 @@ fn intersect_ray_with_triangle_plane(ray: Ray,
 
 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;
+  return dot(cross(p2-p1, p-p1), cross(p2-p1, p0-p1)) >= -0.0;
 }
 
 

src/app/raster/dem.rs

@@ -174,18 +174,20 @@ 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 = -5000.0;
+                let x_min = -500.0;
-                let x_max = 5000.0;
+                let x_max = 500.0;
-                let y_min = -5000.0;
+                let y_min = -500.0;
-                let y_max = 5000.0;
+                let y_max = 500.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(|v| normalized_f32_to_u16(normalize_f32(*v, z_min, z_max)))
+                    .map(|z| z * 0.1)
+                    .map(|v| normalized_f32_to_u16(normalize_f32(v, z_min, z_max)))
                     .collect();
                 Dem {
                     num_cells_x,