Switch to nicer default scale

Implemented basic raycasting but it doesn't work; I'll check this in and
then add GPU unit tests.

Cargo.lock

@@ -1,6 +1,6 @@
 # This file is automatically @generated by Cargo.
 # It is not intended for manual editing.
-version = 3
+version = 4
 
 [[package]]
 name = "ab_glyph"
@@ -1450,9 +1450,9 @@ dependencies = [
 
 [[package]]
 name = "proc-macro2"
-version = "1.0.89"
+version = "1.0.92"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f139b0662de085916d1fb67d2b4169d1addddda1919e696f3252b740b629986e"
+checksum = "37d3544b3f2748c54e147655edb5025752e2303145b5aefb3c3ea2c78b973bb0"
 dependencies = [
  "unicode-ident",
 ]
@@ -1501,6 +1501,7 @@ dependencies = [
  "wasm-bindgen-test",
  "web-sys",
  "wgpu",
+ "wgsl-shader-assembler",
  "winit",
 ]
 
@@ -1804,9 +1805,9 @@ checksum = "6637bab7722d379c8b41ba849228d680cc12d0a45ba1fa2b48f2a30577a06731"
 
 [[package]]
 name = "syn"
-version = "2.0.87"
+version = "2.0.90"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "25aa4ce346d03a6dcd68dd8b4010bcb74e54e62c90c573f394c46eae99aba32d"
+checksum = "919d3b74a5dd0ccd15aeb8f93e7006bd9e14c295087c9896a110f490752bcf31"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -2331,6 +2332,17 @@ dependencies = [
  "web-sys",
 ]
 
+[[package]]
+name = "wgsl-shader-assembler"
+version = "0.1.0"
+source = "git+https://git.gordon.earth/matthew/wgsl-shader-assembler.git#2ce14e25d8f1685bd3486257cba4025587fd0fb9"
+dependencies = [
+ "naga",
+ "quote",
+ "regex",
+ "syn",
+]
+
 [[package]]
 name = "winapi-util"
 version = "0.1.9"

Cargo.toml

@@ -11,6 +11,7 @@ bytemuck = { version = "1.19.0", features = ["derive"] }
 glam = "0.29.2"
 log = "0.4.22"
 tiff = "0.9.1"
+wgsl-shader-assembler = { git = "https://git.gordon.earth/matthew/wgsl-shader-assembler.git" }
 
 [target.'cfg(target_arch = "x86_64")'.dependencies]
 winit = { version = "0.30.3", features = ["rwh_06"] }

rust-toolchain.toml

@@ -0,0 +1,2 @@
+[toolchain]
+channel = "nightly"

src/app/dem_renderer.wgsl

@@ -1,38 +0,0 @@
-struct VertexOutput {
-    @location(0) test: vec2<f32>,
-    @builtin(position) position: vec4<f32>,
-};
-
-struct Uniforms {
-    transform: mat4x4<f32>,
-    dem_texture_size: vec2<u32>,
-    dembvh_texture_size: vec2<u32>
-}
-
-@group(0) @binding(0)
-var<uniform> uniforms: Uniforms;
-
-@group(0) @binding(1)
-var dem_texture: texture_2d<u32>;
-
-@group(0) @binding(2)
-var dembvh_texture: texture_2d<u32>;
-
-
-@vertex
-fn vs_main(
-    @location(0) position: vec4<f32>,
-) -> VertexOutput {
-    var result: VertexOutput;
-    result.position = uniforms.transform * position;
-    result.test = position.xy;
-    return result;
-}
-
-@fragment
-fn fs_solid(vertex: VertexOutput) -> @location(0) vec4<f32> {
-    let texture_index = vec2<i32>(fract(abs(vertex.test)/500.0)*vec2<f32>(uniforms.dem_texture_size));
-    let v_i = vec3<u32>(textureLoad(dem_texture, texture_index, 0).r, textureLoad(dembvh_texture, texture_index, 0).rg);
-    let v = vec3<f32>(v_i) / 65535.0;
-    return vec4<f32>(v, 1.0);
-}

src/app/dem_renderer/dem_renderer.wgsl

@@ -0,0 +1,92 @@
+struct VertexOutput {
+    @location(0) world_space_position: vec4<f32>,
+    @builtin(position) position: vec4<f32>,
+}
+
+struct Uniforms {
+    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>,
+    dem_z_range: vec2<f32>
+}
+
+@group(0) @binding(0)
+var<uniform> uniforms: Uniforms;
+
+@group(0) @binding(1)
+var dem_texture: texture_2d<u32>;
+
+@group(0) @binding(2)
+var dembvh_texture: texture_2d<u32>;
+
+//#include ray_intersection.wgsl
+
+@vertex
+fn vs_main(
+    @location(0) position: vec4<f32>,
+) -> VertexOutput {
+    var result: VertexOutput;
+    result.position = uniforms.world_to_ndc_matrix * position;
+    result.world_space_position = position;
+    return result;
+}
+
+@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);
+
+    var root_node: BoundingNode;
+    root_node.index = vec2<u32>(0);
+    root_node.level = textureNumLevels(dembvh_texture) - 1;
+
+    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;
+    }
+}

src/app/dem_renderer/mod.rs

@@ -1,3 +1,5 @@
+use wgsl_shader_assembler::wgsl_module;
+
 use super::raster::{Dem, DemBvh};
 use {
     bytemuck::{Pod, Zeroable},
@@ -42,7 +44,7 @@ impl Camera {
             std::f32::consts::FRAC_PI_4,
             viewport_aspect_ratio,
             1.0,
-            10000.0,
+            30000.0,
         );
         let view_matrix = glam::Mat4::look_at_rh(glam::Vec3::ZERO, glam::Vec3::ZERO, glam::Vec3::Z);
         Self {
@@ -60,22 +62,30 @@ 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],
-    dem_texture_size: [u32; 2],
+    world_to_ndc_matrix: [f32; 16],
-    dembvh_texture_size: [u32; 2],
+    camera_position: [f32; 4],
+    dem_min_corner: [f32; 2],
+    dem_cell_size: [f32; 2],
+    dem_z_range: [f32; 2],
 }
-// This struct is 72 bytes but WGSL expects the buffer size to be be a multiple
+// WGSL expects the buffer size to be be a multiple of the alignment of the
-// of the alignment of the member with the largest alignment. The mat4x4<f32>
+// member with the largest alignment. (The mat4x4<f32> type has a 16-byte
-// type has a 16-byte alignment, so we round up to 16 * 5 = 80.
+// alignment.) So this value may be rounded up from the actual size of
-const UNIFORM_BUFFER_SIZE: u64 = 80;
+// UniformBufferContents.
+const UNIFORM_BUFFER_SIZE: u64 = 176;
 
 struct UniformBufferManager {
     cpu_buffer: UniformBufferContents,
@@ -96,22 +106,40 @@ 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
+            .camera_to_world_matrix
-            .clone_from_slice(&camera_matrix.to_cols_array())
+            .clone_from_slice(&camera_to_world_matrix.to_cols_array())
     }
 
-    fn set_dem_texture_size(&mut self, dem_texture_size: glam::UVec2) {
+    fn set_world_to_ndc_matrix(&mut self, world_to_ndc_matrix: glam::Mat4) {
         self.cpu_buffer
-            .dem_texture_size
+            .world_to_ndc_matrix
-            .clone_from_slice(&dem_texture_size.to_array());
+            .clone_from_slice(&world_to_ndc_matrix.to_cols_array())
     }
 
-    fn set_dembvh_texture_size(&mut self, dembvh_texture_size: glam::UVec2) {
+    fn set_camera_position(&mut self, camera_position: glam::Vec4) {
         self.cpu_buffer
-            .dembvh_texture_size
+            .camera_position
-            .clone_from_slice(&dembvh_texture_size.to_array());
+            .clone_from_slice(&camera_position.to_array())
+    }
+
+    fn set_dem_min_corner(&mut self, value: glam::Vec2) {
+        self.cpu_buffer
+            .dem_min_corner
+            .clone_from_slice(&value.to_array());
+    }
+
+    fn set_dem_cell_size(&mut self, value: glam::Vec2) {
+        self.cpu_buffer
+            .dem_cell_size
+            .clone_from_slice(&value.to_array());
+    }
+
+    fn set_dem_z_range(&mut self, value: glam::Vec2) {
+        self.cpu_buffer
+            .dem_z_range
+            .clone_from_slice(&value.to_array());
     }
 
     fn get_binding(&self) -> wgpu::BindingResource<'_> {
@@ -150,15 +178,19 @@ impl DemRenderer {
 
         let index_count = index_data.len();
 
-        let (dem_texture_view, dem_texture_size) = load_dem_texture(&source, device, queue);
+        let dem_texture_view = load_dem_texture(&source, device, queue);
-        let (dembvh_texture_view, dembvh_texture_size) =
+        let dembvh_texture_view = create_dembvh_texture(&source, device, queue);
-            create_dembvh_texture(&source, device, queue);
 
         let camera = Camera::new(surface_config.width as f32 / surface_config.height as f32);
 
         let mut uniforms = UniformBufferManager::new(device);
-        uniforms.set_dem_texture_size(dem_texture_size);
+        uniforms.set_dem_min_corner(glam::Vec2::new(source.x_min, source.y_min));
-        uniforms.set_dembvh_texture_size(dembvh_texture_size);
+        uniforms.set_dem_cell_size(
+            (glam::Vec2::new(source.x_max, source.y_max)
+                - glam::Vec2::new(source.x_min, source.y_min))
+                / glam::Vec2::new(source.num_cells_x as f32, source.num_cells_y as f32),
+        );
+        uniforms.set_dem_z_range(glam::Vec2::new(source.z_min, source.z_max));
 
         let bind_group_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
             label: None,
@@ -231,10 +263,7 @@ impl DemRenderer {
             }],
         }];
 
-        let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
+        let shader = device.create_shader_module(wgsl_module!("dem_renderer.wgsl"));
-            label: None,
-            source: wgpu::ShaderSource::Wgsl(Cow::Borrowed(include_str!("dem_renderer.wgsl"))),
-        });
 
         let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
             label: Some("DemRendererPipeline"),
@@ -277,14 +306,16 @@ impl DemRenderer {
         let mut encoder = device.create_command_encoder(&wgpu::CommandEncoderDescriptor {
             label: Some("DemRendererCommandEncoder"),
         });
-        self.camera.set_position(
+        let camera_position = get_animated_camera_position(
-            get_animated_camera_position(
+            self.animation_start.elapsed(),
-                self.animation_start.elapsed(),
+            self.get_max_dem_dimension(),
-                self.get_max_dem_dimension(),
-            ),
-            self.get_dem_centre(),
         );
-        self.uniforms.set_camera_matrix(self.camera.get_matrix());
+        self.camera
+            .set_position(camera_position, self.get_dem_centre());
+        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);
         {
             let mut rpass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor {
@@ -378,11 +409,7 @@ fn create_vertices(dem: &Rc<Dem>) -> (Vec<Vertex>, Vec<u16>) {
     (vertex_data.to_vec(), index_data.to_vec())
 }
 
-fn load_dem_texture(
+fn load_dem_texture(source: &Dem, device: &wgpu::Device, queue: &wgpu::Queue) -> wgpu::TextureView {
-    source: &Dem,
-    device: &wgpu::Device,
-    queue: &wgpu::Queue,
-) -> (wgpu::TextureView, glam::UVec2) {
     let texture_size = wgpu::Extent3d {
         width: source.num_cells_x,
         height: source.num_cells_y,
@@ -415,17 +442,14 @@ fn load_dem_texture(
         texture_size,
     );
 
-    (
+    texture.create_view(&wgpu::TextureViewDescriptor::default())
-        texture.create_view(&wgpu::TextureViewDescriptor::default()),
-        glam::UVec2::new(source.num_cells_x, source.num_cells_y),
-    )
 }
 
 fn create_dembvh_texture(
     dem: &Dem,
     device: &wgpu::Device,
     queue: &wgpu::Queue,
-) -> (wgpu::TextureView, glam::UVec2) {
+) -> wgpu::TextureView {
     let source = DemBvh::new(dem);
     let dembvh_size = source.layers.first().expect("DEM BVH contains layers").size;
 
@@ -466,17 +490,17 @@ fn create_dembvh_texture(
         );
     }
 
-    (
+    texture.create_view(&wgpu::TextureViewDescriptor::default())
-        texture.create_view(&wgpu::TextureViewDescriptor::default()),
-        glam::UVec2::new(dembvh_size, dembvh_size),
-    )
 }
 
 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),
         dem_size * 0.7,
     )
 }
+
+#[cfg(test)]
+mod tests;

src/app/dem_renderer/pack_unpack.wgsl

@@ -0,0 +1,3 @@
+fn scale_u16(value: u32, range: vec2<f32>) -> f32 {
+    return range.x + (f32(value - 1) / 65534.0) * (range.y - range.x);
+}

src/app/dem_renderer/ray_intersection.wgsl

@@ -0,0 +1,139 @@
+//#include pack_unpack.wgsl
+
+struct BoundingNode {
+    index: vec2<u32>,
+    level: u32
+}
+
+struct AABB {
+    min_corner: vec3<f32>,
+    max_corner: vec3<f32>
+}
+
+struct Ray { origin: vec3<f32>, direction: vec3<f32> };
+
+fn invert_ray_direction(v: vec3<f32> ) -> vec3<f32> {
+    return select(vec3<f32>(1.0e+30),
+                  vec3<f32>(1.0) / v,
+                  vec3<bool>(v));
+}
+
+fn intersect_ray_with_aabb(ray: Ray, aabb: AABB) -> f32 {
+  return intersect_ray_with_aabb_optimized(ray.origin,
+                                           invert_ray_direction(ray.direction),
+                                           aabb);
+}
+
+fn intersect_ray_with_aabb_optimized(ray_origin: vec3<f32>, inv_ray_direction: vec3<f32>, aabb: AABB) -> f32 {
+    let t1 = (aabb.min_corner - ray_origin) * inv_ray_direction;
+    let t2 = (aabb.max_corner - ray_origin) * inv_ray_direction;
+    let t_mins = min(t1, t2);
+    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);
+}
+
+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(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<BoundingNode,64>,
+    count: u32
+}
+
+fn push_node_stack(stack: ptr<function,NodeStack>, node: BoundingNode) {
+    (*stack).stack[(*stack).count] = node;
+    (*stack).count++;
+}
+
+fn pop_node_stack(stack: ptr<function,NodeStack>) -> BoundingNode {
+    (*stack).count--;
+    return (*stack).stack[(*stack).count];
+}
+
+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_cell: ptr<function, vec2<u32>>) -> bool {
+    let inv_ray_direction = invert_ray_direction(ray.direction);
+    var node_stack: NodeStack;
+    node_stack.count = 0u;
+    push_node_stack(&node_stack, root_node);
+    var closest_hit_distance = 1.0e+30f;
+    while node_stack.count > 0 {
+        let node = pop_node_stack(&node_stack);
+        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, 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(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 {
+                if hit_distance < closest_hit_distance {
+                    closest_hit_distance = hit_distance;
+                    *hit_cell = node.index;
+                }
+            } else {
+                let next_index = node.index * 2;
+                var next_node: BoundingNode;
+                next_node.index = next_index;
+                next_node.level = node.level - 1;
+                push_node_stack(&node_stack, next_node);
+                next_node.index = next_index + vec2<u32>(1, 0);
+                push_node_stack(&node_stack, next_node);
+                next_node.index = next_index + vec2<u32>(0, 1);
+                push_node_stack(&node_stack, next_node);
+                next_node.index = next_index + vec2<u32>(1, 1);
+                push_node_stack(&node_stack, next_node);
+            }
+        }
+    }
+    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);
+}

src/app/dem_renderer/tests.rs

@@ -0,0 +1,153 @@
+use futures::executor::block_on;
+use std::sync::mpsc::channel;
+
+use wgsl_shader_assembler::wgsl_module;
+
+use super::*;
+
+async fn run_compute_shader_test(
+    shader_module: wgpu::ShaderModuleDescriptor<'_>,
+    test_function: impl AsRef<str>,
+    input: &[u8],
+    output_size: usize,
+) -> Vec<u8> {
+    let wgpu_instance = wgpu::Instance::default();
+    let adapter = wgpu_instance
+        .request_adapter(&wgpu::RequestAdapterOptions::default())
+        .await
+        .unwrap();
+    let (device, queue) = adapter
+        .request_device(
+            &wgpu::DeviceDescriptor {
+                label: None,
+                required_features: wgpu::Features::empty(),
+                required_limits: wgpu::Limits::downlevel_defaults(),
+                memory_hints: wgpu::MemoryHints::MemoryUsage,
+            },
+            None,
+        )
+        .await
+        .unwrap();
+
+    let shader_module = device.create_shader_module(shader_module);
+    let staging_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+        label: None,
+        size: output_size as u64,
+        usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST,
+        mapped_at_creation: false,
+    });
+    let input_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+        label: None,
+        contents: input,
+        usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST,
+    });
+    let output_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+        label: None,
+        size: output_size as u64,
+        usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
+        mapped_at_creation: false,
+    });
+    let compute_pipeline = device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+        label: None,
+        layout: None,
+        module: &shader_module,
+        entry_point: Some(test_function.as_ref()),
+        compilation_options: Default::default(),
+        cache: None,
+    });
+    let bind_group_layout = compute_pipeline.get_bind_group_layout(0);
+    let bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
+        label: None,
+        layout: &bind_group_layout,
+        entries: &[
+            wgpu::BindGroupEntry {
+                binding: 0,
+                resource: input_buffer.as_entire_binding(),
+            },
+            wgpu::BindGroupEntry {
+                binding: 1,
+                resource: output_buffer.as_entire_binding(),
+            },
+        ],
+    });
+    let mut encoder =
+        device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: None });
+    {
+        let mut cpass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+            label: None,
+            timestamp_writes: None,
+        });
+        cpass.set_pipeline(&compute_pipeline);
+        cpass.set_bind_group(0, &bind_group, &[]);
+        cpass.dispatch_workgroups(1, 1, 1);
+    }
+    encoder.copy_buffer_to_buffer(&output_buffer, 0, &staging_buffer, 0, output_size as u64);
+    queue.submit(Some(encoder.finish()));
+
+    let (tx, rx) = channel();
+    let buffer_slice = staging_buffer.slice(..);
+    buffer_slice.map_async(wgpu::MapMode::Read, move |v| tx.send(v).unwrap());
+    device.poll(wgpu::Maintain::wait()).panic_on_timeout();
+    if let Ok(Ok(())) = rx.recv() {
+        buffer_slice.get_mapped_range().to_vec()
+    } else {
+        panic!("failed to run test on GPU")
+    }
+}
+
+#[repr(C, align(16))]
+#[derive(Clone, Copy, Pod, Zeroable)]
+struct Vec3 {
+    elements: [f32; 3],
+    _padding: f32
+}
+
+#[repr(C)]
+#[derive(Clone, Copy, Pod, Zeroable)]
+struct TestInput {
+    ray_origin: Vec3,
+    ray_direction: Vec3,
+    aabb_min_corner: Vec3,
+    aabb_max_corner: Vec3,
+}
+
+fn vec3(x: f32, y: f32, z: f32) -> Vec3 {
+    Vec3 {
+        elements: [x, y, z],
+        _padding: 0.0
+    }
+}
+
+#[test]
+fn test_shaders() {
+    let input_buffer = vec![
+        TestInput {
+            ray_origin: vec3(-5.0, 0.0, 0.0),
+            ray_direction: vec3(1.0, 0.0, 0.0),
+            aabb_min_corner: vec3(-1.0, -1.0, -1.0),
+            aabb_max_corner: vec3(1.0, 1.0, 1.0),
+        },
+        TestInput {
+            ray_origin: vec3(-5.0, 0.0, 0.0),
+            ray_direction: vec3(1.0, 0.0, 0.0),
+            aabb_min_corner: vec3(-1.0, 0.5, -1.0),
+            aabb_max_corner: vec3(1.0, 1.0, 1.0),
+        },
+        TestInput {
+            ray_origin: vec3(-5.0, 0.0, 0.0),
+            ray_direction: vec3(1.0, 0.15, 0.0),
+            aabb_min_corner: vec3(-1.5, 0.5, -1.0),
+            aabb_max_corner: vec3(1.0, 1.0, 1.0),
+        },
+    ];
+    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_eq!(output_buffer[0], 1);
+    assert_eq!(output_buffer[1], 0);
+    assert_eq!(output_buffer[2], 1);
+}

src/app/dem_renderer/tests.wgsl

@@ -0,0 +1,33 @@
+@group(0)
+@binding(0)
+var<storage, read> input_data: array<Input>;
+@group(0)
+@binding(1)
+var<storage, read_write> output_data: array<u32>;
+
+//#include ray_intersection.wgsl
+
+struct Input {
+    ray_origin: vec3<f32>,
+    ray_direction: vec3<f32>,
+    aabb_min_corner: vec3<f32>,
+    aabb_max_corner: vec3<f32>
+}
+
+@compute
+@workgroup_size(1)
+fn test_intersect_ray_with_aabb() {
+    for(var i=0u; i < arrayLength(&input_data); i++) {
+      var ray: Ray;
+      ray.origin = input_data[i].ray_origin;
+      ray.direction = input_data[i].ray_direction;
+      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] = 1u;
+      }else {
+          output_data[i] = 0u;
+      }
+    }
+}

src/app/mod.rs

@@ -46,12 +46,13 @@ impl MvuApp<Model> for App {
     async fn init(&mut self, instance: &Instance, surface: Surface<'static>, size: Size2i) {
         let adapter = instance
             .request_adapter(&wgpu::RequestAdapterOptions {
-                power_preference: wgpu::PowerPreference::default(),
+                power_preference: wgpu::PowerPreference::HighPerformance,
                 force_fallback_adapter: false,
                 compatible_surface: Some(&surface),
             })
             .await
             .expect("Failed to find an appropriate adapter");
+        eprintln!("Using {}", adapter.get_info().name);
 
         let (device, queue) = adapter
             .request_device(

src/app/raster/dem.rs

@@ -174,10 +174,10 @@ 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_min = -5000.0;
-                let x_max = 500.0;
+                let x_max = 5000.0;
-                let y_min = -500.0;
+                let y_min = -5000.0;
-                let y_max = 500.0;
+                let y_max = 5000.0;
                 let (z_min, z_max) = f32_values[1..]
                     .iter()
                     .fold((f32_values[0], f32_values[0]), |(min, max), elem| {