[Makie.jl] 3D voxelized visualization of a cubical matrix

Dear All,
Here is another chapter of the struggle to make cool figures with Makie.
My aim is to generate numerous visualizations of the 3D images that are represented by voxels. While i used the hints from the examples of the volume function, even the simplest solution for a binary (0 and 1) 3D matrix visualization still avoids me.
Here is the example of the solution i would like to achieve:


Now, here is the piece of code i could generate that i expected will produce the image above (but it does not do the trick):

using GLMakie, ColorSchemes
using GeometryBasics: Rect3f
let
    path="d:\\path\\to\\BCC_100.raw" #this is a .raw binary  example for visuazliation
    cd(dirname(path))
    cubesize=trunc(Int64,round(filesize(path)^(1/3)))
    stream = open(path, "r")
    #x, y, z = cubesize
    points3d = Array{UInt8,3}(undef, cubesize, cubesize, cubesize)
    read!(stream, points3d)
    points3d = deepcopy(points3d)
    close(stream)
    for k=1:cubesize, i=1:cubesize, j=1:cubesize
        points3d[i,j,k]=points3d[j,i,k]
    end
    x = y = z = 1:1:cubesize
    vol1 = [ix * iy * iz for ix in x, iy in y, iz in z]
    #points3d = [Point3f(ix, iy, iz) for ix in x, iy in y, iz in z]
    # scale everything to the interval 0,1 (things don't seem to work with colorrange)
    #vol2 = (vol1 .+ 1) ./ 2
    # colormap with transparency in the middle
    cmap = :Hiroshige
    colors = to_colormap(cmap, 101)
    n = length(colors)
    g(x) = x^2
    alphas = [g(x) for x in range(1, cubesize, length = n)]
    cmap_alpha = RGBAf.(colors, alphas)
    # the plot
    fig = Figure(resolution = (1000, 1000))
    #ax1 = Axis3(fig[1, 1], perspectiveness = 0.5, azimuth = 7.19,
     #   elevation = 0.57, aspect = (1, 1, 1))
    #ax2 = Axis3(fig[1, 2], perspectiveness = 0.5, azimuth = 6.62,
     #   elevation = 0.57, aspect = (1, 1, 1))
    #ax3 = Axis3(fig[2, 1], perspectiveness = 0.5, azimuth = 7.38,
     #   elevation = 0.57, aspect = (1, 1, 1))
    ax4 = Axis3(fig[1, 1], perspectiveness = 0.5, azimuth = 6.64,
        elevation = 0.57, aspect = (1, 1, 1))

    #volume!(ax1, x, y, z, vol2; colormap = cmap, transparency = true)
    #contour!(ax2, x, y, z, vol1; colormap = cmap, alpha = 0.05,
    #    levels = [collect(-1:0.01:-0.3)..., collect(0.3:0.01:1)...])
    #meshscatter!(ax3, vec(points3d); color = vec(vol1), colormap = cmap_alpha)
    meshscatter!(ax4, vec(points3d); color = vec(vol1), colormap = cmap_alpha,
        marker = Rect3f(Vec3f(-1), Vec3f(2)))
    limits!(ax4, 0, 301, 0, 301, 0, 301)
    #display(fig)
    Makie.save("test.png", fig)
end

The example raw binary (path=“d:\path\to\BCC_100.raw”) for visualization (a 100^3 matrix of 0 and 1 with the structure as on the image below) can be downloaded here:

Many thanks in advance for providing a direction of how i can make some nice visualizations (which i eventually want to extend to non binary multi-colored 3D structures).
Kirill

@KirillGerke if you are processing geospatial objects often, consider the recipes in MeshViz.jl, they really save a lot of manual adjustments if all you need is a simple command to visualize a voxel model:

Here is a MWE:

using Meshes, MeshViz

# choose a backend
import GLMakie as Mke

# domain where the data lives
grid = CartesianGrid(100, 100, 100)

# visualize a random vector of 0s and 1s
viz(grid, color = rand(0:1, 1000000))

You can pass a vector of Color objects as well instead of simple numbers. So you can use a color with full transparency for the pores and a color without transparency for the matrix.

1 Like

Thank you, Julio!
It is indeed a solution. I managed to get something like this:


Just in case somebody will try to do the same, here is my first try solution in the code:

using GLMakie, ColorSchemes
using GeometryBasics
using Meshes, MeshViz

#reading in the binary image (the .raw format as it is called)
path="you/path/binary_image.raw"
cd(dirname(path))
cubesize=trunc(Int64,round(filesize(path)^(1/3)))
stream = open(path, "r")
#x, y, z = cubesize
points3d = Array{UInt8,3}(undef, cubesize, cubesize, cubesize)
read!(stream, points3d)
close(stream)

colorvec=Vector{typeof((:red,0.5))}([])
#ok, this part is lame, but i used it to find the correct x-y-z orientation
for k=1:cubesize, i=1:cubesize, j=1:cubesize
    if points3d[i,j,k]==0
        push!(colorvec,(:black,0.01))
    else
        push!(colorvec,(:yellow,0.7))
    end
end

grid = CartesianGrid(cubesize, cubesize, cubesize)
viz(grid, color = colorvec)

Now will try to automatically save the figure with the visualization.
Again, many thanks for fast reply and solution!

You are welcome @KirillGerke! Keep in mind that you don’t need to load GeometryBasics.jl if you are relying on Meshes.jl/MeshViz.jl.

Also consider doing import GLMakie as Mke instead of loading all symbols into scope. This is to avoid conflict with GeometryBasics.jl primitives used internally in Makie.jl

1 Like

And another tip to make your code look more Julian:

colorvec = replace(points3d, 0 => (:black,0.01), 1 => (:yellow,0.7)) |> vec

You don’t need to loop manually, just use the high-level replace

1 Like

One round of cleanup:

using Meshes, MeshViz

import GLMakie as Mke

path = "/path/to/image.raw"
nvox = round(Int, filesize(path)^(1/3))
buff = zeros(UInt8, nvox, nvox, nvox)

stream = open(path, "r")
read!(stream, buff)
close(stream)

colors = replace(buff, 0 => (:black,0.01), 1 => (:yellow,0.7))

grid = CartesianGrid(nvox, nvox, nvox)
viz(grid, color = vec(colors))
1 Like

This fixed the undefined Rect3 problem i got previously (before adding GeometryBasics).
Will now update the code as you suggested.
Thank you again for all advice!

1 Like

Hello Jullio!
Thank you again a lot for help in this thread: [Makie.jl] 3D voxelized visualization of a cubical matrix - #5 by juliohm
However, just today i can’t make the same 3D visualizations - i guess due to some updates. Here are the errors i got:

ERROR: LoadError: MethodError: Cannot convert an object of type 
  Tuple{Symbol, Float64} to an object of type 
  ColorTypes.Colorant
Closest candidates are:
  convert(::Type{C}, ::C) where C<:ColorTypes.Colorant at C:\Users\M K\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:72
  convert(::Type{C}, ::ColorTypes.Colorant) where C<:ColorTypes.Colorant at C:\Users\M K\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:73
  convert(::Type{C}, ::Number) where C<:ColorTypes.Colorant at C:\Users\M K\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:74

The error is for the code in our thread:

using Meshes, MeshViz

import GLMakie as Mke

path = "/path/to/image.raw"
nvox = round(Int, filesize(path)^(1/3))
buff = zeros(UInt8, nvox, nvox, nvox)

stream = open(path, "r")
read!(stream, buff)
close(stream)

colors = replace(buff, 0 => (:black,0.01), 1 => (:yellow,0.7))

grid = CartesianGrid(nvox, nvox, nvox)
viz(grid, color = vec(colors))

Any idea how this error could be fixed?
Thank you in advance

Hi Kirill we recently updated the color handling in the recipes. Can you try passing the vector of colors with the symbols :black etc and a vector of alpha with entries in 0 and 1? For example color=[:black, :yellow] and alpha=[0.5,0.8].

I’ll check in the ColorTypes.jl specs if a color spec (:black,0.5) is officially supported and will add a conversion rule when I’m back on Monday.

that’s a makie conversion rule, which you can use with Makie.to_color(x)

Hello Simon,
Thank you for a suggestion! Well, i tried to apply it. e.g.,

fig = viz(grid, color = Mke.to_color(vec(colors))

and got something like this:

ERROR: LoadError: UndefVarError: Rect3 not defined
Stacktrace:
 [1] plot!(plot::MakieCore.Combined{MeshViz.viz, Tuple{CartesianGrid{3, Float64}}})
   @ MeshViz C:\Users\M K\.julia\packages\MeshViz\jcuPy\src\cartesiangrid.jl:38
 [2] plot!(scene::Makie.Scene, P::Type{MakieCore.Combined{MeshViz.viz, Tuple{CartesianGrid{3, Float64}}}}, attributes::MakieCore.Attributes, input::Tuple{Observables.Observable{CartesianGrid{3, Float64}}}, args::Observables.Observable{Tuple{CartesianGrid{3, Float64}}})
   @ Makie C:\Users\M K\.julia\packages\Makie\NL7Xw\src\interfaces.jl:428
 [3] plot!(scene::Makie.Scene, P::Type{MakieCore.Combined{MeshViz.viz}}, attributes::MakieCore.Attributes, args::CartesianGrid{3, Float64}; kw_attributes::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:show_axis,), Tuple{Bool}}})
   @ Makie C:\Users\M K\.julia\packages\Makie\NL7Xw\src\interfaces.jl:339
 [4] plot(P::Type{MakieCore.Combined{MeshViz.viz}}, args::CartesianGrid{3, Float64}; axis::NamedTuple{(), Tuple{}}, figure::NamedTuple{(), Tuple{}}, kw_attributes::Base.Pairs{Symbol, Vector{ColorTypes.RGBA{Float32}}, Tuple{Symbol}, NamedTuple{(:color,), Tuple{Vector{ColorTypes.RGBA{Float32}}}}})
   @ Makie C:\Users\M K\.julia\packages\Makie\NL7Xw\src\figureplotting.jl:28
 [5] #viz#1
   @ C:\Users\M K\.julia\packages\MakieCore\S8PkO\src\recipes.jl:31 [inlined]

Or did you mean something else while suggesting to use Makie.to_color?
Many thanks in advance!

MeshViz.jl has its own color conversion mechanism for any Julia object

Now that I am back to my computer, let me explain what I meant by conversion mechanism.

The following file in MeshViz.jl defines the conversion mechanism:

Users can plot any Julia object as colors by adding a new method to the ascolor trait function:

# tell how to convert object into color
ascolor(obj::MyObj, scheme) = RGB(0,0,0)

We provide implementations for vector of numbers, colors, symbols, strings, categorical values, etc. So you can simply pass

color = [:red, :green, :blue, ...]

or

color = [RGB(255,0,0), RGB(0,255,0), RGB(0,0,255), ...]

or whatever object with well-defined ascolor method. Additionally, you can pass a vector of transparencies alpha = [0.5, 0.2, 0.8, ...], which are always floats between 0 and 1 that override the transparency channel of specified colors.