diff --git a/src/app/dem_renderer/dem_renderer.wgsl b/src/app/dem_renderer/dem_renderer.wgsl index c800575..7d1a5f3 100644 --- a/src/app/dem_renderer/dem_renderer.wgsl +++ b/src/app/dem_renderer/dem_renderer.wgsl @@ -44,48 +44,37 @@ fn fs_solid(vertex: VertexOutput) -> @location(0) vec4 { root_node.level = textureNumLevels(dembvh_texture) - 1; var hit_index = vec2(0); + var hit_location: vec3; + var hit_normal: vec3; 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(0,1), 0).r; - let v10 = textureLoad(dem_texture, hit_index + vec2(1,0), 0).r; - let v11 = textureLoad(dem_texture, hit_index + vec2(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(hit_index); - let p00 = vec3(xy00, z00); - let p01 = vec3(xy00 + vec2(0, 1), z01); - let p10 = vec3(xy00 + vec2(1, 0), z10); - let p11 = vec3(xy00 + vec2(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 + &hit_index, + &hit_location, + &hit_normal) { // Calculate light let sun_direction = (uniforms.camera_to_world_matrix * vec4(normalize(vec3(1.0, 1.0, 1.0)), 0.0)).xyz; - let l = dot(normal, sun_direction); - return vec4(l, l, l, 1.0); - //return vec4(normal.xy, 0.0, 1.0); + 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; + var dummy1: vec3; + 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(vec3(l)*shadow_value, 1.0); } else { discard; } diff --git a/src/app/dem_renderer/ray_intersection.wgsl b/src/app/dem_renderer/ray_intersection.wgsl index 2dd4725..51c3893 100644 --- a/src/app/dem_renderer/ray_intersection.wgsl +++ b/src/app/dem_renderer/ray_intersection.wgsl @@ -31,7 +31,7 @@ fn intersect_ray_with_aabb_optimized(ray_origin: vec3, 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, @@ -67,7 +67,9 @@ fn intersect_ray_with_node(tree_texture: texture_2d, dem_z_range: vec2, ray: Ray, root_node: BoundingNode, - hit_cell: ptr>) -> bool { + hit_cell: ptr>, + hit_location: ptr>, + hit_normal: ptr>) -> 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, 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; + var normal: vec3; 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, return closest_hit_distance < 1.0e+20; } +fn intersect_ray_with_grid_cell(ray: Ray, + cell_index: vec2, + dem_texture: texture_2d, + dem_min_corner: vec2, + dem_cell_size: vec2, + dem_z_range: vec2, + location: ptr>, + normal: ptr>) -> 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(0,1), 0).r; + let v10 = textureLoad(dem_texture, cell_index + vec2(1,0), 0).r; + let v11 = textureLoad(dem_texture, cell_index + vec2(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(cell_index); + let p00 = vec3(xy00, z00); + let p01 = vec3(xy00 + vec2(0, 1), z01); + let p10 = vec3(xy00 + vec2(1, 0), z10); + let p11 = vec3(xy00 + vec2(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, @@ -127,7 +191,7 @@ fn intersect_ray_with_triangle_plane(ray: Ray, fn point_is_inside_triangle_side(p: vec3, p0: vec3, p1: vec3, p2: vec3) -> 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; } diff --git a/src/app/raster/dem.rs b/src/app/raster/dem.rs index e932f49..f6033ed 100644 --- a/src/app/raster/dem.rs +++ b/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_max = 5000.0; - let y_min = -5000.0; - let y_max = 5000.0; + let x_min = -500.0; + let x_max = 500.0; + let y_min = -500.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,