Block-rendering based slicing

tests/test_yt_idv.py

@@ -61,6 +61,14 @@ def test_snapshots(osmesa_fake_amr, image_store):
     image_store(osmesa_fake_amr)
 
 
+def test_slice(osmesa_fake_amr, image_store):
+    osmesa_fake_amr.scene.components[0].render_method = "slice"
+    image_store(osmesa_fake_amr)
+    osmesa_fake_amr.scene.components[0].slice_normal = (1.0, 1.0, 0.0)
+    osmesa_fake_amr.scene.components[0].slice_position = (0.5, 0.25, 0.5)
+    image_store(osmesa_fake_amr)
+
+
 def test_annotate_boxes(osmesa_empty, image_store):
     """Check the box annotation."""
     osmesa_empty.scene.add_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])

yt_idv/scene_components/blocks.py

@@ -4,6 +4,7 @@
 import traitlets
 from OpenGL import GL
 
+from yt_idv.gui_support import add_popup_help
 from yt_idv.opengl_support import TransferFunctionTexture
 from yt_idv.scene_components.base_component import SceneComponent
 from yt_idv.scene_data.block_collection import BlockCollection
@@ -26,6 +27,8 @@ class BlockRendering(SceneComponent):
     tf_min = traitlets.CFloat(0.0)
     tf_max = traitlets.CFloat(1.0)
     tf_log = traitlets.Bool(True)
+    slice_position = traitlets.Tuple((0.5, 0.5, 0.5), trait=traitlets.CFloat())
+    slice_normal = traitlets.Tuple((1.0, 0.0, 0.0), trait=traitlets.CFloat())
 
     priority = 10
 
@@ -103,6 +106,21 @@ def render_gui(self, imgui, renderer, scene):
                     data[xb1:xb2, 0, i] = np.mgrid[yv1 : yv2 : (xb2 - xb1) * 1j]
             if update:
                 self.transfer_function.data = (data * 255).astype("u1")
+
+        elif self.render_method == "slice":
+            imgui.text("Set slicing parameters:")
+
+            _, self.slice_position = imgui.input_float3(
+                "Position", *self.slice_position
+            )
+            changed = changed or _
+            _ = add_popup_help(imgui, "The position of a point on the slicing plane.")
+            changed = changed or _
+            _, self.slice_normal = imgui.input_float3("Normal", *self.slice_normal)
+            changed = changed or _
+            _ = add_popup_help(imgui, "The normal vector of the slicing plane.")
+            changed = changed or _
+
         return changed
 
     @traitlets.default("transfer_function")
@@ -129,6 +147,8 @@ def _set_uniforms(self, scene, shader_program):
         shader_program._set_uniform("tf_min", self.tf_min)
         shader_program._set_uniform("tf_max", self.tf_max)
         shader_program._set_uniform("tf_log", float(self.tf_log))
+        shader_program._set_uniform("slice_normal", np.array(self.slice_normal))
+        shader_program._set_uniform("slice_position", np.array(self.slice_position))
 
     def _get_sanitized_iso_layers(self):
         # return the sanitized layers

yt_idv/shaders/known_uniforms.inc.glsl

@@ -27,6 +27,10 @@ uniform int channel;
 // Mesh rendering
 uniform mat4 model_to_clip;
 
+// Slicing
+uniform vec3 slice_position;
+uniform vec3 slice_normal;
+
 // Matrices for projection and positions
 uniform mat4 modelview;
 uniform mat4 projection;

yt_idv/shaders/shaderlist.yaml

@@ -125,6 +125,13 @@ shader_definitions:
       blend_func:
       - one
       - one
+    slice_sample:
+      info: Slice through a block collection
+      source: slice_sample.frag.glsl
+      depth_test: less
+      blend_func:
+      - one
+      - zero
       blend_equation: func add
   vertex:
     default:
@@ -207,6 +214,13 @@ component_shaders:
       first_fragment: isocontour
       second_vertex: passthrough
       second_fragment: apply_colormap
+    slice:
+      description: Slice
+      first_vertex: grid_position
+      first_geometry: grid_expand
+      first_fragment: slice_sample
+      second_vertex: passthrough
+      second_fragment: apply_colormap
   octree_block_rendering:
     default_value: max_intensity
     max_intensity:

yt_idv/shaders/slice_sample.frag.glsl

@@ -0,0 +1,66 @@
+in vec4 v_model;
+flat in vec3 dx;
+flat in vec3 left_edge;
+flat in vec3 right_edge;
+flat in mat4 inverse_proj;
+flat in mat4 inverse_mvm;
+flat in mat4 inverse_pmvm;
+out vec4 output_color;
+
+bool within_bb(vec3 pos)
+{
+    bvec3 left =  greaterThanEqual(pos, left_edge);
+    bvec3 right = lessThanEqual(pos, right_edge);
+    return all(left) && all(right);
+}
+
+bool sample_texture(vec3 tex_curr_pos, inout vec4 curr_color, float tdelta,
+                    float t, vec3 dir);
+vec4 cleanup_phase(in vec4 curr_color, in vec3 dir, in float t0, in float t1);
+
+// This main() function will call a function called sample_texture at every
+// step along the ray.  It must be of the form
+//   void (vec3 tex_curr_pos, inout vec4 curr_color, float tdelta, float t,
+//         vec3 direction);
+
+void main()
+{
+    // Obtain screen coordinates
+    // https://www.opengl.org/wiki/Compute_eye_space_from_window_space#From_gl_FragCoord
+    vec4 ndcPos;
+    ndcPos.xy = ((2.0 * gl_FragCoord.xy) - (2.0 * viewport.xy)) / (viewport.zw) - 1;
+    ndcPos.z = (2.0 * gl_FragCoord.z - 1.0);
+    ndcPos.w = 1.0;
+
+    vec4 clipPos = ndcPos / gl_FragCoord.w;
+    vec4 eyePos = inverse_proj * clipPos;
+    eyePos /= eyePos.w;
+
+    vec3 ray_position = (inverse_pmvm * clipPos).xyz;
+
+    // Five samples
+    vec3 dir = normalize(camera_pos.xyz - ray_position);
+    dir = max(abs(dir), 0.0001) * sign(dir);
+    vec4 curr_color = vec4(0.0);
+
+    // We'll compute the t at which this ray intersects the slice. If that t
+    // results in a position that is within this box, we'll sample and return.
+    // For a nice, rust-y walkthrough: https://samsymons.com/blog/math-notes-ray-plane-intersection/
+    float t_intersect = dot((slice_position - ray_position), slice_normal) / dot(dir, slice_normal);
+    if (abs(t_intersect) < 1e-5) discard;
+    ray_position += t_intersect * dir;
+    if (!within_bb(ray_position)) discard;
+
+    vec3 range = (right_edge + dx/2.0) - (left_edge - dx/2.0);
+    vec3 nzones = range / dx;
+    vec3 ndx = 1.0/nzones;
+
+    vec3 tex_curr_pos = (ray_position - left_edge) / range;  // Scale from 0 .. 1
+    // But, we actually need it to be 0 + normalized dx/2 to 1 - normalized dx/2
+    tex_curr_pos = (tex_curr_pos * (1.0 - ndx)) + ndx/2.0;
+
+    float map_sample = texture(bitmap_tex, tex_curr_pos).r;
+    if (!(map_sample > 0.0)) discard;
+
+    output_color = texture(ds_tex, tex_curr_pos);
+}