Cleanup volume (#4726)

* add indexedabsorption and additive volume algorithms to WGLMakie

* fix mip for negative values

* fix RGB(A) volume plots no accepting 3(4) component data

* use absorption in absorptionrgba, indexedabsorption

* document what algorithm do

* update clip planes to test volume with more fitting data

* avoid newline in docstirng list

* fix automatic x/y/z values for Color data

* fix and test absorption

* update attribute docs

* add rgba, index examples

* update changelog

* don't show algorithm docs twice

* fix typo

* skip volume test in CairoMakie

* fix merge error

Author: ffreyer

Diff for: CHANGELOG.md

@@ -16,7 +16,8 @@
 ## [0.22.1] - 2025-01-17
 
 - Allow volume textures for mesh color, to e.g. implement a performant volume slice display [#2274](https://github.com/MakieOrg/Makie.jl/pull/2274).
-- Fix `alpha` use in legends and some CairoMakie cases [#4721](https://github.com/MakieOrg/Makie.jl/pull/4721).
+- Fixed `alpha` use in legends and some CairoMakie cases [#4721](https://github.com/MakieOrg/Makie.jl/pull/4721).
+- Cleaned up `volume` plots: Added `:indexedabsorption` and `:additive` to WGLMakie, generalized `:mip` to include negative values, fixed missing conversions for rgba algorithms (`:additive`, `:absorptionrgba`), fixed missing conversion for `absorption` attribute & extended it to `:indexedabsorption` and `absorptionrgba`, added tests and improved docs. [#4726](https://github.com/MakieOrg/Makie.jl/pull/4726)
 
 ## [0.22.0] - 2024-12-12
 

Diff for: CairoMakie/test/runtests.jl

@@ -194,6 +194,7 @@ excludes = Set([
     "picking", # Not implemented
     "MetaMesh (Sponza)", # makes little sense without per pixel depth order
     "Mesh with 3d volume texture", # Not implemented yet
+    "Volume absorption",
 ])
 
 functions = [:volume, :volume!, :uv_mesh]

Diff for: GLMakie/assets/shader/volume.frag

@@ -35,45 +35,30 @@ const float max_distance = 1.3;
 const int num_samples = 200;
 const float step_size = max_distance / float(num_samples);
 
-float _normalize(float val, float from, float to)
-{
-    return (val-from) / (to - from);
-}
+float _normalize(float val, float from, float to) { return (val-from) / (to - from);}
 
-vec4 color_lookup(float intensity, Nothing color_map, Nothing norm, vec4 color)
-{
+vec4 color_lookup(float intensity, Nothing color_map, Nothing norm, vec4 color) {
     return color;
 }
-
-vec4 color_lookup(float intensity, samplerBuffer color_ramp, vec2 norm, Nothing color)
-{
+vec4 color_lookup(float intensity, samplerBuffer color_ramp, vec2 norm, Nothing color) {
     return texelFetch(color_ramp, int(_normalize(intensity, norm.x, norm.y)*textureSize(color_ramp)));
 }
-
-vec4 color_lookup(float intensity, samplerBuffer color_ramp, Nothing norm, Nothing color)
-{
+vec4 color_lookup(float intensity, samplerBuffer color_ramp, Nothing norm, Nothing color) {
     return vec4(0);  // stub method
 }
-
-vec4 color_lookup(float intensity, sampler1D color_ramp, vec2 norm, Nothing color)
-{
+vec4 color_lookup(float intensity, sampler1D color_ramp, vec2 norm, Nothing color) {
     return texture(color_ramp, _normalize(intensity, norm.x, norm.y));
 }
-
-vec4 color_lookup(samplerBuffer colormap, int index)
-{
-    return texelFetch(colormap, index);
+vec4 color_lookup(vec4 data_color, Nothing color_ramp, Nothing norm, Nothing color) {
+    return data_color;  // stub method
 }
-
-vec4 color_lookup(sampler1D colormap, int index)
-{
-    return texelFetch(colormap, index, 0);
+vec4 color_lookup(float intensity, Nothing color_ramp, Nothing norm, Nothing color) {
+    return vec4(0);  // stub method
 }
 
-vec4 color_lookup(Nothing colormap, int index)
-{
-    return vec4(0);
-}
+vec4 color_lookup(samplerBuffer colormap, int index) { return texelFetch(colormap, index); }
+vec4 color_lookup(sampler1D colormap, int index) { return texelFetch(colormap, index, 0); }
+vec4 color_lookup(Nothing colormap, int index) { return vec4(0); }
 
 vec3 gennormal(vec3 uvw, float d, vec3 o)
 {
@@ -184,7 +169,7 @@ vec4 absorptionrgba(vec3 front, vec3 dir)
     int i = 0;
     for (i; i < num_samples ; ++i) {
         vec4 density = texture(volumedata, pos);
-        float opacity = step_size * density.a;
+        float opacity = step_size * density.a * absorption;
         T *= 1.0-opacity;
         if (T <= 0.01)
             break;
@@ -204,7 +189,7 @@ vec4 volumeindexedrgba(vec3 front, vec3 dir)
     for (i; i < num_samples; ++i) {
         int index = int(texture(volumedata, pos).x) - 1;
         vec4 density = color_lookup(color_map, index);
-        float opacity = step_size*density.a;
+        float opacity = step_size*density.a * absorption;
         Lo += (T*opacity)*density.rgb;
         T *= 1.0 - opacity;
         if (T <= 0.01)
@@ -284,9 +269,9 @@ vec4 isosurface(vec3 front, vec3 dir)
 
 vec4 mip(vec3 front, vec3 dir)
 {
-    vec3 pos = front;
-    int i = 0;
-    float maximum = 0.0;
+    vec3 pos = front + dir;
+    int i = 1;
+    float maximum = texture(volumedata, front).x;
     for (i; i < num_samples; ++i, pos += dir){
         float density = texture(volumedata, pos).x;
         if(maximum < density)
@@ -315,7 +300,7 @@ bool process_clip_planes(inout vec3 p1, inout vec3 p2)
             p2 = p1;
             return true;
         }
-        
+
         // one outside - shorten segment
         else if (d1 < 0.0)
             // solve 0 = m * t + b = (d2 - d1) * t + d1 with t in (0, 1)

Diff for: GLMakie/src/drawing_primitives.jl

@@ -959,6 +959,9 @@ function draw_atomic(screen::Screen, scene::Scene, plot::Volume)
         if haskey(gl_attributes, :intensity)
             intensity = pop!(gl_attributes, :intensity)
             return draw_volume(screen, Tex(intensity), gl_attributes)
+        elseif haskey(gl_attributes, :color)
+            color = pop!(gl_attributes, :color)
+            return draw_volume(screen, Tex(color), gl_attributes)
         else
             return draw_volume(screen, Tex(plot[4]), gl_attributes)
         end

Diff for: GLMakie/src/glshaders/image_like.jl

@@ -62,9 +62,9 @@ function draw_volume(screen, main::VolumeTypes, data::Dict)
         volumedata = main => Texture
         model = Mat4f(I)
         modelinv = const_lift(inv, model)
-        color_map = default(Vector{RGBA}, s) => Texture
-        color_norm = color_map === nothing ? nothing : const_lift(extrema2f0, main)
-        color = color_map === nothing ? default(RGBA, s) : nothing
+        color_map = nothing => Texture
+        color_norm = nothing
+        color = nothing => Texture
 
         algorithm = MaximumIntensityProjection
         absorption = 1f0

Diff for: MakieCore/src/basic_plots.jl

@@ -275,32 +275,38 @@ end
     volume(volume_data)
     volume(x, y, z, volume_data)
 
-Plots a volume, with optional physical dimensions `x, y, z`.
-Available algorithms are:
-* `:iso` => IsoValue
-* `:absorption` => Absorption
-* `:mip` => MaximumIntensityProjection
-* `:absorptionrgba` => AbsorptionRGBA
-* `:additive` => AdditiveRGBA
-* `:indexedabsorption` => IndexedAbsorptionRGBA
+Plots a volume with optional physical dimensions `x, y, z`.
+
+All volume plots are derived from casting rays for each drawn pixel. These rays
+intersect with the volume data to derive some color, usually based on the given
+colormap. How exactly the color is derived depends on the algorithm used.
 """
 @recipe Volume (
         x::EndPoints,
         y::EndPoints,
         z::EndPoints,
-        volume::AbstractArray{Float32,3}
+        # TODO: consider using RGB{N0f8}, RGBA{N0f8} instead of Vec/RGB(A){Float32}
+        volume::AbstractArray{<: Union{Float32, Vec3f, RGB{Float32}, Vec4f, RGBA{Float32}}, 3}
     ) begin
-    "Sets the volume algorithm that is used."
+    """
+    Sets the volume algorithm that is used. Available algorithms are:
+    * `:iso`: Shows an isovalue surface within the given float data. For this only samples within `isovalue - isorange .. isovalue + isorange` are included in the final color of a pixel.
+    * `:absorption`: Accumulates color based on the float values sampled from volume data. At each ray step (starting from the front) a value is sampled from the volume data and then used to sample the colormap. The resulting color is weighted by the ray step size and blended the previously accumulated color. The weight of each step can be adjusted with the multiplicative `absorption` attribute.
+    * `:mip`: Shows the maximum intensity projection of the given float data. This derives the color of a pixel from the largest value sampled from the respective ray.
+    * `:absorptionrgba`: This algorithm matches :absorption, but samples colors directly from RGBA volume data. For each ray step a color is sampled from the data, weighted by the ray step size and blended with the previously accumulated color. Also considers `absorption`.
+    * `:additive`: Accumulates colors using `accumulated_color = 1 - (1 - accumulated_color) * (1 - sampled_color)` where `sampled_color` is a sample of volume data at the current ray step.
+    * `:indexedabsorption`: This algorithm acts the same as :absorption, but interprets the volume data as indices. They are used as direct indices to the colormap. Also considers `absorption`.
+    """
     algorithm = :mip
-    "Sets the target value for the IsoValue algorithm."
+    "Sets the target value for the :iso algorithm. `accepted = isovalue - isorange < value < isovalue + isorange`"
     isovalue = 0.5
-    "Sets the range of values picked up by the IsoValue algorithm."
+    "Sets the maximum accepted distance from the isovalue for the :iso algorithm. `accepted = isovalue - isorange < value < isovalue + isorange`"
     isorange = 0.05
     "Sets whether the volume data should be sampled with interpolation."
     interpolate = true
-    "Enables depth write for Volume, so that volume correctly occludes other objects."
+    "Enables depth write for :iso so that volume correctly occludes other objects."
     enable_depth = true
-    "Absorption multiplier for algorithm=:absorption. This changes how much light each voxel absorbs."
+    "Absorption multiplier for algorithm = :absorption, :absorptionrgba and :indexedabsorption. This changes how much light each voxel absorbs."
     absorption = 1f0
     mixin_generic_plot_attributes()...
     mixin_shading_attributes()...

Diff for: ReferenceTests/src/tests/examples3d.jl

@@ -77,6 +77,17 @@ end
     f
 end
 
+@reference_test "Volume absorption" begin
+    f = Figure(size = (600, 300))
+    r = range(-5, 5, length=31)
+    data = [cos(x*x + y*y + z*z)^2 for x in r, y in r, z in r]
+    absorption = 5.0
+    volume(f[1, 1], data, algorithm = :absorption; absorption)
+    volume(f[1, 2], 128 .+ 120 .* data, algorithm = :indexedabsorption; absorption)
+    volume(f[1, 3], HSV.(180 .* data, 0.8, 0.9), algorithm = :absorptionrgba; absorption)
+    f
+end
+
 @reference_test "Textured meshscatter" begin
     catmesh = loadasset("cat.obj")
     img = loadasset("diffusemap.png")
@@ -661,11 +672,15 @@ end
 end
 
 @reference_test "Clip planes - volume" begin
-    f = Figure(size = (600, 400))
+    f = Figure(size = (600, 400), backgroundcolor = :black)
     r = -10:10
     data = [1 - (1 + cos(x^2) + cos(y^2) + cos(z^2)) for x in r, y in r, z in r]
-    clip_planes = [Plane3f(Vec3f(-1), 0.0)]
+    index_data = round.(Int, 10 .* abs.(data))
+    N = maximum(index_data)
+    density_data = 0.005 .* abs.(data)
+    rgba_data = [RGBAf(cos(x^2)^2, cos(y^2)^2, cos(z^2)^2, 0.5 + 0.5 * sin(x^2 + y^2 + z^2)) for x in r, y in r, z in r]
 
+    clip_planes = [Plane3f(Vec3f(-1), 0.0)]
     attr = (clip_planes = clip_planes, axis = (show_axis = false,))
 
     volume(f[1, 1], -10..10, -10..10, -10..10, data; attr...,
@@ -675,13 +690,14 @@ end
 
     volume(f[1, 2], -10..10, -10..10, -10..10, data; attr...,
         algorithm = :mip)
-    volume(f[2, 2], -10..10, -10..10, -10..10, data; attr...,
+    volume(f[2, 2], -10..10, -10..10, -10..10, rgba_data; attr...,
         algorithm = :absorptionrgba)
 
-    volume(f[1, 3], -10..10, -10..10, -10..10, data; attr...,
-        algorithm = :additive)
-    volume(f[2, 3], -10..10, -10..10, -10..10, data; attr...,
-        algorithm = :indexedabsorption)
+    # TODO: doesn't work as intended anymore?
+    volume(f[1, 3], -10..10, -10..10, -10..10, rgba_data; attr...,
+        algorithm = :additive, alpha = 0.01)
+    volume(f[2, 3], -10..10, -10..10, -10..10, index_data; attr...,
+        algorithm = :indexedabsorption, colormap = Makie.resample(to_colormap(:viridis), N))
 
     f
 end

Diff for: WGLMakie/assets/volume.frag

@@ -13,15 +13,27 @@ const float step_size = max_distance / float(num_samples);
 
 uniform vec4 clip_planes[8];
 
-float _normalize(float val, float from, float to)
-{
-    return (val-from) / (to - from);
-}
+float _normalize(float val, float from, float to) { return (val-from) / (to - from); }
 
-vec4 color_lookup(float intensity, sampler2D color_ramp, vec2 norm)
-{
+vec4 color_lookup(float intensity, sampler2D color_ramp, vec2 norm) {
     return texture(color_ramp, vec2(_normalize(intensity, norm.x, norm.y), 0.0));
 }
+vec4 color_lookup(vec4 color, bool color_ramp, bool norm) {
+    return color; // stub method
+}
+vec4 color_lookup(float intensity, bool color_ramp, bool norm) {
+    return vec4(0); // stub method
+}
+
+vec4 color_lookup(sampler2D colormap, int index) {
+    return texelFetch(colormap, ivec2(index, 0), 0);
+}
+vec4 color_lookup(bool colormap, vec4 color) {
+    return color; // stub method
+}
+vec4 color_lookup(bool colormap, int index) {
+    return vec4(0); // stub method
+}
 
 vec3 gennormal(vec3 uvw, float d)
 {
@@ -110,15 +122,15 @@ vec4 volume(vec3 front, vec3 dir)
 }
 
 
-vec4 volumergba(vec3 front, vec3 dir)
+vec4 absorptionrgba(vec3 front, vec3 dir)
 {
     vec3  pos = front;
     float T = 1.0;
     vec3 Lo = vec3(0.0);
     int i = 0;
     for (i; i < num_samples ; ++i) {
         vec4 density = texture(volumedata, pos);
-        float opacity = step_size * density.a;
+        float opacity = step_size * density.a * absorption;
         T *= 1.0 - opacity;
         if (T <= 0.01)
             break;
@@ -179,9 +191,9 @@ vec4 isosurface(vec3 front, vec3 dir)
 
 vec4 mip(vec3 front, vec3 dir)
 {
-    vec3 pos = front;
-    int i = 0;
-    float maximum = 0.0;
+    vec3 pos = front + dir;
+    int i = 1;
+    float maximum = texture(volumedata, front).x;
     for (i; i < num_samples; ++i, pos += dir){
         float density = texture(volumedata, pos).x;
         if(maximum < density)
@@ -190,6 +202,38 @@ vec4 mip(vec3 front, vec3 dir)
     return color_lookup(maximum, colormap, colorrange);
 }
 
+vec4 additivergba(vec3 front, vec3 dir)
+{
+    vec3 pos = front;
+    vec4 integrated_color = vec4(0., 0., 0., 0.);
+    int i = 0;
+    for (i; i < num_samples ; ++i) {
+        vec4 density = texture(volumedata, pos);
+        integrated_color = 1.0 - (1.0 - integrated_color) * (1.0 - density);
+        pos += dir;
+    }
+    return integrated_color;
+}
+
+vec4 volumeindexedrgba(vec3 front, vec3 dir)
+{
+    vec3 pos = front;
+    float T = 1.0;
+    vec3 Lo = vec3(0.0);
+    int i = 0;
+    for (i; i < num_samples; ++i) {
+        int index = int(texture(volumedata, pos).x) - 1;
+        vec4 density = color_lookup(colormap, index);
+        float opacity = step_size * density.a * absorption;
+        Lo += (T*opacity)*density.rgb;
+        T *= 1.0 - opacity;
+        if (T <= 0.01)
+            break;
+        pos += dir;
+    }
+    return vec4(Lo, 1.0 - T);
+}
+
 uniform uint objectid;
 
 const float typemax = 100000000000000000000000000000000000000.0;
@@ -290,10 +334,11 @@ void main()
     else if(algorithm == uint(2))
         color = mip(start, step_in_dir);
     else if(algorithm == uint(3))
-        color = volumergba(start, step_in_dir);
+        color = absorptionrgba(start, step_in_dir);
     else if(algorithm == uint(4))
-        color = vec4(0.0);
-        // color = volumeindexedrgba(start, step_in_dir);
+        color = additivergba(start, step_in_dir);
+    else if(algorithm == uint(5))
+        color = volumeindexedrgba(start, step_in_dir);
     else
         color = contours(start, step_in_dir);
 

Diff for: WGLMakie/src/meshes.jl

@@ -44,7 +44,7 @@ function handle_color!(plot, uniforms, buffers, uniform_color_name = :uniform_co
             ) do uvt, pv, res, model, pattern
             return Makie.pattern_uv_transform(uvt, pv * model, res, pattern, true)
         end
-    elseif color[] isa AbstractMatrix
+    elseif color[] isa Union{AbstractMatrix, AbstractArray{<: Any, 3}}
         uniforms[uniform_color_name] = Sampler(convert_texture(color); minfilter=minfilter)
     elseif color[] isa Makie.ColorMapping
         if color[].color_scaled[] isa AbstractVector