Add voxel plot type (#3527)

* create voxel rendering prototype * enable lighting * prototype voxel id generation & color handling * add is_air attribute * prototype texture mapping * fix shader reloading * fix texture mapping * implement local updates * optimize render order (depthsorting = false) * add depthsorting = true * render z planes first * add lowclip and highclip * add refimg tests + some fixes * fix colorrange * fix local chunk update * handle colorrange more efficiently * handle voxel id data more efficiently * docstring & formatting * switch back to lrbt order for uvmap * add docs * try fix tests * fix show * fix test? * add missing dimensions * add arguments for placement and scale * allow Colon * add Colon() to local_update * minor cleanup * prototype WGLMakie version * add fallback in CairoMakie * add RPRMakie fallback * skip invisible voxels * fix typo * rename voxel -> voxels * update docs, fix placement * update news * fix Colorbar for voxels * enable tests * fix texture rotation * cleanup print * cleanup comment * generalize array access * debug WGLMakie * get voxels rendering in WGLMakie * fix texture mapping * activate tests * fix moving planes, cleanup prints * add unit tests * add gap attribute * tests & docs * mention potential issues with picking * fix WGLMakie picking * fix depthsorting/gap handling * switch to integer mod * fix render order * use RNG * fix 1.6 3d array syntax * fix refimage * Update CHANGELOG.md * fix julia 1.6 --------- Co-authored-by: Simon <[email protected]>
MakieOrg · Mar 8, 2024 · 274df26 · 274df26
1 parent 6681722
commit 274df26
Show file tree

Hide file tree

Showing 32 changed files with 1,601 additions and 24 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,10 +1,13 @@
 # Changelog
 
 ## [Unreleased]
+
 - Added supported markers hint to unsupported marker warn message.
+- Add `voxels` plot [#3527](https://github.com/MakieOrg/Makie.jl/pull/3527)
 
-- Remove StableHashTraits in favor of calculating hashes directly with CRC32c [#3667](https://github.com/MakieOrg/Makie.jl/pull/3667).
+## [0.21.0] - 2024-03-0X
 
+- Remove StableHashTraits in favor of calculating hashes directly with CRC32c [#3667](https://github.com/MakieOrg/Makie.jl/pull/3667).
 - **Breaking (sort of)** Added a new `@recipe` variant which allows documenting attributes directly where they are defined and validating that all attributes are known whenever a plot is created. This is not breaking in the sense that the API changes, but user code is likely to break because of misspelled attribute names etc. that have so far gone unnoticed.
 - **Breaking** Reworked line shaders in GLMakie and WGLMakie [#3558](https://github.com/MakieOrg/Makie.jl/pull/3558)
   - GLMakie: Removed support for per point linewidths

diff --git a/CairoMakie/src/primitives.jl b/CairoMakie/src/primitives.jl
@@ -1206,3 +1206,42 @@ function draw_atomic(scene::Scene, screen::Screen, @nospecialize(primitive::Maki
 
     return nothing
 end
+
+
+
+################################################################################
+#                                    Voxel                                     #
+################################################################################
+
+
+function draw_atomic(scene::Scene, screen::Screen, @nospecialize(primitive::Makie.Voxels))
+    pos = Makie.voxel_positions(primitive)
+    scale = Makie.voxel_size(primitive)
+    colors = Makie.voxel_colors(primitive)
+    marker = normal_mesh(Rect3f(Point3f(-0.5), Vec3f(1)))
+
+    # For correct z-ordering we need to be in view/camera or screen space
+    model = copy(primitive.model[])
+    view = scene.camera.view[]
+
+    zorder = sortperm(pos, by = p -> begin
+        p4d = to_ndim(Vec4f, to_ndim(Vec3f, p, 0f0), 1f0)
+        cam_pos = view * model * p4d
+        cam_pos[3] / cam_pos[4]
+    end, rev=false)
+
+    submesh = Attributes(
+        model=model,
+        shading=primitive.shading, diffuse=primitive.diffuse,
+        specular=primitive.specular, shininess=primitive.shininess,
+        faceculling=get(primitive, :faceculling, -10),
+        transformation=Makie.transformation(primitive)
+    )
+
+    for i in zorder
+        submesh[:calculated_colors] = colors[i]
+        draw_mesh3D(scene, screen, submesh, marker, pos = pos[i], scale = scale)
+    end
+
+    return nothing
+end
diff --git a/CairoMakie/test/runtests.jl b/CairoMakie/test/runtests.jl
@@ -185,6 +185,7 @@ excludes = Set([
     "scatter with stroke",
     "heatmaps & surface",
     "Textured meshscatter", # not yet implemented
+    "Voxel - texture mapping", # not yet implemented
     "Miter Joints for line rendering", # CairoMakie does not show overlap here and extrudes lines a little more
 ])
 

diff --git a/GLMakie/assets/shader/voxel.frag b/GLMakie/assets/shader/voxel.frag
@@ -0,0 +1,133 @@
+#version 330 core
+// {{GLSL VERSION}}
+// {{GLSL_EXTENSIONS}}
+
+// debug FLAGS
+// #define DEBUG_RENDER_ORDER 0 // (0, 1, 2) - dimensions
+
+struct Nothing{ //Nothing type, to encode if some variable doesn't contain any data
+    bool _; //empty structs are not allowed
+};
+
+// Sets which shading procedures to use
+{{shading}}
+
+flat in vec3 o_normal;
+in vec3 o_uvw;
+flat in int o_side;
+in vec2 o_tex_uv;
+
+#ifdef DEBUG_RENDER_ORDER
+flat in float plane_render_idx; // debug
+flat in int plane_dim;
+flat in int plane_front;
+#endif
+
+uniform lowp usampler3D voxel_id;
+uniform uint objectid;
+uniform float gap;
+
+{{uv_map_type}} uv_map;
+{{color_map_type}} color_map;
+{{color_type}} color;
+
+vec4 debug_color(uint id) {
+    return vec4(
+        float((id & uint(225)) >> uint(5)) / 5.0,
+        float((id & uint(25)) >> uint(3)) / 3.0,
+        float((id & uint(7)) >> uint(1)) / 3.0,
+        1.0
+    );
+}
+vec4 debug_color(int id) { return debug_color(uint(id)); }
+
+// unused but compilation requires it
+vec4 get_lrbt(Nothing uv_map, int id, int side) {
+    return vec4(0,0,1,1);
+}
+vec4 get_lrbt(sampler1D uv_map, int id, int side) {
+    return texelFetch(uv_map, id-1, 0);
+}
+vec4 get_lrbt(sampler2D uv_map, int id, int side) {
+    return texelFetch(uv_map, ivec2(id-1, side), 0);
+}
+
+vec4 get_color_from_texture(sampler2D color, int id) {
+    vec4 lrbt = get_lrbt(uv_map, id, o_side);
+    // compute uv normalized to voxel
+    // TODO: float precision causes this to wrap sometimes (e.g. 5.999..7.0002)
+    vec2 voxel_uv = mod(o_tex_uv, 1.0);
+    voxel_uv = mix(lrbt.xz, lrbt.yw, voxel_uv);
+    return texture(color, voxel_uv);
+}
+
+
+vec4 get_color(Nothing color, Nothing color_map, int id) {
+    return debug_color(id);
+}
+vec4 get_color(Nothing color, sampler1D color_map, int id) {
+    return texelFetch(color_map, id-1, 0);
+}
+vec4 get_color(sampler1D color, sampler1D color_map, int id) {
+    return texelFetch(color, id-1, 0);
+}
+vec4 get_color(sampler1D color, Nothing color_map, int id) {
+    return texelFetch(color, id-1, 0);
+}
+vec4 get_color(sampler2D color, sampler1D color_map, int id) {
+    return get_color_from_texture(color, id);
+}
+vec4 get_color(sampler2D color, Nothing color_map, int id) {
+    return get_color_from_texture(color, id);
+}
+
+
+void write2framebuffer(vec4 color, uvec2 id);
+
+#ifndef NO_SHADING
+vec3 illuminate(vec3 normal, vec3 base_color);
+#endif
+
+void main()
+{
+    vec2 voxel_uv = mod(o_tex_uv, 1.0);
+    if (voxel_uv.x < 0.5 * gap || voxel_uv.x > 1.0 - 0.5 * gap ||
+        voxel_uv.y < 0.5 * gap || voxel_uv.y > 1.0 - 0.5 * gap)
+        discard;
+
+    // grab voxel id
+    int id = int(texture(voxel_id, o_uvw).x);
+
+    // id is invisible so we simply discard
+    if (id == 0) {
+        discard;
+    }
+
+    // otherwise we draw. For now just some color...
+    vec4 voxel_color = get_color(color, color_map, id);
+
+#ifdef DEBUG_RENDER_ORDER
+    if (plane_dim != DEBUG_RENDER_ORDER)
+        discard;
+    voxel_color = vec4(
+        plane_front * plane_render_idx,
+        -plane_front * plane_render_idx,
+        0,
+        id == 0 ? 0.1 : 1.0
+    );
+    // voxel_color = vec4(o_normal, id == 0 ? 0.1 : 1.0);
+    // voxel_color = vec4(plane_front, 0, 0, 1.0);
+#endif
+
+#ifndef NO_SHADING
+    voxel_color.rgb = illuminate(o_normal, voxel_color.rgb);
+#endif
+
+    // TODO: index into 3d array
+    ivec3 size = ivec3(textureSize(voxel_id, 0).xyz);
+    ivec3 idx = ivec3(o_uvw * size);
+    int lin = 1 + idx.x + size.x * (idx.y + size.y * idx.z);
+
+    // draw
+    write2framebuffer(voxel_color, uvec2(objectid, lin));
+}
diff --git a/GLMakie/assets/shader/voxel.vert b/GLMakie/assets/shader/voxel.vert
@@ -0,0 +1,184 @@
+#version 330 core
+// {{GLSL_VERSION}}
+// {{GLSL_EXTENSIONS}}
+
+// debug FLAGS
+// #define DEBUG_RENDER_ORDER
+
+in vec2 vertices;
+
+flat out vec3 o_normal;
+out vec3 o_uvw;
+flat out int o_side;
+out vec2 o_tex_uv;
+
+#ifdef DEBUG_RENDER_ORDER
+flat out float plane_render_idx;
+flat out int plane_dim;
+flat out int plane_front;
+#endif
+
+out vec3 o_camdir;
+out vec3 o_world_pos;
+
+uniform mat4 model;
+uniform mat3 world_normalmatrix;
+uniform mat4 projectionview;
+uniform vec3 eyeposition;
+uniform vec3 view_direction;
+uniform lowp usampler3D voxel_id;
+uniform float depth_shift;
+uniform bool depthsorting;
+uniform float gap;
+
+const vec3 unit_vecs[3] = vec3[]( vec3(1, 0, 0), vec3(0, 1, 0), vec3(0, 0, 1) );
+const mat2x3 orientations[3] = mat2x3[](
+    mat2x3(0, 1, 0, 0, 0, 1), // xy -> _yz (x normal)
+    mat2x3(1, 0, 0, 0, 0, 1), // xy -> x_z (y normal)
+    mat2x3(1, 0, 0, 0, 1, 0)  // xy -> xy_ (z normal)
+);
+
+void main() {
+    /* How this works:
+    To simplify lets consider a 2d grid of pixel where the voxel surface would
+    be the square outline of around a data point x.
+        +---+---+---+
+        | x | x | x |
+        +---+---+---+
+        | x | x | x |
+        +---+---+---+
+        | x | x | x |
+        +---+---+---+
+    Naively we would draw 4 lines for each point x, coloring them based on the
+    data attached to x. This would result in 4 * N^2 lines with N^2 = number of
+    pixels. We can do much better though by drawing a line for each column and
+    row of pixels:
+      1 +---+---+---+
+        | x | x | x |
+      2 +---+---+---+
+        | x | x | x |
+      3 +---+---+---+
+        | x | x | x |
+      4 +---+---+---+
+        5   6   7   8
+    This results in 2 * (N+1) lines. We can adjust the color of the line by
+    sampling a Texture containing the information previously attached to vertices.
+
+    Generalized to 3D voxels, lines become planes and the texture becomes 3D.
+    We draw the planes through instancing. So first we will need to map the
+    instance id to a dimension (xy, xz or yz plane) and an offset (in z, y or
+    x direction respectively).
+    */
+
+    // TODO: might be better for transparent rendering to alternate xyz?
+    // But how would we do this for non-cubic chunks?
+
+    // Map instance id to dimension and index along dimension (0..N+1 or 0..2N)
+    ivec3 size = textureSize(voxel_id, 0);
+    int dim, id = gl_InstanceID, front = 1;
+    if (gap > 0.01) {
+        front = 1 - 2 * int(gl_InstanceID & 1);
+        if (id < 2 * size.z) {
+            dim = 2;
+            id = id;
+        } else if (id < 2 * (size.z + size.y)) {
+            dim = 1;
+            id = id - 2 * size.z;
+        } else { // if (id > 2 * (size.z + size.y)) {
+            dim = 0;
+            id = id - 2 * (size.z + size.y);
+        }
+    } else {
+        if (id < size.z + 1) {
+            dim = 2;
+            id = id;
+        } else if (id < size.z + size.y + 2) {
+            dim = 1;
+            id = id - (size.z + 1);
+        } else {
+            dim = 0;
+            id = id - (size.z + size.y + 2);
+        }
+    }
+
+#ifdef DEBUG_RENDER_ORDER
+    plane_render_idx = float(id) / float(size[dim]-1);
+    plane_dim = dim;
+    plane_front = front;
+#endif
+
+    // plane placement
+    // Figure out which plane to start with
+    vec3 normal = world_normalmatrix * unit_vecs[dim];
+    int dir = int(sign(dot(view_direction, normal))), start;
+    if (depthsorting) {
+        // TODO: depthsorted should start far away from viewer so every plane draws
+        start = int((0.5 + 0.5 * dir) * size[dim]);
+        dir *= -1;
+    } else {
+        // otherwise we should start at viewer and expand in view direction so
+        // that the depth test can quickly eliminate unnecessary fragments
+        // Note that model can have rotations and (uneven) scaling
+        vec4 origin = model * vec4(0, 0, 0, 1);
+        vec4 back   = model * vec4(size, 1);
+        float front_dist = dot(origin.xyz / origin.w, normal);
+        float back_dist  = dot(back.xyz / back.w, normal);
+        float cam_dist   = dot(eyeposition, normal);
+        float dist01 = (cam_dist - front_dist) / (back_dist - front_dist);
+
+        // index of voxel closest to (and in front of) the camera
+        start = clamp(int(float(size[dim]) * dist01), 0, size[dim]);
+    }
+
+    vec3 displacement;
+    if (gap > 0.01) {
+        // planes are doubled
+        // 2 * start + min(dir, 0)                  closest (camera facing) plane
+        // dir * id                                 iterate away from first plane
+        // (idx + 2 * size[dim]) % 2 * size[dim]    normalize to [0, 2size[dim])
+        int plane_idx = (2 * start + min(dir, 0) + dir * id + 2 * size[dim]) % (2 * size[dim]);
+        // (plane_idx + 1) / 2          map to idx 0, 1, 2, 3, 4 -> displacements 0, 1, 1, 2, 2, ...
+        // 0.5 * dir * gap * front      gap based offset from space filling placements
+        displacement = ((plane_idx + 1) / 2 + 0.5 * dir * front * gap) * unit_vecs[dim];
+    } else {
+        // similar to above but with N+1 indices around N voxels
+        displacement = ((start + dir * id + size[dim] + 1) % (size[dim] + 1)) * unit_vecs[dim];
+    }
+
+    // place plane vertices
+    vec3 plane_vertex = size * (orientations[dim] * vertices) + displacement;
+    vec4 world_pos = model * vec4(plane_vertex, 1.0f);
+    o_world_pos = world_pos.xyz;
+    gl_Position = projectionview * world_pos;
+    gl_Position.z += gl_Position.w * depth_shift;
+
+    // For each plane the normal is constant in
+    // +- normal = +- world_normalmatrix * unit_vecs[dim] direction. We just need
+    // the correct prefactor
+    o_camdir = eyeposition - world_pos.xyz / world_pos.w;
+    float normal_dir;
+    if (gap > 0.01) {
+        // With a gap we render front and back faces. Since the gap always takes
+        // away from a voxel the normal should go in the opposite direction.
+        normal_dir = -float(dir * front);
+    } else {
+        // Without a gap we skip back facing faces so every normal faces the camera
+        normal_dir = sign(dot(o_camdir, normal));
+    }
+    o_normal = normal_dir * normalize(normal);
+
+    // The quad we render here is directly between two slices of voxels in our
+    // chunk. Each `plane_vertex` of the quad is in a 0 .. size scale, so they
+    // can be mapped directly to texture indices. We just need to figure out if
+    // the quad is associated with the "previous" or "next" slice of voxels. We
+    // can derive that from the normal direction, as the normal always points
+    // away from the voxel center.
+    o_uvw = (plane_vertex - 0.5 * (1.0 - gap) * o_normal) / size;
+
+    // normal in: -x -y -z +x +y +z direction
+    o_side = dim + 3 * int(0.5 + 0.5 * normal_dir);
+
+    // map plane_vertex (-w/2 .. w/2 scale) back to 2d (scaled 0 .. w)
+    // if the normal is negative invert range (w .. 0)
+    o_tex_uv = transpose(orientations[dim]) * (vec3(-normal_dir, normal_dir, 1.0) * plane_vertex);
+}