create voxel rendering prototype

GLMakie/assets/shader/voxel.frag

@@ -0,0 +1,38 @@
+#version 330 core
+// {{GLSL VERSION}}
+// {{GLSL_EXTENSIONS}}
+
+in vec3 f_normal;
+in vec3 f_uvw;
+
+uniform isampler3D voxel_id;
+uniform uint objectid;
+
+void write2framebuffer(vec4 color, uvec2 id);
+
+void main()
+{
+    // grab voxel id
+    uint id = uint(texture(voxel_id, f_uvw).x);
+
+    // id is invisible so we simply discard
+    if (id == uint(0)) {
+        discard;
+    }
+
+    // otherwise we draw. For now just some color...
+    vec4 color = 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
+    );
+
+    // TODO: index into 3d array
+    uvec3 size = uvec3(textureSize(voxel_id, 0).xyz);
+    uvec3 idx = uvec3(f_uvw * size);
+    uint lin = uint(1) + idx.x + size.x * (idx.y + size.y * idx.z);
+
+    // draw
+    write2framebuffer(color, uvec2(objectid, lin));
+}

GLMakie/assets/shader/voxel.vert

@@ -0,0 +1,106 @@
+#version 330 core
+// {{GLSL_VERSION}}
+// {{GLSL_EXTENSIONS}}
+
+in vec2 vertices;
+
+flat out vec3 f_normal;
+out vec3 f_uvw;
+
+uniform mat4 model;
+uniform mat3 world_normalmatrix;
+uniform mat4 projectionview;
+uniform vec3 eyeposition;
+uniform isampler3D voxel_id;
+uniform float depth_shift;
+
+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).
+    Note that the render order of planes can make a significant impact on
+    performance. It may be worth it to adjust this based on eyeposition.
+
+    For now we alternate x, y, z planes and start from the center.
+    */
+    ivec3 size = textureSize(voxel_id, 0);
+    int dim = 0, id = gl_InstanceID;
+    if (gl_InstanceID > size.x + size.y + 1) {
+        dim = 2;
+        id = gl_InstanceID - (size.x + size.y + 2);
+    } else if (gl_InstanceID > size.x) {
+        dim = 1;
+        id = gl_InstanceID - (size.x + 1);
+    }
+
+    // start rendering from center
+    ivec3 offset = size / 2;
+    id = id <= (size[dim] >> 1) ? -id : id - (size[dim] >> 1);
+
+    // plane placement
+    vec3 displacement = id * unit_vecs[dim];
+    vec3 voxel_pos = size * (orientations[dim] * vertices) + displacement;
+    vec4 world_pos = model * vec4(voxel_pos, 1.0f);
+    gl_Position = projectionview * world_pos;
+    gl_Position.z += gl_Position.w * depth_shift;
+
+    // For each plane the normal is constant and its direction is given by the
+    // `displacement` direction, i.e. `n = unit_vecs[dim]`. We just need to derive
+    // whether it's +n or -n.
+    // If we assume the viewer to be outside of a voxel, the normal direction
+    // should always be facing them. Thus:
+    vec3 n = world_normalmatrix * unit_vecs[dim];
+    vec3 camdir = eyeposition - world_pos.xyz / world_pos.w;
+    f_normal = normalize(sign(dot(camdir, n)) * n);
+
+    // The texture coordinate can also be derived. `voxel_pos` effectively gives
+    // an integer index into the chunk, shifted to be centered. We can convert
+    // this to a float index into the voxel_id texture by normalizing.
+    // The minor ceveat here is that because planes are drawn between voxels we
+    // would be sampling between voxels like this. To fix this we want to shift
+    // the uvw coordinate to the relevant voxel center, which we can do using the
+    // normal direction.
+    // Here we want to shift in -normal direction to get a front face. Consider
+    // this example with 1, 2 solid, 0 air and v the viewer:
+    // | 1 | 2 | 0 |   v
+    // If we shift in +normal direction (towards viewer) the planes would sample
+    // from the id closer to the viewer, drawing a backface.
+    f_uvw = (voxel_pos - 0.5 * f_normal) / size + 0.5;
+
+    // TODO:
+    // - verify no missing faces (we are missing 3 planes I think, but that should be possible...)
+    // - verify idx
+}

GLMakie/src/GLAbstraction/GLShader.jl

@@ -124,7 +124,7 @@ function compile_shader(source::ShaderSource, template_src::String)
     glShaderSource(shaderid, template_src)
     glCompileShader(shaderid)
     if !GLAbstraction.iscompiled(shaderid)
-        GLAbstraction.print_with_lines(String(source))
+        GLAbstraction.print_with_lines(String(source.source))
         @warn("shader $(name) didn't compile. \n$(GLAbstraction.getinfolog(shaderid))")
     end
     return Shader(name, Vector{UInt8}(template_src), source.typ, shaderid)

GLMakie/src/drawing_primitives.jl

@@ -815,3 +815,11 @@ function draw_atomic(screen::Screen, scene::Scene, plot::Volume)
         end
     end
 end
+
+function draw_atomic(screen::Screen, scene::Scene, plot::Voxel)
+    return cached_robj!(screen, scene, plot) do gl_attributes
+        @assert to_value(plot[1]) isa Array{UInt8, 3}
+        tex = Texture(plot[1], minfilter = :nearest)
+        return draw_voxels(screen, tex, gl_attributes)
+    end
+end

GLMakie/src/gl_backend.jl

@@ -69,6 +69,7 @@ include("glshaders/image_like.jl")
 include("glshaders/mesh.jl")
 include("glshaders/particles.jl")
 include("glshaders/surface.jl")
+include("glshaders/voxel.jl")
 
 include("picking.jl")
 include("rendering.jl")

GLMakie/src/glshaders/voxel.jl

@@ -0,0 +1,39 @@
+@nospecialize
+function draw_voxels(screen, main::VolumeTypes, data::Dict)
+    geom = Rect2f(Point2f(-0.5), Vec2f(1.0))
+    to_opengl_mesh!(data, const_lift(GeometryBasics.triangle_mesh, geom))
+    shading = pop!(data, :shading, FastShading)
+
+    # Not applicable
+    for key in (
+            :color, :lowclip,
+            :color_map, :nan_color,
+            :color_norm, :highclip,
+            :px_per_unit, :pixel_space,
+            :alpha, :intensity, :backlight
+        )
+        pop!(data, key)
+    end
+    @gen_defaults! data begin
+        voxel_id = main => Texture
+        instances = const_lift(t -> sum(size(t)) + 3, voxel_id)
+        model = Mat4f(I)
+        # modelinv = const_lift(inv, model)
+        transparency = false
+        shader = GLVisualizeShader(
+            screen,
+            "voxel.vert",
+            "fragment_output.frag", "voxel.frag", # "lighting.frag", "volume.frag",
+            view = Dict(
+                # "shading" => light_calc(shading),
+                # "MAX_LIGHTS" => "#define MAX_LIGHTS $(screen.config.max_lights)",
+                # "MAX_LIGHT_PARAMETERS" => "#define MAX_LIGHT_PARAMETERS $(screen.config.max_light_parameters)",
+                "buffers" => output_buffers(screen, to_value(transparency)),
+                "buffer_writes" => output_buffer_writes(screen, to_value(transparency))
+            )
+        )
+    end
+
+    return assemble_shader(data)
+end
+@specialize

MakieCore/src/basic_plots.jl

@@ -498,6 +498,14 @@ $(Base.Docs.doc(MakieCore.generic_plot_attributes!))
     return colormap_attributes!(attr, theme(scene, :colormap))
 end
 
+@recipe(Voxel, chunk) do scene
+    attr = Attributes()
+    shading_attributes!(attr)
+    generic_plot_attributes!(attr)
+    return colormap_attributes!(attr, theme(scene, :colormap))
+end
+
+
 """
     poly(vertices, indices; kwargs...)
     poly(points; kwargs...)

src/Makie.jl

@@ -89,8 +89,8 @@ using MakieCore: not_implemented_for
 import MakieCore: plot, plot!, theme, plotfunc, plottype, merge_attributes!, calculated_attributes!,
                   get_attribute, plotsym, plotkey, attributes, used_attributes
 import MakieCore: create_axis_like, create_axis_like!, figurelike_return, figurelike_return!
-import MakieCore: arrows, heatmap, image, lines, linesegments, mesh, meshscatter, poly, scatter, surface, text, volume
-import MakieCore: arrows!, heatmap!, image!, lines!, linesegments!, mesh!, meshscatter!, poly!, scatter!, surface!, text!, volume!
+import MakieCore: arrows, heatmap, image, lines, linesegments, mesh, meshscatter, poly, scatter, surface, text, volume, voxel
+import MakieCore: arrows!, heatmap!, image!, lines!, linesegments!, mesh!, meshscatter!, poly!, scatter!, surface!, text!, volume!, voxel!
 import MakieCore: convert_arguments, convert_attribute, default_theme, conversion_trait
 
 export @L_str, @colorant_str
@@ -354,9 +354,9 @@ include("basic_recipes/text.jl")
 include("basic_recipes/raincloud.jl")
 include("deprecated.jl")
 
-export Arrows  , Heatmap  , Image  , Lines  , LineSegments  , Mesh  , MeshScatter  , Poly  , Scatter  , Surface  , Text  , Volume  , Wireframe
-export arrows  , heatmap  , image  , lines  , linesegments  , mesh  , meshscatter  , poly  , scatter  , surface  , text  , volume  , wireframe
-export arrows! , heatmap! , image! , lines! , linesegments! , mesh! , meshscatter! , poly! , scatter! , surface! , text! , volume! , wireframe!
+export Arrows  , Heatmap  , Image  , Lines  , LineSegments  , Mesh  , MeshScatter  , Poly  , Scatter  , Surface  , Text  , Volume  , Wireframe,  Voxel
+export arrows  , heatmap  , image  , lines  , linesegments  , mesh  , meshscatter  , poly  , scatter  , surface  , text  , volume  , wireframe,  voxel
+export arrows! , heatmap! , image! , lines! , linesegments! , mesh! , meshscatter! , poly! , scatter! , surface! , text! , volume! , wireframe!, voxel!
 
 export AmbientLight, PointLight, DirectionalLight, SpotLight, EnvironmentLight, RectLight, SSAO
 

src/interfaces.jl

@@ -54,6 +54,11 @@ function calculated_attributes!(::Type{<:Volume}, plot)
     return
 end
 
+function calculated_attributes!(::Type{<:Voxel}, plot)
+    color_and_colormap!(plot, plot[1])
+    return
+end
+
 function calculated_attributes!(::Type{<:Text}, plot)
     color_and_colormap!(plot)
     return
@@ -101,7 +106,7 @@ end
 
 const atomic_functions = (
     text, meshscatter, scatter, mesh, linesegments,
-    lines, surface, volume, heatmap, image
+    lines, surface, volume, heatmap, image, voxel
 )
 const Atomic{Arg} = Union{map(x-> Plot{x, Arg}, atomic_functions)...}
 

src/layouting/data_limits.jl

@@ -57,6 +57,15 @@ function data_limits(text::Text)
     return data_limits(text.plots[1])
 end
 
+function point_iterator(plot::Voxel)
+    s = size(plot[1][])
+    r = Rect3f(Point3f(-0.5 .* s), Vec3f(s))
+    return map(corners(r)) do p
+        p4d = plot.model[] * Point4f(p[1], p[2], p[3], 1.0)
+        return Point3f(p4d) / p4d[4]
+    end
+end
+
 point_iterator(mesh::GeometryBasics.Mesh) = decompose(Point, mesh)
 
 function point_iterator(list::AbstractVector)