diff --git a/CHANGELOG.md b/CHANGELOG.md
index d996546ab8e..4b87b83f46e 100644
--- 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
index 09499f2aa9b..f880494ce8f 100644
--- 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
index 30d050497bb..c066b345543 100644
--- 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
new file mode 100644
index 00000000000..cc5caf52c34
--- /dev/null
+++ 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));
+}
\ No newline at end of file
diff --git a/GLMakie/assets/shader/voxel.vert b/GLMakie/assets/shader/voxel.vert
new file mode 100644
index 00000000000..fa1ace9f348
--- /dev/null
+++ 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);
+}
\ No newline at end of file
diff --git a/GLMakie/src/GLAbstraction/GLTexture.jl b/GLMakie/src/GLAbstraction/GLTexture.jl
index 028c23db012..39d3cf3bb60 100644
--- a/GLMakie/src/GLAbstraction/GLTexture.jl
+++ b/GLMakie/src/GLAbstraction/GLTexture.jl
@@ -381,21 +381,25 @@ function gpu_resize!(t::Texture{T, ND}, newdims::NTuple{ND, Int}) where {T, ND}
     return t
 end
 
-texsubimage(t::Texture{T, 1}, newvalue::Array{T, 1}, xrange::UnitRange, level=0) where {T} = glTexSubImage1D(
-    t.texturetype, level, first(xrange)-1, length(xrange), t.format, t.pixeltype, newvalue
-)
-function texsubimage(t::Texture{T, 2}, newvalue::Array{T, 2}, xrange::UnitRange, yrange::UnitRange, level=0) where T
+function texsubimage(t::Texture{T, 1}, newvalue::Array{T}, xrange::UnitRange, level=0) where {T}
+    glTexSubImage1D(
+        t.texturetype, level, first(xrange)-1, length(xrange), t.format, t.pixeltype, newvalue
+    )
+end
+function texsubimage(t::Texture{T, 2}, newvalue::Array{T}, xrange::UnitRange, yrange::UnitRange, level=0) where T
     glTexSubImage2D(
         t.texturetype, level,
         first(xrange)-1, first(yrange)-1, length(xrange), length(yrange),
         t.format, t.pixeltype, newvalue
     )
 end
-texsubimage(t::Texture{T, 3}, newvalue::Array{T, 3}, xrange::UnitRange, yrange::UnitRange, zrange::UnitRange, level=0) where {T} = glTexSubImage3D(
-    t.texturetype, level,
-    first(xrange)-1, first(yrange)-1, first(zrange)-1, length(xrange), length(yrange), length(zrange),
-    t.format, t.pixeltype, newvalue
-)
+function texsubimage(t::Texture{T, 3}, newvalue::Array{T}, xrange::UnitRange, yrange::UnitRange, zrange::UnitRange, level=0) where {T}
+    glTexSubImage3D(
+        t.texturetype, level,
+        first(xrange)-1, first(yrange)-1, first(zrange)-1, length(xrange), length(yrange), length(zrange),
+        t.format, t.pixeltype, newvalue
+    )
+end
 
 Base.iterate(t::TextureBuffer{T}) where {T} = iterate(t.buffer)
 function Base.iterate(t::TextureBuffer{T}, state::Tuple{Ptr{T}, Int}) where T
diff --git a/GLMakie/src/GLMakie.jl b/GLMakie/src/GLMakie.jl
index 0c6244f24a5..49296d13556 100644
--- a/GLMakie/src/GLMakie.jl
+++ b/GLMakie/src/GLMakie.jl
@@ -57,9 +57,9 @@ end
 const GL_ASSET_DIR = RelocatableFolders.@path joinpath(@__DIR__, "..", "assets")
 const SHADER_DIR = RelocatableFolders.@path joinpath(GL_ASSET_DIR, "shader")
 const LOADED_SHADERS = Dict{String, Tuple{Float64, ShaderSource}}()
+
 function loadshader(name)
-    # Turns out, joinpath is so slow, that it actually makes sense
-    # To memoize it :-O
+    # Turns out, loading shaders is so slow, that it actually makes sense to memoize it :-O
     # when creating 1000 plots with the PlotSpec API, timing drop from 1.5s to 1s just from this change:
     # Note that we need to check if the file is still valid to enable hot reloading of shaders
     path = joinpath(SHADER_DIR, name)
diff --git a/GLMakie/src/drawing_primitives.jl b/GLMakie/src/drawing_primitives.jl
index 14c00d8479e..63760bf35f1 100644
--- a/GLMakie/src/drawing_primitives.jl
+++ b/GLMakie/src/drawing_primitives.jl
@@ -798,3 +798,74 @@ function draw_atomic(screen::Screen, scene::Scene, plot::Volume)
         end
     end
 end
+
+function draw_atomic(screen::Screen, scene::Scene, plot::Voxels)
+    return cached_robj!(screen, scene, plot) do gl_attributes
+        @assert to_value(plot.converted[end]) isa Array{UInt8, 3}
+
+        # voxel ids
+        tex = Texture(plot.converted[end], minfilter = :nearest)
+
+        # local update
+        buffer = Vector{UInt8}(undef, 1)
+        on(plot, pop!(gl_attributes, :_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)
+            return
+        end
+
+        # adjust model matrix according to x/y/z limits
+        gl_attributes[:model] = map(
+                plot, plot.converted...,  pop!(gl_attributes, :model)
+            ) do xs, ys, zs, chunk, model
+            mini = minimum.((xs, ys, zs))
+            width = maximum.((xs, ys, zs)) .- mini
+            return model *
+                Makie.translationmatrix(Vec3f(mini)) *
+                Makie.scalematrix(Vec3f(width ./ size(chunk)))
+        end
+
+        # color attribute adjustments
+        # TODO:
+        pop!(gl_attributes, :lowclip, nothing)
+        pop!(gl_attributes, :highclip, nothing)
+        # Invalid:
+        pop!(gl_attributes, :nan_color, nothing)
+        pop!(gl_attributes, :alpha, nothing) # Why is this even here?
+        pop!(gl_attributes, :intensity, nothing)
+        pop!(gl_attributes, :color_norm, nothing)
+        # cleanup
+        pop!(gl_attributes, :_limits)
+        pop!(gl_attributes, :is_air)
+
+        # make sure these exist
+        get!(gl_attributes, :color, nothing)
+        get!(gl_attributes, :color_map, nothing)
+
+        # process texture mapping
+        uv_map = pop!(gl_attributes, :uvmap)
+        if !isnothing(to_value(uv_map))
+            gl_attributes[:uv_map] = Texture(uv_map, minfilter = :nearest)
+
+            interp = to_value(pop!(gl_attributes, :interpolate))
+            interp = interp ? :linear : :nearest
+            color = gl_attributes[:color]
+            gl_attributes[:color] = Texture(color, minfilter = interp)
+        elseif !isnothing(to_value(gl_attributes[:color]))
+            gl_attributes[:color] = Texture(gl_attributes[:color], minfilter = :nearest)
+        end
+
+        # for depthsorting
+        gl_attributes[:view_direction] = camera(scene).view_direction
+
+        return draw_voxels(screen, tex, gl_attributes)
+    end
+end
diff --git a/GLMakie/src/gl_backend.jl b/GLMakie/src/gl_backend.jl
index b8ed7a35dfd..ee4bffe3402 100644
--- a/GLMakie/src/gl_backend.jl
+++ b/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")
diff --git a/GLMakie/src/glshaders/voxel.jl b/GLMakie/src/glshaders/voxel.jl
new file mode 100644
index 00000000000..f0a3c679405
--- /dev/null
+++ b/GLMakie/src/glshaders/voxel.jl
@@ -0,0 +1,35 @@
+@nospecialize
+function draw_voxels(screen, main::VolumeTypes, data::Dict)
+    geom = Rect2f(Point2f(0), Vec2f(1.0))
+    to_opengl_mesh!(data, const_lift(GeometryBasics.triangle_mesh, geom))
+    shading = pop!(data, :shading, FastShading)
+    @gen_defaults! data begin
+        voxel_id = main => Texture
+        gap = 0f0
+        instances = const_lift(gap, voxel_id) do gap, chunk
+            N = sum(size(chunk))
+            ifelse(gap > 0.01, 2 * N, N + 3)
+        end
+        model = Mat4f(I)
+        transparency = false
+        backlight = 0f0
+        color = nothing => Texture
+        color_map = nothing => Texture
+        uv_map = nothing => Texture
+        shader = GLVisualizeShader(
+            screen,
+            "voxel.vert",
+            "fragment_output.frag", "voxel.frag", "lighting.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
\ No newline at end of file
diff --git a/MakieCore/src/basic_plots.jl b/MakieCore/src/basic_plots.jl
index 113575f6fef..a6b29650696 100644
--- a/MakieCore/src/basic_plots.jl
+++ b/MakieCore/src/basic_plots.jl
@@ -479,6 +479,73 @@ function deprecated_attributes(::Type{<:Text})
     )
 end
 
+"""
+    voxels(x, y, z, chunk::Array{<:Real, 3})
+    voxels(chunk::Array{<:Real, 3})
+
+Plots a chunk of voxels centered at 0. Optionally the placement and scaling of
+the chunk can be given as range-like x, y and z. (Only the extrema are
+considered here. Voxels are always uniformly sized.)
+
+Internally voxels are represented as 8 bit unsigned integer, with `0x00` always
+being an invisible "air" voxel. Passing a chunk with matching type will directly
+set those values. Note that color handling is specialized for the internal
+representation and may behave a bit differently than usual.
+
+## Attributes
+
+### Specific to `Voxel`
+
+- `color = nothing` sets colors per voxel id, skipping `0x00`. This means that a voxel with id 1 will grab `plot.colors[1]`
+  and so on up to id 255. This can also be set to a Matrix of colors, i.e. an image for texture mapping.
+- `is_air = x -> isnothing(x) || ismissing(x) || isnan(x)` controls which values in the input data are mapped to invisible
+  (air) voxels.
+- `depthsorting = false` controls the render order of voxels. If set to `false` voxels close to the viewer are rendered
+  first which should reduce overdraw and yield better performance. If set to `true` voxels are rendered back to front
+  enabling correct order for transparent voxels.
+- `uvmap = nothing` defines a map from voxel ids (and optionally sides) to uv coordinates. These uv coordinates are then
+  used to sample a 2D texture passed through `color` for texture mapping.
+- `interpolate = false` controls whether the texture map is sampled with interpolation (i.e. smoothly) or not (i.e. pixelated).
+- `gap = 0.0` sets the gap between adjacent voxels in units of the voxel size. This needs to be larger than 0.01 to take effect.
+- `_limits`: Internal attribute for keeping track of `extrema(chunk)`.
+- `_local_update`: Internal attribute for communicating updates to the backend.
+
+$(Base.Docs.doc(shading_attributes!))
+
+### Color attributes
+
+- `colormap::Union{Symbol, Vector{<:Colorant}} = :viridis` sets the colormap that voxels sample to get their color.
+  Internally the colormap is always represented by 253 colors, matching voxel ids 2..254. Ids 1 and 255 are reserved
+  for `lowclip` and `highclip`. `PlotUtils.cgrad(...)`, `Makie.Reverse(any_colormap)` can be used as well, or any
+  symbol from ColorBrewer or PlotUtils. To see all available color gradients, you can call `Makie.available_gradients()`.
+- `colorscale::Function = identity` color transform function. Can be any function, but only works well together with
+  `Colorbar` for `identity`, `log`, `log2`, `log10`, `sqrt`, `logit`, `Makie.pseudolog10` and `Makie.Symlog10`.
+- `colorrange::Tuple{<:Real, <:Real}` sets the values representing the start and end points of `colormap`.
+  Internally this effects the voxel ids generated from the input data so that the colormap resolution remains at 253 colors.
+- `lowclip::Union{Nothing, Symbol, <:Colorant} = nothing` sets a color for any value below the colorrange. (voxel id = 1)
+- `highclip::Union{Nothing, Symbol, <:Colorant} = nothing` sets a color for any value above the colorrange. (voxel id = 255)
+- `alpha = 1.0` sets the alpha value of the colormap or color attribute. Multiple alphas like in `plot(alpha=0.2, color=(:red, 0.5)` will get multiplied.
+
+$(Base.Docs.doc(MakieCore.generic_plot_attributes!))
+"""
+@recipe(Voxels, chunk) do scene
+    attr = Attributes(
+        color = nothing,
+        is_air = x -> isnothing(x) || ismissing(x) || isnan(x),
+        uvmap = nothing,
+        interpolate = false, # texture
+        depthsorting = false, # false = fast, true = correct transparency
+        gap = 0.0,
+        # Internal:
+        _limits = (0.0, 1.0),
+        _local_update = (0:0, 0:0, 0:0),
+    )
+    shading_attributes!(attr)
+    generic_plot_attributes!(attr)
+    return colormap_attributes!(attr, theme(scene, :colormap))
+end
+
+
 """
     poly(vertices, indices; kwargs...)
     poly(points; kwargs...)
diff --git a/RPRMakie/src/meshes.jl b/RPRMakie/src/meshes.jl
index 43bb3862d17..e2d43bf074c 100644
--- a/RPRMakie/src/meshes.jl
+++ b/RPRMakie/src/meshes.jl
@@ -129,6 +129,40 @@ function to_rpr_object(context, matsys, scene, plot::Makie.MeshScatter)
 end
 
 
+function to_rpr_object(context, matsys, scene, plot::Makie.Voxels)
+    # Potentially per instance attributes
+    positions = Makie.voxel_positions(plot)
+    m_mesh = normal_mesh(Rect3f(Point3f(-0.5), Vec3f(1)))
+    marker = RPR.Shape(context, m_mesh)
+    instances = [marker]
+    n_instances = length(positions)
+    RPR.rprShapeSetObjectID(marker, 0)
+    material = extract_material(matsys, plot)
+    set!(marker, material)
+    for i in 1:(n_instances-1)
+        inst = RPR.Shape(context, marker)
+        RPR.rprShapeSetObjectID(inst, i)
+        push!(instances, inst)
+    end
+
+    color_from_num = Makie.voxel_colors(plot)
+    object_id = RPR.InputLookupMaterial(matsys)
+    object_id.value = RPR.RPR_MATERIAL_NODE_LOOKUP_OBJECT_ID
+    uv = object_id * Vec3f(0, 1/n_instances, 0)
+    tex = RPR.Texture(matsys, collect(color_from_num'); uv = uv)
+    material.color = tex
+
+    scales = Iterators.repeated(Makie.voxel_size(plot), n_instances)
+
+    for (instance, position, scale) in zip(instances, positions, scales)
+        mat = Makie.transformationmatrix(position, scale)
+        transform!(instance, mat)
+    end
+
+    return instances
+end
+
+
 function to_rpr_object(context, matsys, scene, plot::Makie.Surface)
     x = plot[1]
     y = plot[2]
diff --git a/ReferenceTests/src/tests/primitives.jl b/ReferenceTests/src/tests/primitives.jl
index b369406cd58..a161de9cbb3 100644
--- a/ReferenceTests/src/tests/primitives.jl
+++ b/ReferenceTests/src/tests/primitives.jl
@@ -523,6 +523,83 @@ end
     fig
 end
 
+@reference_test "Voxel - texture mapping" begin
+    texture = let
+        w = RGBf(1, 1, 1)
+        r = RGBf(1, 0, 0)
+        g = RGBf(0, 1, 0)
+        b = RGBf(0, 0, 1)
+        o = RGBf(1, 1, 0)
+        c = RGBf(0, 1, 1)
+        m = RGBf(1, 0, 1)
+        k = RGBf(0, 0, 0)
+        [r w g w b w k w;
+         w w w w w w w w;
+         r k g k b k w k;
+         k k k k k k k k]
+    end
+
+    # Use same uvs/texture-sections for every side of one voxel id
+    flat_uv_map = [Vec4f(x, x + 1 / 2, y, y + 1 / 4)
+                   for x in range(0.0, 1.0; length=3)[1:(end - 1)]
+                   for y in range(0.0, 1.0; length=5)[1:(end - 1)]]
+
+    flat_voxels = UInt8[1, 0, 3,
+                        0, 0, 0,
+                        2, 0, 4,
+                        0, 0, 0,
+                        0, 0, 0,
+                        0, 0, 0,
+                        5, 0, 7,
+                        0, 0, 0,
+                        6, 0, 8]
+
+    # Reshape the flat vector into a 3D array of dimensions 3x3x3.
+    voxels_3d = reshape(flat_voxels, (3, 3, 3))
+
+    fig = Figure(; size=(800, 400))
+    a1 = LScene(fig[1, 1]; show_axis=false)
+    p1 = voxels!(a1, voxels_3d; uvmap=flat_uv_map, color=texture)
+
+    # Use red for x, green for y, blue for z
+    sided_uv_map = Matrix{Vec4f}(undef, 1, 6)
+    sided_uv_map[1, 1:3] .= flat_uv_map[1:3]
+    sided_uv_map[1, 4:6] .= flat_uv_map[5:7]
+
+    sided_voxels = reshape(UInt8[1], 1, 1, 1)
+    a2 = LScene(fig[1, 2]; show_axis=false)
+    p2 = voxels!(a2, sided_voxels; uvmap=sided_uv_map, color=texture)
+
+    fig
+end
+
+@reference_test "Voxel - colors and colormap" begin
+    # test direct mapping of ids to colors & upsampling of vector colormap
+    fig = Figure(size = (800, 400))
+    chunk = reshape(UInt8[1, 0, 2, 0, 0, 0, 3, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 6, 0, 0, 0, 7, 0, 8],
+                    (3, 3, 3))
+
+    cs = [:white, :red, :green, :blue, :black, :orange, :cyan, :magenta]
+    voxels(fig[1, 1], chunk, color = cs, axis=(show_axis = false,))
+    a, p = voxels(fig[1, 2], Float32.(chunk), colormap = [:red, :blue], is_air = x -> x == 0.0, axis=(show_axis = false,))
+    fig
+end
+
+@reference_test "Voxel - lowclip and highclip" begin
+    # test that lowclip and highclip are visible for values just outside the colorrange
+    fig = Figure(size = (800, 400))
+    chunk = reshape(collect(1:900), 30, 30, 1)
+    a1, _ = voxels(fig[1,1], chunk, lowclip = :red, highclip = :red, colorrange = (1.0, 900.0), shading = NoShading)
+    a2, _ = voxels(fig[1,2], chunk, lowclip = :red, highclip = :red, colorrange = (1.1, 899.9), shading = NoShading)
+    foreach(a -> update_cam!(a.scene, Vec3f(0, 0, 40), Vec3f(0), Vec3f(0,1,0)), (a1, a2))
+    fig
+end
+
+@reference_test "Voxel - gap attribute" begin
+    # test direct mapping of ids to colors & upsampling of vector colormap
+    voxels(RNG.rand(3,3,3), gap = 0.3)
+end
+
 @reference_test "Plot transform overwrite" begin
     # Tests that (primitive) plots can have different transform function to their
     # parent scene (identity in this case)
@@ -548,4 +625,4 @@ end
     meshscatter!(ax, [0.1, 0.9], [0.6, 0.7], markersize = 0.05, color = :red, transformation = Transformation())
 
     fig
-end
\ No newline at end of file
+end
diff --git a/WGLMakie/assets/voxel.frag b/WGLMakie/assets/voxel.frag
new file mode 100644
index 00000000000..3b2c4fead11
--- /dev/null
+++ b/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
+flat in int plane_dim;
+flat in int plane_front;
+#endif
+
+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()
+{
+    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, uv_map, id);
+
+#ifdef DEBUG_RENDER_ORDER
+    if (plane_dim != 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 = idx.x + size.x * (idx.y + size.y * idx.z);
+        fragment_color = pack_int(object_id, lin);
+        return;
+    }
+
+    fragment_color = voxel_color;
+}
\ No newline at end of file
diff --git a/WGLMakie/assets/voxel.vert b/WGLMakie/assets/voxel.vert
new file mode 100644
index 00000000000..967e440264e
--- /dev/null
+++ b/WGLMakie/assets/voxel.vert
@@ -0,0 +1,171 @@
+// debug FLAGS
+// #define DEBUG_RENDER_ORDER
+
+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;
+
+uniform mat4 projection, view;
+
+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() {
+    get_dummy(); // otherwise this doesn't render :)
+
+    /* 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;
+    float gap = get_gap();
+    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 = get_normalmatrix() * unit_vecs[dim];
+    int dir = int(sign(dot(get_view_direction(), normal))), start;
+    if (depthsorting) {
+        // TODO: depthsorted should start far away from viewer so every plane draws
+        start = int((0.5 + 0.5 * float(dir)) * float(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 = (float((plane_idx + 1) / 2) + 0.5 * float(dir * front) * gap) * unit_vecs[dim];
+    } else {
+        // similar to above but with N+1 indices around N voxels
+        displacement = float((start + dir * id + size[dim] + 1) % (size[dim] + 1)) * unit_vecs[dim];
+    }
+
+    // place plane vertices
+    vec3 plane_vertex = vec3(size) * (orientations[dim] * get_position()) + displacement;
+    vec4 world_pos = get_model() * vec4(plane_vertex, 1.0f);
+    gl_Position = projection * view * world_pos;
+    gl_Position.z += gl_Position.w * get_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) / vec3(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);
+}
\ No newline at end of file
diff --git a/WGLMakie/src/WGLMakie.jl b/WGLMakie/src/WGLMakie.jl
index 65240bfb962..5d683b3011e 100644
--- a/WGLMakie/src/WGLMakie.jl
+++ b/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)
 
diff --git a/WGLMakie/src/serialization.jl b/WGLMakie/src/serialization.jl
index 0bd6ad2d261..3769e3b1f53 100644
--- a/WGLMakie/src/serialization.jl
+++ b/WGLMakie/src/serialization.jl
@@ -71,13 +71,20 @@ function serialize_three(p::Makie.AbstractPattern)
     return serialize_three(Makie.to_image(p))
 end
 
+three_format(::Type{<:Integer}) = "RedIntegerFormat"
 three_format(::Type{<:Real}) = "RedFormat"
 three_format(::Type{<:RGB}) = "RGBFormat"
 three_format(::Type{<:RGBA}) = "RGBAFormat"
 
+three_format(::Type{<: Makie.VecTypes{1}}) = "RedFormat"
+three_format(::Type{<: Makie.VecTypes{2}}) = "RGFormat"
+three_format(::Type{<: Makie.VecTypes{3}}) = "RGBFormat"
+three_format(::Type{<: Makie.VecTypes{4}}) = "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"
diff --git a/WGLMakie/src/voxel.jl b/WGLMakie/src/voxel.jl
new file mode 100644
index 00000000000..04290d7c48c
--- /dev/null
+++ b/WGLMakie/src/voxel.jl
@@ -0,0 +1,103 @@
+function create_shader(scene::Scene, plot::Makie.Voxels)
+    uniform_dict = Dict{Symbol, Any}()
+    uniform_dict[:voxel_id] = Sampler(plot.converted[end], minfilter = :nearest)
+    # for plane sorting
+    uniform_dict[:depthsorting] = plot.depthsorting
+    uniform_dict[:eyeposition] = Vec3f(1)
+    uniform_dict[:view_direction] = camera(scene).view_direction
+    # lighting
+    uniform_dict[:diffuse] = lift(x -> convert_attribute(x, Key{:diffuse}()), plot, plot.diffuse)
+    uniform_dict[:specular] = lift(x -> convert_attribute(x, Key{:specular}()), plot, plot.specular)
+    uniform_dict[:shininess] = lift(x -> convert_attribute(x, Key{:shininess}()), plot, plot.shininess)
+    uniform_dict[:depth_shift] = get(plot, :depth_shift, Observable(0.0f0))
+    uniform_dict[:light_direction] = Vec3f(1)
+    uniform_dict[:light_color] = Vec3f(1)
+    uniform_dict[:ambient] = Vec3f(1)
+    # picking
+    uniform_dict[:picking] = false
+    uniform_dict[:object_id] = UInt32(0)
+    # other
+    uniform_dict[: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
+    uniform_dict[:shading] = to_value(get(plot, :shading, NoShading)) != NoShading
+    uniform_dict[:gap] = lift(x -> convert_attribute(x, Key{:gap}(), Key{:voxels}()), plot.gap)
+
+    # 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
+    dummy_data = Observable(Float32[])
+    onany(plot, plot.gap, plot.converted[end]) do gap, chunk
+        N = sum(size(chunk))
+        N_instances = ifelse(gap > 0.01, 2 * N, N + 3)
+        if N_instances != length(dummy_data[]) # avoid updating unneccesarily
+            dummy_data[] = [0f0 for _ in 1:N_instances]
+        end
+        return
+    end
+    notify(plot.gap)
+
+    instance = GeometryBasics.mesh(Rect2(0f0, 0f0, 1f0, 1f0))
+
+    return InstancedProgram(WebGL(), lasset("voxel.vert"), lasset("voxel.frag"),
+                        instance, VertexArray(dummy = dummy_data), uniform_dict)
+end
diff --git a/assets/voxel_spritesheet.png b/assets/voxel_spritesheet.png
new file mode 100644
index 00000000000..b1030b995d0
Binary files /dev/null and b/assets/voxel_spritesheet.png differ
diff --git a/docs/reference/plots/voxel.md b/docs/reference/plots/voxel.md
new file mode 100644
index 00000000000..90985ff2976
--- /dev/null
+++ b/docs/reference/plots/voxel.md
@@ -0,0 +1,220 @@
+# voxels
+
+{{doc voxels}}
+
+### Examples
+
+
+
+#### Basic Example
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+# Same as volume example
+r = LinRange(-1, 1, 100)
+cube = [(x.^2 + y.^2 + z.^2) for x = r, y = r, z = r]
+cube_with_holes = cube .* (cube .> 1.4)
+
+# To match the volume example with isovalue=1.7 and isorange=0.05 we map all
+# values outside the range (1.65..1.75) to invisible air blocks with is_air
+f, a, p = voxels(-1..1, -1..1, -1..1, cube_with_holes, is_air = x -> !(1.65 <= x <= 1.75))
+```
+\end{examplefigure}
+
+
+#### Gap Attribute
+
+The `gap` attribute allows you to specify a gap size between adjacent voxels.
+It is given in units of the voxel size (at `gap = 0`) so that `gap = 0` creates no gaps and `gap = 1` reduces the voxel size to 0.
+Note that this attribute only takes effect at values `gap > 0.01`.
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+chunk = reshape(collect(1:27), 3, 3, 3)
+voxels(chunk, gap = 0.33)
+```
+\end{examplefigure}
+
+
+#### Color and the internal representation
+
+Voxels are represented as an `Array{UInt8, 3}` of voxel ids internally.
+In this representation the voxel id `0x00` is defined as an invisible air block.
+All other ids (0x01 - 0xff or 1 - 255) are visible and derive their color from the various color attributes.
+For `plot.color` specifically the voxel id acts as an index into an array of colors:
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+chunk = UInt8[
+    1 0 2; 0 0 0; 3 0 4;;;
+    0 0 0; 0 0 0; 0 0 0;;;
+    5 0 6; 0 0 0; 7 0 8;;;
+]
+f, a, p = voxels(chunk, color = [:white, :red, :green, :blue, :black, :orange, :cyan, :magenta])
+```
+\end{examplefigure}
+
+
+#### Colormaps
+
+With non `UInt8` inputs, colormap attributes (colormap, colorrange, highclip, lowclip and colorscale) work as usual, with the exception of `nan_color` which is not applicable:
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+chunk = reshape(collect(1:512), 8, 8, 8)
+
+f, a, p = voxels(chunk,
+    colorrange = (65, 448), colorscale = log10,
+    lowclip = :red, highclip = :orange,
+    colormap = [:blue, :green]
+)
+```
+\end{examplefigure}
+
+When passing voxel ids directly (i.e. an `Array{UInt8, 3}`) they are used to index a vector `[lowclip; sampled_colormap; highclip]`.
+This means id 1 maps to lowclip, 2..254 to colors of the colormap and 255 to highclip.
+`colorrange` and `colorscale` are ignored in this case.
+
+
+#### Texturemaps
+
+You can also map a texture to voxels based on their id (and optionally the direction the face is facing).
+For this `plot.color` needs to be an image (matrix of colors) and `plot.uvmap` needs to be defined.
+The `uvmap` can take two forms here.
+The first is a `Vector{Vec4f}` which maps voxel ids (starting at 1) to normalized uv coordinates, formatted left-right-bottom-top.
+
+\begin{examplefigure}{}
+```julia
+using GLMakie, FileIO
+GLMakie.activate!() # hide
+
+# load a sprite sheet with 10 x 9 textures
+texture = FileIO.load(Makie.assetpath("voxel_spritesheet.png"))
+
+# create a map idx -> LRBT coordinate of the textures, normalized to a 0..1 range
+uv_map = [
+    Vec4f(x, x+1/10, y, y+1/9)
+    for x in range(0.0, 1.0, length = 11)[1:end-1]
+    for y in range(0.0, 1.0, length = 10)[1:end-1]
+]
+
+# Define which textures/uvs apply to which voxels (0 is invisible/air)
+chunk = UInt8[
+    1 0 2; 0 0 0; 3 0 4;;;
+    0 0 0; 0 0 0; 0 0 0;;;
+    5 0 6; 0 0 0; 7 0 9;;;
+]
+
+# draw
+f, a, p = voxels(chunk, uvmap = uv_map, color = texture)
+```
+\end{examplefigure}
+
+The second format allows you define sides in the second dimension of the uvmap.
+The order of sides is: -x, -y, -z, +x, +y, +z.
+
+\begin{examplefigure}{}
+```julia
+using GLMakie, FileIO
+GLMakie.activate!() # hide
+
+texture = FileIO.load(Makie.assetpath("voxel_spritesheet.png"))
+
+# idx -> uv LRBT map for convenience. Note the change in order loop order
+uvs = [
+    Vec4f(x, x+1/10, y, y+1/9)
+    for y in range(0.0, 1.0, length = 10)[1:end-1]
+    for x in range(0.0, 1.0, length = 11)[1:end-1]
+]
+
+# Create uvmap with sides (-x -y -z x y z) in second dimension
+uv_map = Matrix{Vec4f}(undef, 4, 6)
+uv_map[1, :] = [uvs[9],  uvs[9],  uvs[8],  uvs[9],  uvs[9],  uvs[8]]  # 1 -> birch
+uv_map[2, :] = [uvs[11], uvs[11], uvs[10], uvs[11], uvs[11], uvs[10]] # 2 -> oak
+uv_map[3, :] = [uvs[2],  uvs[2],  uvs[2],  uvs[2],  uvs[2],  uvs[18]] # 3 -> crafting table
+uv_map[4, :] = [uvs[1],  uvs[1],  uvs[1],  uvs[1],  uvs[1],  uvs[1]]  # 4 -> planks
+
+chunk = UInt8[
+    1 0 1; 0 0 0; 1 0 1;;;
+    0 0 0; 0 0 0; 0 0 0;;;
+    2 0 2; 0 0 0; 3 0 4;;;
+]
+
+f, a, p = voxels(chunk, uvmap = uv_map, color = texture)
+```
+\end{examplefigure}
+
+The textures used in these examples are from [Kenney's Voxel Pack](https://www.kenney.nl/assets/voxel-pack).
+
+
+
+#### Updating Voxels
+
+The voxel plot is a bit different from other plot types which affects how you can and should update its data.
+
+First you *can* pass your data as an `Observable` and update that observable as usual:
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+chunk = Observable(ones(8,8,8))
+f, a, p = voxels(chunk, colorrange = (0, 1))
+chunk[] = rand(8,8,8)
+f
+```
+\end{examplefigure}
+
+You can also update the data contained in the plot object.
+For this you can't index into the plot though, since that will return the converted voxel id data.
+Instead you need to index into `p.args`.
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+f, a, p = voxels(ones(8,8,8), colorrange = (0, 1))
+p.args[end][] = rand(8,8,8)
+f
+```
+\end{examplefigure}
+
+Both of these solutions triggers a full replacement of the input array (i.e. `chunk`), the internal representation (`plot.converted[4]`) and the texture on gpu.
+This can be quite slow and wasteful if you only want to update a small section of a large chunk.
+In that case you should instead update your input data without triggering an update (using `obs.val`) and then call `local_update(plot, is, js, ks)` to process the update:
+
+\begin{examplefigure}{}
+```julia
+using GLMakie
+GLMakie.activate!() # hide
+
+chunk = Observable(rand(64, 64, 64))
+f, a, p = voxels(chunk, colorrange = (0, 1))
+chunk.val[30:34, :, :] .= NaN # or p.args[end].val
+Makie.local_update(p, 30:34, :, :)
+f
+```
+\end{examplefigure}
+
+
+
+#### Picking Voxels
+
+The `pick` function is able to pick individual voxels in a voxel plot.
+The returned index is a flat index into the array passed to `voxels`, i.e. `plt.args[end][][idx]` will return the relevant data.
+One important thing to note here is that the returned index is a `UInt32` internally and thus has limited range.
+Very large voxel plots (~4.3 billion voxels or 2048 x 2048 x 1024) can reach this limit and trigger an integer overflow.
diff --git a/src/Makie.jl b/src/Makie.jl
index 0859d41f837..8d41d76dcee 100644
--- a/src/Makie.jl
+++ b/src/Makie.jl
@@ -88,8 +88,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, voxels
+import MakieCore: arrows!, heatmap!, image!, lines!, linesegments!, mesh!, meshscatter!, poly!, scatter!, surface!, text!, volume!, voxels!
 import MakieCore: convert_arguments, convert_attribute, default_theme, conversion_trait
 
 export @L_str, @colorant_str
@@ -168,6 +168,7 @@ include("basic_recipes/tricontourf.jl")
 include("basic_recipes/triplot.jl")
 include("basic_recipes/volumeslices.jl")
 include("basic_recipes/voronoiplot.jl")
+include("basic_recipes/voxels.jl")
 include("basic_recipes/waterfall.jl")
 include("basic_recipes/wireframe.jl")
 include("basic_recipes/tooltip.jl")
@@ -356,9 +357,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,  Voxels
+export arrows  , heatmap  , image  , lines  , linesegments  , mesh  , meshscatter  , poly  , scatter  , surface  , text  , volume  , wireframe,  voxels
+export arrows! , heatmap! , image! , lines! , linesegments! , mesh! , meshscatter! , poly! , scatter! , surface! , text! , volume! , wireframe!, voxels!
 
 export AmbientLight, PointLight, DirectionalLight, SpotLight, EnvironmentLight, RectLight, SSAO
 
diff --git a/src/basic_recipes/voxels.jl b/src/basic_recipes/voxels.jl
new file mode 100644
index 00000000000..a015eb43715
--- /dev/null
+++ b/src/basic_recipes/voxels.jl
@@ -0,0 +1,233 @@
+function convert_arguments(T::Type{<:Voxels}, chunk::Array)
+    X, Y, Z = to_ndim(Vec3{Int}, size(chunk), 1)
+    return convert_arguments(T, -0.5X..0.5X, -0.5Y..0.5Y, -0.5Z..0.5Z, chunk)
+end
+function convert_arguments(T::Type{<:Voxels}, xs, ys, zs, chunk::Array)
+    xi = Float32(minimum(xs))..Float32(maximum(xs))
+    yi = Float32(minimum(ys))..Float32(maximum(ys))
+    zi = Float32(minimum(zs))..Float32(maximum(zs))
+    return convert_arguments(T, xi, yi, zi, chunk)
+end
+function convert_arguments(::Type{<:Voxels}, xs::ClosedInterval{Float32}, ys::ClosedInterval{Float32}, zs::ClosedInterval{Float32}, chunk::Array)
+    return (xs, ys, zs, Array{UInt8, 3}(undef, to_ndim(Vec3{Int}, size(chunk), 1)...))
+end
+function convert_arguments(::Type{<:Voxels}, xs::ClosedInterval{Float32}, ys::ClosedInterval{Float32}, zs::ClosedInterval{Float32}, chunk::Array{UInt8, 3})
+    return (xs, ys, zs, chunk)
+end
+
+function calculated_attributes!(::Type{<:Voxels}, plot)
+    if !isnothing(plot.color[])
+        cc = lift(plot, plot.color, plot.alpha) do color, a
+            if color isa AbstractVector
+                output = Vector{RGBAf}(undef, 255)
+                @inbounds for i in 1:min(255, length(color))
+                    c = to_color(color[i])
+                    output[i] = RGBAf(Colors.color(c), Colors.alpha(c) * a)
+                end
+                for i in min(255, length(color))+1 : 255
+                    output[i] = RGBAf(0,0,0,0)
+                end
+            elseif color isa AbstractArray
+                output = similar(color, RGBAf)
+                @inbounds for i in eachindex(color)
+                    c = to_color(color[i])
+                    output[i] = RGBAf(Colors.color(c), Colors.alpha(c) * a)
+                end
+            else
+                c = to_color(color)
+                output .= RGBAf(Colors.color(c), Colors.alpha(c) * a)
+            end
+            return output
+        end
+        attributes(plot.attributes)[:calculated_colors] = cc
+
+    else
+
+        # ...
+        dummy_data = Observable(UInt8[1, 255])
+
+        # Always sample N colors
+        cmap = map(plot.colormap, plot.lowclip, plot.highclip) do cmap, lowclip, highclip
+            cm = if cmap isa Vector && length(cmap) != 255
+                resample_cmap(cmap, 253)
+            else
+                categorical_colors(cmap, 253)
+            end
+            lc = lowclip === automatic ? first(cm) : to_color(lowclip)
+            hc = highclip === automatic ? last(cm) : to_color(highclip)
+            return [lc; cm; hc]
+        end
+
+        # always use 1..N
+        colorrange = Observable(Vec2f(1, 255))
+
+        # Needs to happen in voxel id generation
+        colorscale = Observable(identity)
+
+        # We always treat nan as air, invalid
+        nan_color = Observable(:transparent)
+
+        # TODO: categorical?
+        attributes(plot.attributes)[:calculated_colors] = ColorMapping(
+            dummy_data[], dummy_data, cmap, colorrange, colorscale,
+            plot.alpha, plot.lowclip, plot.highclip, nan_color
+        )
+
+    end
+
+    return nothing
+end
+
+"""
+    local_update(p::Voxels, i, j, k)
+
+Updates a section of the Voxel plot given by indices i, j, k (Integer, UnitRange
+or Colon()) according to the data present in `p.args[end]`.
+
+This is used to avoid updating the whole chunk with a pattern such as
+```
+p.args[end].val[20:30, 7:10, 8] = new_data
+local_update(plot, 20:30, 7:10, 8)
+```
+"""
+function local_update(plot::Voxels, is, js, ks)
+    to_range(N, i::Integer) = i:i
+    to_range(N, r::UnitRange) = r
+    to_range(N, ::Colon) = 1:N
+    to_range(N, x::Any) = throw(ArgumentError("Indices can't be converted to a range representation ($x)"))
+
+    _size = size(plot.converted[end].val)
+    is, js, ks = to_range.(_size, (is, js, ks))
+
+    mini, maxi = apply_scale(plot.colorscale[], plot._limits[])
+    input = plot.args[end][]
+    for k in ks, j in js, i in is
+        idx = i + _size[1] * ((j-1) + _size[2] * (k-1))
+        _update_voxel(plot.converted[end].val, input, idx, plot.is_air[], plot.colorscale[], mini, maxi)
+    end
+    plot._local_update[] = (is, js, ks)
+    return nothing
+end
+
+Base.@propagate_inbounds function _update_voxel(
+        output::Array{UInt8, 3}, input::Array, i::Integer,
+        is_air::Function, scale, mini::Real, maxi::Real
+    )
+
+    @boundscheck checkbounds(Bool, output, i) && checkbounds(Bool, input, i)
+    # Rescale data to UInt8 range for voxel ids
+    c = 252.99998
+    @inbounds begin
+        x = input[i]
+        if is_air(x)
+            output[i] = 0x00
+        else
+            lin = clamp(c * (apply_scale(scale, x) - mini) / (maxi - mini) + 2, 1, 255)
+            output[i] = trunc(UInt8, lin)
+        end
+    end
+    return nothing
+end
+
+Base.@propagate_inbounds function _update_voxel(
+        output::Array{UInt8, 3}, input::Array{UInt8, 3}, i::Integer,
+        is_air::Function, scale, mini::Real, maxi::Real
+    )
+    return nothing
+end
+
+function plot!(plot::Voxels)
+    # Disconnect automatic mapping
+    # I want to avoid recalculating limits every time the input is updated.
+    # Maybe this can be done with conversion kwargs...?
+    off(plot.args[end], plot.args[end].listeners[1][2])
+
+    # If a UInt8 Array is passed we don't do any mapping between plot.args and
+    # plot.converted. Instead we just set plot.converted = plot.args in
+    # convert_arguments
+    if eltype(plot.args[end][]) == UInt8
+        plot._limits[] = (1, 255)
+        return
+    end
+
+
+    # Use new mapping that doesn't recalculate limits
+    onany(plot, plot._limits, plot.is_air, plot.colorscale) do lims, is_air, scale
+        # _limits always triggers after plot.args[1]
+        chunk = plot.args[end][]
+        output = plot.converted[end]
+
+        # TODO: Julia doesn't allow this
+        # maybe resize
+        # if size(chunk) != size(output.val)
+        #     resize!(output.val, size(chunk))
+        # end
+
+        # update voxel ids
+        mini, maxi = apply_scale(scale, lims)
+        maxi = max(mini + 10eps(float(mini)), maxi)
+        @inbounds for i in eachindex(chunk)
+            _update_voxel(output.val, chunk, i, is_air, scale, mini, maxi)
+        end
+
+        # trigger converted
+        notify(output)
+
+        return
+    end
+
+    # Initial limits
+    map!(plot, plot._limits, plot.args[end], plot.colorrange) do data, colorrange
+        if colorrange !== automatic
+            return colorrange
+        end
+
+        mini, maxi = (Inf, -Inf)
+        for elem in data
+            plot.is_air[](elem) && continue
+            mini = min(mini, elem)
+            maxi = max(maxi, elem)
+        end
+        if !(isfinite(mini) && isfinite(maxi) && isa(mini, Real))
+            throw(ArgumentError("Voxel Chunk contains invalid data, resulting in invalid limits ($mini, $maxi)."))
+        end
+        return (mini, maxi)
+    end
+
+    return
+end
+
+function voxel_size(p::Voxels)
+    mini = minimum.(to_value.(p.converted[1:3]))
+    maxi = maximum.(to_value.(p.converted[1:3]))
+    _size = size(p.converted[4][])
+    return Vec3f((maxi .- mini) ./ _size .- convert_attribute(p.gap[], key"gap"(), key"voxels"()))
+end
+
+function voxel_positions(p::Voxels)
+    mini = minimum.(to_value.(p.converted[1:3]))
+    maxi = maximum.(to_value.(p.converted[1:3]))
+    voxel_id = p.converted[4][]
+    _size = size(voxel_id)
+    step = (maxi .- mini) ./ _size
+    return [
+        Point3f(mini .+ step .* (i-0.5, j-0.5, k-0.5))
+        for k in 1:_size[3] for j in 1:_size[2] for i in 1:_size[1]
+        if voxel_id[i, j, k] !== 0x00
+    ]
+end
+
+function voxel_colors(p::Voxels)
+    voxel_id = p.converted[4][]
+    colormapping = p.calculated_colors[]
+    uv_map = p.uvmap[]
+    if !isnothing(uv_map)
+        @warn "Voxel textures are not implemented in this backend!"
+    elseif colormapping isa ColorMapping
+        color = colormapping.colormap[]
+    else
+        color = colormapping
+    end
+
+    return [color[id] for id in voxel_id if id !== 0x00]
+end
\ No newline at end of file
diff --git a/src/camera/camera.jl b/src/camera/camera.jl
index 2d75f3dc9ca..3d70c677bd3 100644
--- a/src/camera/camera.jl
+++ b/src/camera/camera.jl
@@ -18,8 +18,8 @@ function Base.show(io::IO, camera::Camera)
     println(io, "  projection: ", camera.projection[])
     println(io, "  projectionview: ", camera.projectionview[])
     println(io, "  resolution: ", camera.resolution[])
-    println(io, "  lookat: ", camera.lookat[])
     println(io, "  eyeposition: ", camera.eyeposition[])
+    println(io, "  view direction: ", camera.view_direction[])
 end
 
 function disconnect!(c::Camera)
@@ -84,7 +84,7 @@ function Camera(viewport)
         proj,
         proj_view,
         lift(a-> Vec2f(widths(a)), viewport),
-        Observable(Vec3f(0)),
+        Observable(Vec3f(0, 0, -1)),
         Observable(Vec3f(1)),
         ObserverFunction[],
         Dict{Symbol, Observable}()
diff --git a/src/camera/camera2d.jl b/src/camera/camera2d.jl
index d9d354a180f..b73b2a000ce 100644
--- a/src/camera/camera2d.jl
+++ b/src/camera/camera2d.jl
@@ -38,6 +38,7 @@ function cam2d!(scene::SceneLike; kw_args...)
     cam = from_dict(Camera2D, cam_attributes)
     # remove previously connected camera
     disconnect!(camera(scene))
+    camera(scene).view_direction[] = Vec3f(0, 0, -1)
     add_zoom!(scene, cam)
     add_pan!(scene, cam)
     correct_ratio!(scene, cam)
@@ -343,6 +344,7 @@ controls.
 """
 function campixel!(scene::Scene; nearclip=-10_000f0, farclip=10_000f0)
     disconnect!(camera(scene))
+    camera(scene).view_direction[] = Vec3f(0, 0, -1)
     update_once = Observable(false)
     closure = UpdatePixelCam(camera(scene), nearclip, farclip)
     on(closure, camera(scene), viewport(scene))
@@ -364,6 +366,8 @@ Creates a camera for the given `scene` which maps the scene area to a 0..1 by
 0..1 range. This camera does not feature controls.
 """
 function cam_relative!(scene::Scene; nearclip=-10_000f0, farclip=10_000f0)
+    disconnect!(camera(scene))
+    camera(scene).view_direction[] = Vec3f(0, 0, -1)
     projection = orthographicprojection(0f0, 1f0, 0f0, 1f0, nearclip, farclip)
     set_proj_view!(camera(scene), projection, Mat4f(I))
     cam = RelativeCamera()
diff --git a/src/camera/camera3d.jl b/src/camera/camera3d.jl
index 6c8d72d0647..3f8ebe68e31 100644
--- a/src/camera/camera3d.jl
+++ b/src/camera/camera3d.jl
@@ -750,7 +750,7 @@ function update_cam!(scene::Scene, cam::Camera3D)
     set_proj_view!(camera(scene), proj, view)
 
     scene.camera.eyeposition[] = cam.eyeposition[]
-    scene.camera.lookat[] = cam.lookat[]
+    scene.camera.view_direction[] = normalize(cam.lookat[] - cam.eyeposition[])
 end
 
 
diff --git a/src/camera/old_camera3d.jl b/src/camera/old_camera3d.jl
index 33d381d2e45..4d1ace33bcf 100644
--- a/src/camera/old_camera3d.jl
+++ b/src/camera/old_camera3d.jl
@@ -334,6 +334,7 @@ function update_cam!(scene::Scene, cam::OldCamera3D)
     view = Makie.lookat(eyeposition, lookat, upvector)
     set_proj_view!(camera(scene), proj, view)
     scene.camera.eyeposition[] = cam.eyeposition[]
+    scene.camera.view_direction[] = normalize(cam.lookat[] - cam.eyeposition[])
 end
 
 function update_cam!(scene::Scene, camera::OldCamera3D, area3d::Rect)
diff --git a/src/conversions.jl b/src/conversions.jl
index 7b20101c507..dfe806e73bd 100644
--- a/src/conversions.jl
+++ b/src/conversions.jl
@@ -1694,6 +1694,7 @@ end
 convert_attribute(value, ::key"marker", ::key"scatter") = to_spritemarker(value)
 convert_attribute(value, ::key"isovalue", ::key"volume") = Float32(value)
 convert_attribute(value, ::key"isorange", ::key"volume") = Float32(value)
+convert_attribute(value, ::key"gap", ::key"voxels") = ifelse(value <= 0.01, 0f0, Float32(value))
 
 function convert_attribute(value::Symbol, ::key"marker", ::key"meshscatter")
     if value === :Sphere
diff --git a/src/interfaces.jl b/src/interfaces.jl
index bc72bc46ed6..13dfeeac88e 100644
--- a/src/interfaces.jl
+++ b/src/interfaces.jl
@@ -101,7 +101,7 @@ end
 
 const atomic_functions = (
     text, meshscatter, scatter, mesh, linesegments,
-    lines, surface, volume, heatmap, image
+    lines, surface, volume, heatmap, image, voxels
 )
 const Atomic{Arg} = Union{map(x-> Plot{x, Arg}, atomic_functions)...}
 
diff --git a/src/layouting/data_limits.jl b/src/layouting/data_limits.jl
index 92d4e94a38b..4980bd7caa3 100644
--- a/src/layouting/data_limits.jl
+++ b/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::Voxels)
+    xyz = to_value.(plot.converted[1:3])
+    r = Rect3f(minimum.(xyz), maximum.(xyz) .- minimum.(xyz))
+    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)
diff --git a/src/makielayout/blocks/colorbar.jl b/src/makielayout/blocks/colorbar.jl
index 9032a15a888..20887a36605 100644
--- a/src/makielayout/blocks/colorbar.jl
+++ b/src/makielayout/blocks/colorbar.jl
@@ -75,6 +75,17 @@ function extract_colormap(plot::Union{Contourf,Tricontourf})
                     elow, ehigh, plot.nan_color)
 end
 
+function extract_colormap(plot::Voxels)
+    limits = plot._limits
+    # TODO: does this need padding for lowclip and highclip?
+    discretized_values = map(lims -> range(lims[1], lims[2], length = 253), plot, limits)
+
+    return ColorMapping(
+        discretized_values[], discretized_values, plot.colormap, limits, plot.colorscale,
+        plot.alpha, plot.lowclip, plot.highclip, Observable(:transparent)
+    )
+end
+
 
 function extract_colormap_recursive(@nospecialize(plot::T)) where {T <: AbstractPlot}
     cmap = extract_colormap(plot)
diff --git a/src/types.jl b/src/types.jl
index a3144183a7b..7c8748a6044 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -241,9 +241,9 @@ struct Camera
     resolution::Observable{Vec2f}
 
     """
-    Focal point of the camera, used for e.g. camera synchronized light direction.
+    Direction in which the camera looks.
     """
-    lookat::Observable{Vec3f}
+    view_direction::Observable{Vec3f}
 
     """
     Eye position of the camera, used for e.g. ray tracing.
diff --git a/test/primitives.jl b/test/primitives.jl
index ad5a74ad1e6..2a38916cc1e 100644
--- a/test/primitives.jl
+++ b/test/primitives.jl
@@ -22,6 +22,42 @@ end
     @test cbar.limits[] ≈ Vec2f(0.5, 0.6)
 end
 
+@testset "voxels" begin
+    data = reshape(collect(range(0.3, 1.8, length=6*5*4)), 6, 5, 4)
+    f, a, p = voxels(
+        data,
+        lowclip = RGBf(1, 0, 1), highclip = RGBf(0, 1, 0),
+        colormap = [RGBf(0, 0, 0), RGBf(1, 1, 1)], gap = 0.1
+    )
+
+    # data conversion pipeline
+    @test p.args[1][] === data
+    @test p.converted[1][] ≈ -3.0..3.0
+    @test p.converted[2][] ≈ -2.5..2.5
+    @test p.converted[3][] ≈ -2.0..2.0
+
+    @test p.colorrange[] == Makie.automatic # otherwise no auto _limits
+    @test all(p._limits[] .≈ (0.3, 1.8)) # controls conversion to voxel ids
+    ids = map(data) do val
+        trunc(UInt8, clamp(2 + 253 * (val - 0.3) / (1.8 - 0.3), 2, 254))
+    end
+    @test p.converted[4][] == ids
+
+    # colormap
+    cc = p.calculated_colors[]
+    @test length(cc.colormap[]) == 255
+    @test cc.colormap[][1] == RGBAf(1,0,1,1)
+    @test cc.colormap[][2] == RGBAf(0,0,0,1)
+    @test cc.colormap[][2:end-1] == resample_cmap([RGBAf(0, 0, 0, 1), RGBAf(1, 1, 1, 1)], 253)
+    @test cc.colormap[][end-1] == RGBAf(1,1,1,1)
+    @test cc.colormap[][end] == RGBAf(0,1,0,1)
+
+    # voxels-as-meshscatter helpers
+    @test Makie.voxel_size(p) ≈ Vec3f(0.9)
+    ps = [Point3f(x - 2.5, y - 2.0, z - 1.5) for z in 0:3 for y in 0:4 for x in 0:5]
+    @test Makie.voxel_positions(p) ≈ ps
+    @test Makie.voxel_colors(p) == cc.colormap[][p.converted[end][][:]]
+end
 
 # TODO, test all primitives and argument conversions