[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

1 Like

@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:

https://github.com/JuliaGeometry/MeshViz.jl

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.

2 Likes

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))
2 Likes

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.

Julio,

Finally i set down to rework that. Before i just used an older version of the environment and the old code worked…

Now the working code looks like this:

using Meshes, MeshViz

import GLMakie as Mke

path = "/path/to/file.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, 1 => :orange)
alpha = replace(buff, 0 => 0.05, 1 => 0.5)

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

Btw, a question about crediting yours and Simon’s work in creating the packages. For Simon i could google this one:
Danisch, S., & Krumbiegel, J. (2021). Makie. jl: Flexible high-performance data visualization for Julia. Journal of Open Source Software, 6(65), 3349.
But to credit your input, i did not find anything better than to refer to this one:
Hoffimann, J. (2018). GeoStats. jl–High-performance geostatistics in Julia. Journal of Open Source Software, 3(24), 692.

Again, thank you for all the answers!

3 Likes

Hi Kirill, thanks for considering giving credit to us, that is very kind of you. On my side there is active work in Meshes.jl and MeshViz.jl. The latter is just a set of visualization recipes optimized for the geometries and meshes of the former. These recipes are built with the Makie project that you already have a reference for.

I should probably add a DOI for Meshes.jl at some point. The reference you found on GeoStats.jl is for a more ambitious project on top of Meshes.jl.

1 Like