Update 3D contours

CHANGELOG.md

@@ -5,6 +5,7 @@
 - Widened types for axis keys [#5243](https://github.com/MakieOrg/Makie.jl/pull/5243)
 - Fixed `getlimits(::Axis3)` error related to unchecked access of `:visible` attribute.
 - Add simple compression for arrays containing only the same value in WGLMakie [#5252](https://github.com/MakieOrg/Makie.jl/pull/5252).
+- Fixed 3D `contour` plots not rendering the correct isosurfaces when `colorrange` is given. Also fixed `isorange` not working, tweaked default `isorange`, colormap resolution, and changed colormap extractor for `Colorbar` to ignore alpha. [#5213](https://github.com/MakieOrg/Makie.jl/pull/5213)
 
 ## [0.24.5] - 2025-08-06
 

Makie/src/basic_recipes/contours.jl

@@ -45,6 +45,10 @@ If only `z::Matrix` is supplied, the indices of the elements in `z` will be used
     labelformatter = contour_label_formatter
     "Font size of the contour labels"
     labelsize = 10 # arbitrary
+    """
+    Sets the tolerance for sampling of a `level` in 3D contour plots.
+    """
+    isorange = automatic
     mixin_colormap_attributes()...
     mixin_generic_plot_attributes()...
 end
@@ -111,40 +115,82 @@ conversion_trait(::Type{<:Contour}) = VertexGrid()
 conversion_trait(::Type{<:Contour}, x, y, z, ::Union{Function, AbstractArray{<:Number, 3}}) = VolumeLike()
 conversion_trait(::Type{<:Contour}, ::AbstractArray{<:Number, 3}) = VolumeLike()
 
+# 3D Contour
+
 function plot!(plot::Contour{<:Tuple{X, Y, Z, Vol}}) where {X, Y, Z, Vol}
-    x, y, z, volume = plot[1:4]
-    @extract plot (colormap, levels, linewidth, alpha)
-    valuerange = lift(nan_extrema, plot, volume)
-    cliprange = map((v, default) -> ifelse(v === automatic, default, v), plot, plot.colorrange, valuerange)
-    cmap = lift(plot, colormap, levels, alpha, cliprange, valuerange) do _cmap, l, alpha, cliprange, vrange
-        levels = to_levels(l, vrange)
-        nlevels = length(levels)
-        N = 50 * nlevels
-
-        iso_eps = if haskey(plot, :isorange)
-            plot.isorange[]
+    map!(nan_extrema, plot, :converted_4, :value_range)
+    map!(default_automatic, plot, [:colorrange, :value_range], :tight_colorrange)
+
+    map!(to_levels, plot, [:levels, :value_range], :value_levels)
+
+    # the default isorange should be smaller than the gap between levels, but not
+    # so small that surfaces disappear/get skipped
+    map!(plot, [:isorange, :value_levels, :value_range], :computed_isorange) do isorange, value_levels, (min, max)
+        if isorange === automatic
+            if length(value_levels) > 1
+                minstep = minimum(value_levels[2:end] .- value_levels[1:(end - 1)])
+                return 0.1 * minstep
+            else
+                return 0.1 * (max - min)
+            end
         else
-            nlevels * ((vrange[2] - vrange[1]) / N) # TODO calculate this
+            return isorange
+        end
+    end
+
+    # The colorrange and colormap needs to be padded with RGBAf(..., 0) so that
+    # samples outside the colorrange are not drawn
+    map!(plot, [:tight_colorrange, :computed_isorange], :padded_colorrange) do (min, max), isorange
+        return (min - 2isorange, max + 2isorange)
+    end
+
+    map!(plot, [:value_levels, :tight_colorrange], :clamped_levels) do levels, (min, max)
+        return filter(lvl -> min <= lvl <= max, levels)
+    end
+
+    map!(to_colormap, plot, :colormap, :input_colormap)
+
+    map!(
+        plot,
+        [:clamped_levels, :tight_colorrange, :padded_colorrange, :computed_isorange, :alpha, :input_colormap],
+        :computed_colormap
+    ) do levels, tight_colorrange, (min, max), isorange, alpha, cmap
+        # We need colormap values for the full color range (with padding)
+        # We also need enough color values to have samples in
+        # `level - isorange .. level + isorange`, otherwise we might skip over
+        # isosurfaces
+        # GLMakie texture size is typically limited 8192+
+        # WGLMakie texture size may be limited to 4096+
+        N = ceil(Int, 2.5 * (max - min) / isorange)
+        if N > 4096
+            min_isorange = (max - min) / 4096
+            @warn "Isorange maybe too small to resolve iso surfaces. Try `isorange > $min_isorange`"
         end
-        cmap = to_colormap(_cmap)
-        v_interval = cliprange[1] .. cliprange[2]
-        # resample colormap and make the empty area between iso surfaces transparent
-        map(1:N) do i
-            i01 = (i - 1) / (N - 1)
-            c = Makie.interpolated_getindex(cmap, i01)
-            isoval = vrange[1] + (i01 * (vrange[2] - vrange[1]))
-            line = reduce(levels, init = false) do v0, level
-                isoval in v_interval || return false
-                v0 || abs(level - isoval) <= iso_eps
+        N = clamp(N, 100, 4096)
+
+        clip_range = tight_colorrange[1] - isorange .. tight_colorrange[2] + isorange
+        return map(1:N) do i
+            isoval = min + (i - 1) / (N - 1) * (max - min)
+            c = Colors.color(interpolated_getindex(cmap, isoval, tight_colorrange))
+            if isoval in clip_range && any(lvl -> lvl - isorange < isoval < lvl + isorange, levels)
+                return RGBAf(c, alpha)
+            else
+                return RGBAf(c, 0.0)
             end
-            RGBAf(Colors.color(c), line ? alpha : 0.0)
         end
     end
 
-    return volume!(
-        plot, Attributes(plot), x, y, z, volume, alpha = 1.0, # don't apply alpha 2 times
-        algorithm = 7, colorrange = cliprange, colormap = cmap
+    volume!(
+        plot, Attributes(plot),
+        plot.converted_1, plot.converted_2, plot.converted_3, plot.converted_4,
+        alpha = 1.0, # don't apply alpha 2 times
+        algorithm = 7, # contour algorithm
+        colorrange = plot.padded_colorrange,
+        colormap = plot.computed_colormap,
+        isorange = 0.0 # unused, but needs to be a float
     )
+
+    return plot
 end
 
 color_per_level(color, args...) = color_per_level(to_color(color), args...)

Makie/src/makielayout/blocks/colorbar.jl

@@ -72,6 +72,21 @@ function extract_colormap(plot::Union{Contourf, Tricontourf})
     )
 end
 
+function extract_colormap(plot::Contour{<:Tuple{X, Y, Z, Vol}}) where {X, Y, Z, Vol}
+    levels = ComputePipeline.get_observable!(plot.value_levels)
+    # Users may use transparency to make layered isosurfaces visible. Because
+    # 3D contours often accumulate the color of an isosurface over multiple
+    # samples one typically needs very low alpha values for this, which would
+    # make the colors in the colormap very faint. To keep the Colorbar useful,
+    # we remove user alpha here. (The recipe also uses `alpha = 0` to remove
+    # samples outside of isosurfaces. This is preserved here)
+    colormap = map(cm -> RGBAf.(Colors.color.(cm), Colors.alpha.(cm) .> 0.0f0), plot.computed_colormap)
+    return ColorMapping(
+        levels[], levels, colormap, plot.padded_colorrange, plot.colorscale,
+        Observable(1.0), Observable(automatic), Observable(automatic), plot.nan_color
+    )
+end
+
 function extract_colormap(plot::Voxels)
     limits = plot.value_limits
     # TODO: does this need padding for lowclip and highclip?

ReferenceTests/src/tests/examples3d.jl

@@ -136,7 +136,7 @@ end
     limits!(ax, Rect3f(-1, -1, -1, 2, 2, 2))
 
     ax = Axis3(f[2, 1])
-    contour!(ax, -5 .. 5, -5 .. 5, -5 .. 5, r, levels = [0.5, 0.9, 1.9])
+    contour!(ax, -5 .. 5, -5 .. 5, -5 .. 5, r, levels = [0.5, 0.9, 1.8])
     limits!(ax, Rect3f(-1, -1, -1, 2, 2, 2))
 
     ax = Axis3(f[3, 1])
@@ -239,14 +239,22 @@ end
     end
     x = range(-2pi, stop = 2pi, length = 100)
     # c[4] == fourth argument of the above plotting command
-    fig, ax, c = contour(x, x, x, test, levels = 6, alpha = 0.3, transparency = true)
+    fig = Figure(size = (400, 700))
+    ax, c = contour(fig[1, 1], x, x, x, test, levels = 6, alpha = 0.03, colormap = [:white, :black], transparency = true)
 
     xm, ym, zm = minimum(data_limits(c))
     contour!(ax, x, x, map(v -> v[1, :, :], c[4]), transformation = (:xy, zm), linewidth = 2)
     heatmap!(ax, x, x, map(v -> v[:, 1, :], c[4]), transformation = (:xz, ym))
     contourf!(ax, x, x, map(v -> v[:, :, 1], c[4]), transformation = (:yz, xm))
     # reorder plots for transparency
     ax.scene.plots[:] = ax.scene.plots[[1, 3, 4, 5, 2]]
+
+    contour(
+        fig[2, 1], x, x, x, (x, y, z) -> sqrt(x * x + y * y) / (10 + z * z),
+        levels = [0.01, 0.1, 0.2, 0.5, 1.0],
+        colorrange = (0.1, 0.5), # this should clip 0.01 and 1.0
+    )
+
     fig
 end
 

WGLMakie/test/runtests.jl

@@ -27,7 +27,6 @@ excludes = Set(
         "Image on Surface Sphere", # TODO: texture rotated 180°
         "Array of Images Scatter", # scatter does not support texture images
         "Order Independent Transparency",
-        "3D Contour with 2D contour slices", # looks like a z-fighting issue
         "Mesh with 3d volume texture", # Not implemented yet
         "matcap", # not yet implemented
     ]

docs/src/reference/plots/contour.md

@@ -76,6 +76,44 @@ ctr = contour!(ax, xs, ys, zs; color = :orange, levels = levels, labels = true,
 fig
 ```
 
+### 3D contours
+
+3D contour plots exist in two variants.
+`contour` implements a variant showing multiple isosurfaces, i.e. surfaces that sample the same value from a 3D array.
+[contour3d](@ref) computes the same isolines as a 2D `contour` plot but renders them in 3D at z values equal to their level.
+
+```@figure backend=GLMakie
+r = range(-pi, pi, length = 21)
+data2d = [cos(x) + cos(y) for x in r, y in r]
+data3d = [cos(x) + cos(y) + cos(z) for x in r, y in r, z in r]
+
+f = Figure(size = (700, 400))
+a1 = Axis3(f[1, 1], title = "3D contour()")
+contour!(a1, -pi .. pi, -pi .. pi, -pi .. pi, data3d)
+
+a2 = Axis3(f[1, 2], title = "contour3d()")
+contour3d!(a2, r, r, data2d, linewidth = 3, levels = 10)
+f
+```
+
+```@figure backend=GLMakie
+r = range(-pi, pi, length = 21)
+data3d = [cos(x) + cos(y) + cos(z) for x in r, y in r, z in r]
+
+f = Figure(size = (700, 300))
+
+# isorange controls the thickness of isosurfaces
+# Note that artifacts may appear if isorange becomes too small (< 0.03 here)
+a1 = Axis3(f[1, 1])
+contour!(a1, -pi .. pi, -pi .. pi, -pi .. pi, data3d, isorange = 0.04)
+
+# small alpha can be used to see into the contour plot
+a2 = Axis3(f[1, 2])
+contour!(a2, -pi .. pi, -pi .. pi, -pi .. pi, data3d, alpha = 0.05)
+f
+```
+
+
 ## Attributes
 
 ```@attrdocs

docs/src/reference/plots/contour3d.md

@@ -7,6 +7,24 @@ contour3d
 
 ## Examples
 
+3D contour plots exist in two variants.
+[contour](@ref) implements a variant showing multiple isosurfaces, i.e. surfaces that sample the same value from a 3D array.
+`contour3d` computes the same isolines as a 2D `contour` plot but renders them in 3D at z values equal to their level.
+
+```@figure backend=GLMakie
+r = range(-pi, pi, length = 21)
+data2d = [cos(x) + cos(y) for x in r, y in r]
+data3d = [cos(x) + cos(y) + cos(z) for x in r, y in r, z in r]
+
+f = Figure(size = (700, 400))
+a1 = Axis3(f[1, 1], title = "3D contour()")
+contour!(a1, -pi .. pi, -pi .. pi, -pi .. pi, data3d)
+
+a2 = Axis3(f[1, 2], title = "contour3d()")
+contour3d!(a2, r, r, data2d, linewidth = 3, levels = 10)
+f
+```
+
 ```@figure backend=GLMakie
 f = Figure()
 Axis3(f[1, 1], aspect=(0.5,0.5,1), perspectiveness=0.75)