prototype WGLMakie version

WGLMakie/assets/voxel.frag

@@ -0,0 +1,130 @@
+// debug FLAGS
+// #define DEBUG_RENDER_ORDER 2 // (0, 1, 2) - dimensions
+
+flat in vec3 o_normal;
+in vec3 o_uvw;
+flat in int o_side;
+in vec2 o_tex_uv;
+
+in vec3 o_camdir;
+
+#ifdef DEBUG_RENDER_ORDER
+flat in float plane_render_idx; // debug
+#endif
+
+// uniform isampler3D voxel_id;
+// uniform uint objectid;
+
+// {{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)); }
+
+vec4 get_color(bool color, bool color_map, bool uv_map, int id) {
+    return debug_color(id);
+}
+vec4 get_color(bool color, sampler2D color_map, bool uv_map, int id) {
+    return texelFetch(color_map, ivec2(id-1, 0), 0);
+}
+vec4 get_color(sampler2D color, sampler2D color_map, bool uv_map, int id) {
+    return texelFetch(color, ivec2(id-1, 0), 0);
+}
+vec4 get_color(sampler2D color, bool color_map, bool uv_map, int id) {
+    return texelFetch(color, ivec2(id-1, 0), 0);
+}
+vec4 get_color(sampler2D color, sampler2D color_map, sampler2D uv_map, int id) {
+    vec4 lrbt = texelFetch(uv_map, ivec2(id-1, o_side), 0);
+    // 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(sampler2D color, bool color_map, sampler2D uv_map, int id) {
+    vec4 lrbt = texelFetch(uv_map, ivec2(id-1, o_side), 0);
+    // 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);
+}
+
+// Smoothes out edge around 0 light intensity, see GLMakie
+float smooth_zero_max(float x) {
+    const float c = 0.00390625, xswap = 0.6406707120152759, yswap = 0.20508383900190955;
+    const float shift = 1.0 + xswap - yswap;
+    float pow8 = x + shift;
+    pow8 = pow8 * pow8; pow8 = pow8 * pow8; pow8 = pow8 * pow8;
+    return x < yswap ? c * pow8 : x;
+}
+
+vec3 blinnphong(vec3 N, vec3 V, vec3 L, vec3 color){
+    float diff_coeff = smooth_zero_max(dot(L, -N));
+
+    // specular coefficient
+    vec3 H = normalize(L + V);
+
+    float spec_coeff = pow(max(dot(H, -N), 0.0), get_shininess());
+    if (diff_coeff <= 0.0)
+        spec_coeff = 0.0;
+
+    // final lighting model
+    return get_light_color() * vec3(
+        get_diffuse() * diff_coeff * color +
+        get_specular() * spec_coeff
+    );
+}
+
+flat in uint frag_instance_id;
+vec4 pack_int(uint id, uint index) {
+    vec4 unpack;
+    unpack.x = float((id & uint(0xff00)) >> 8) / 255.0;
+    unpack.y = float((id & uint(0x00ff)) >> 0) / 255.0;
+    unpack.z = float((index & uint(0xff00)) >> 8) / 255.0;
+    unpack.w = float((index & uint(0x00ff)) >> 0) / 255.0;
+    return unpack;
+}
+void main()
+{
+    // 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, get_uv_map(), id);
+
+#ifdef DEBUG_RENDER_ORDER
+    if (mod(o_side, 3) != DEBUG_RENDER_ORDER)
+        discard;
+    voxel_color = vec4(plane_render_idx, 0, 0, id == 0 ? 0.01 : 1.0);
+#endif
+
+    if(get_shading()){
+        vec3 L = get_light_direction();
+        vec3 light = blinnphong(o_normal, normalize(o_camdir), L, voxel_color.rgb);
+        voxel_color.rgb = get_ambient() * voxel_color.rgb + light;
+    }
+
+    if (picking) {
+        uvec3 size = uvec3(textureSize(voxel_id, 0).xyz);
+        uvec3 idx = uvec3(o_uvw * vec3(size));
+        uint lin = uint(1) + idx.x + size.x * (idx.y + size.y * idx.z);
+        fragment_color = pack_int(object_id, lin);
+        return;
+    }
+
+    fragment_color = voxel_color;
+}

WGLMakie/assets/voxel.vert

@@ -0,0 +1,134 @@
+// 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;
+#endif
+
+out vec3 o_camdir;
+
+uniform mat4 projection, view;
+
+// uniform mat4 model;
+// uniform mat3 world_normalmatrix;
+// uniform vec3 eyeposition;
+// uniform vec3 view_direction;
+// uniform float depth_shift;
+// uniform bool depthsorting;
+
+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?
+    // How do we do this for non-cubic chunks?
+    ivec3 size = textureSize(voxel_id, 0);
+    int dim = 2, id = gl_InstanceID;
+    if (gl_InstanceID > size.z + size.y + 1) {
+        dim = 0;
+        id = gl_InstanceID - (size.z + size.y + 2);
+    } else if (gl_InstanceID > size.z) {
+        dim = 1;
+        id = gl_InstanceID - (size.z + 1);
+    }
+
+#ifdef DEBUG_RENDER_ORDER
+    plane_render_idx = float(id) / float(size[dim]-1);
+#endif
+
+    // plane placement
+    // Figure out which plane to start with
+    vec3 normal = get_normalmatrix() * unit_vecs[dim];
+    float dir = sign(dot(get_view_direction(), normal));
+    vec3 displacement;
+    if (depthsorting) {
+        // depthsorted should start far away from viewer so every plane draws
+        displacement = ((0.5 + 0.5 * dir) * float(size[dim]) - dir * float(id)) * unit_vecs[dim];
+    } else {
+        // no sorting should start at viewer and expand in view direction so
+        // that depth test can quickly eliminate unnecessary fragments
+        vec4 origin = get_model() * vec4(0, 0, 0, 1);
+        float dist = dot(get_eyeposition() - origin.xyz / origin.w, normal) / dot(normal, normal);
+        float start = clamp(dist, 0.0, float(size[dim]));
+        // this should work better with integer modulo...
+        displacement = mod(start + dir * float(id), float(size[dim]) + 0.001) * unit_vecs[dim];
+    }
+
+    // place plane vertices
+    vec3 voxel_pos = vec3(size) * (orientations[dim] * vertices) + displacement;
+    vec4 world_pos = get_model() * vec4(voxel_pos, 1.0f);
+    gl_Position = projection * view * world_pos;
+    gl_Position.z += gl_Position.w * get_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:
+    o_camdir = get_eyeposition() - world_pos.xyz / world_pos.w;
+    float normal_dir = sign(dot(o_camdir, normal));
+    o_normal = normalize(normal_dir * normal);
+
+    // 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.
+    o_uvw = (voxel_pos - 0.5 * o_normal) / vec3(size);
+
+    // normal in: -x -y -z +x +y +z direction
+    o_side = dim + 3 * int(0.5 + 0.5 * normal_dir);
+
+    // map voxel_pos (-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]) * (normal_dir * voxel_pos);
+}

WGLMakie/src/WGLMakie.jl

@@ -42,6 +42,7 @@ include("lines.jl")
 include("meshes.jl")
 include("imagelike.jl")
 include("picking.jl")
+include("voxel.jl")
 
 const LAST_INLINE = Base.RefValue{Union{Automatic, Bool}}(Makie.automatic)
 

WGLMakie/src/serialization.jl

@@ -81,6 +81,7 @@ three_format(::Type{<:RGBA}) = "RGBAFormat"
 three_type(::Type{Float16}) = "FloatType"
 three_type(::Type{Float32}) = "FloatType"
 three_type(::Type{N0f8}) = "UnsignedByteType"
+three_type(::Type{UInt8}) = "UnsignedByteType"
 
 function three_filter(sym::Symbol)
     sym === :linear && return "LinearFilter"

WGLMakie/src/voxel.jl

@@ -0,0 +1,96 @@
+function create_shader(scene::Scene, plot::Makie.Voxel)
+
+    uniform_dict = Dict{Symbol, Any}(
+        :voxel_id => Sampler(plot.converted[end], minfilter = :nearest),
+        # for plane sorting
+        :depthsorting => plot.depthsorting,
+        :eyeposition => Vec3f(1),
+        :view_direction => camera(scene).view_direction,
+        # lighting
+        :diffuse => lift(x -> convert_attribute(x, Key{:diffuse}()), plot, plot.diffuse),
+        :specular => lift(x -> convert_attribute(x, Key{:specular}()), plot, plot.specular),
+        :shininess => lift(x -> convert_attribute(x, Key{:shininess}()), plot, plot.shininess),
+        :depth_shift => get(plot, :depth_shift, Observable(0.0f0)),
+        :light_direction => Vec3f(1),
+        :light_color => Vec3f(1),
+        :ambient => Vec3f(1),
+        # picking
+        :picking => false,
+        :object_id => UInt32(0),
+        # other
+        :normalmatrix => map(plot.model) do m
+            # should be fine to ignore placement matrix here because
+            # translation is ignored and scale shouldn't matter
+            i = Vec(1, 2, 3)
+            return transpose(inv(m[i, i]))
+        end,
+        :shading => to_value(get(plot, :shading, NoShading)) != NoShading,
+    )
+
+    # TODO: localized update
+    # buffer = Vector{UInt8}(undef, 1)
+    on(plot, plot._local_update) do (is, js, ks)
+        # required_length = length(is) * length(js) * length(ks)
+        # if length(buffer) < required_length
+        #     resize!(buffer, required_length)
+        # end
+        # idx = 1
+        # for k in ks, j in js, i in is
+        #     buffer[idx] = plot.converted[end].val[i, j, k]
+        #     idx += 1
+        # end
+        # GLAbstraction.texsubimage(tex, buffer, is, js, ks)
+        notify(plot.converted[end])
+        return
+    end
+
+    # adjust model matrix with placement matrix
+    uniform_dict[:model] = map(
+            plot, plot.converted...,  plot.model
+        ) do xs, ys, zs, chunk, model
+        mini = minimum.((xs, ys, zs))
+        width = maximum.((xs, ys, zs)) .- mini
+        return model *
+            Makie.scalematrix(Vec3f(width ./ size(chunk))) *
+            Makie.translationmatrix(Vec3f(mini))
+    end
+
+    maybe_color_mapping = plot.calculated_colors[]
+    uv_map = plot.uvmap
+    if maybe_color_mapping isa Makie.ColorMapping
+        uniform_dict[:color_map] = Sampler(maybe_color_mapping.colormap, minfilter = :nearest)
+        uniform_dict[:uv_map] = false
+        uniform_dict[:color] = false
+    elseif !isnothing(to_value(uv_map))
+        uniform_dict[:color_map] = false
+        # WebGL doesn't have sampler1D so we need to pad id -> uv mappings to
+        # (id, side) -> uv mappings
+        wgl_uv_map = map(plot, uv_map) do uv_map
+            if uv_map isa Vector
+                new_map = Matrix{Vec4f}(undef, length(uv_map), 6)
+                for col in 1:6
+                    new_map[:, col] .= uv_map
+                end
+                return new_map
+            else
+                return uv_map
+            end
+        end
+        uniform_dict[:uv_map] = Sampler(wgl_uv_map, minfilter = :nearest)
+        interp = to_value(plot.interpolate) ? :linear : :nearest
+        uniform_dict[:color] = Sampler(maybe_color_mapping, minfilter = interp)
+    else
+        uniform_dict[:color_map] = false
+        uniform_dict[:uv_map] = false
+        uniform_dict[:color] = Sampler(maybe_color_mapping, minfilter = :nearest)
+    end
+
+    # TODO: this is a waste
+    N_instances = sum(size(plot.converted[end][])) + 3
+    dummy_data = [0f0 for _ in 1:N_instances]
+
+    instance = uv_mesh(Rect2(0f0, 0f0, 1f0, 1f0))
+
+    return InstancedProgram(WebGL(), lasset("voxel.vert"), lasset("voxel.frag"),
+                        instance, VertexArray(dummy = dummy_data), uniform_dict)
+end