Makie: Isosurface colored by values from another array

Specific Domains Visualization
,
1

I have two 3D arrays. I would like to make an isosurface from the first one and color it using the values of the second array. As far as I understand, a possible solution would be to construct a mesh for the isosurface (e.g. by marching cubes) and then, using interpolation, calculate the values of the second array at the mesh vertices.

So far I found Makie: Colouring 3D Contours by another Array discussion from where I adopted the following code:

import GLMakie as mak
import Interpolations: linear_interpolation
import MarchingCubes: MC, march
import Meshes
import MeshViz

function mesh_data(vals, isovalue)
    mc = MC(vals, Int)
    march(mc, isovalue)
    return mc
end

xmin, xmax, Nx = -10, 10, 101
ymin, ymax, Ny = -10, 10, 101
zmin, zmax, Nz = -10, 10, 101

x = range(xmin, xmax, Nx)
y = range(ymin, ymax, Ny)
z = range(zmin, zmax, Nz)

G = [zi <= 0 || sqrt(xi^2+yi^2+zi^2) <= 5 ? 1.0 : 0.0 for xi=x, yi=y, zi=z]

F = [cos(sqrt(xi^2+zi^2)) for xi=x, yi=y, zi=z]

# Extract an isosurface of G at isovalue=1
meshes = mesh_data(G, 1)
mesh_points = Meshes.Point3[pt for pt=meshes.vertices]
mesh_tris = Meshes.connect.([Tuple(tri) for tri in meshes.triangles], Meshes.Triangle)
to_SimpleMesh = Meshes.SimpleMesh(mesh_points, mesh_tris)

# Interpolate F to colour by onto the vertex locations
# Marching cubes uses cell indices starting at 0
itp = linear_interpolation((0:Nx-1, 0:Ny-1, 0:Nz-1), F)
F_at_verts = [itp(pt[1], pt[2], pt[3]) for pt=Meshes.coordinates.(to_SimpleMesh.vertices)]

fig = mak.Figure()
ax = mak.Axis3(fig[1,1]; aspect=:data)
mak.plot!(ax, to_SimpleMesh; color=F_at_verts)

mak.display(fig)

which gives the following:

out
out800×600 76.7 KB

However, this solution uses a deprecated package MeshViz and I do not understand how I can change the properties of the final plot.

Another related examples are:

  1. Visualizing a 3D volume in Makie: Lighting, smoothness
  2. ANN: Meshing 0.5.0

However, these examples do not work because Makie does not want to plot the meshes generated there.

Do you know how to realize this isosurface coloring using Makie and currently maintained packages?

2

Didn’t have the time to read the details of the code, but please take a look at our recent book, chapters 2, 10 and 11:

https://juliaearth.github.io/geospatial-data-science-with-julia

You can interpolate fields over arbitrary meshes with the Interpolate and InterpolateNeighbors transforms.

3

Thank you for the link. The book seems to be interesting. However, I can not find there any examples on how to take a 3D array and construct a mesh for its isosurface. I hope, the interpolation is the easiest step here.

4

After you construct the isosurface mesh with any method of interest, you can simply call the viz plot recipe on it. The MeshViz.jl package was converted into a package extension, which is available whenever you load Meshes.jl and Makie.jl in the same session.

The book explains how you can visualize different fields over the mesh (geometry column).

5

Oh, thank you. MeshViz.jl held the older version of Meshes.jl. When I removed it, the viz function became available.
The updated version of the code:

import GLMakie
import Interpolations
import MarchingCubes
import Meshes

function mesh_data(x, y, z, vals; isovalue=1)
    mc = MarchingCubes.MC(vals, Int; x, y, z)
    MarchingCubes.march(mc, isovalue)
    return mc
end

# Prepare data:
x = collect(range(-10, 10, 101))
y = collect(range(-10, 10, 101))
z = collect(range(-10, 10, 101))
G = [zi <= 0 || sqrt(xi^2+yi^2+zi^2) <= 5 ? 1.0 : 0.0 for xi=x, yi=y, zi=z]
F = [cos(sqrt(xi^2+zi^2)) for xi=x, yi=y, zi=z]

# Extract the isosurface of G at isovalue=1:
mc = mesh_data(x, y, z, G; isovalue=1)
mesh_points = [Meshes.Point3(vt) for vt=mc.vertices]
mesh_triangles = Meshes.connect.([Tuple(tri) for tri=mc.triangles], Meshes.Triangle)
simple_mesh = Meshes.SimpleMesh(mesh_points, mesh_triangles)

# Interpolate F at vertex locations:
itp = Interpolations.linear_interpolation((x, y, z), F)
Fcolors = [itp(pt[1], pt[2], pt[3]) for pt=Meshes.coordinates.(simple_mesh.vertices)]

# Visualize:
fig = GLMakie.Figure()
ax = GLMakie.Axis3(fig[1,1]; aspect=:data)
Meshes.viz!(ax, simple_mesh; color=Fcolors)
GLMakie.display(fig)

and the output:

out
out800×600 78.1 KB

However, I would prefer the original Makie visualization functions. Do you know what kind of mesh object I have to provide to Makie.mesh function in order to make it work? Do you know where I can find the sources of viz function? In Meshes.jl I found only its template.

And slightly unrelated questions. Do you know any good resource where I can find how to control the light sources in Makie? There are not so much information in the documentation.

6

The viz recipe is defined in Meshes/ext for different mesh types using dispatch. If you want to keep things at a high level I would stick to it and follow the recommendations in the book. Feel free to go low level and use Makie.jl directly with points and connec lists of triangles.

7

A volume plot like volume(x, y, z, G, isovalue = 1.0) should give you the same shape but it doesn’t allow you to paint F on the surface. (I think that’s what’s happening here?) This also includes lighting. I’d guess your examples don’t because they don’t generate normals or because the normals face the other direction.

With this specific example you could also reorganize your data to have z(x, y) and F(x, y) and plot it as a surface(x, y, z_xy, color = F_xy). If that’s not an option you’ll have to generate a mesh like in your examples.

8

Yes, the main issue is to paint the surface using an external array.

I managed to find out how to control the lights by set_theme! function:

GLMakie.set_theme!(
    lightposition = GLMakie.Vec3f(1,1,1),
    ambient = GLMakie.RGBf(1,1,1),
    # backlight = 1,
    # shininess = 128f0,
    # specular = GLMakie.RGBf(1,1,0),
    # diffuse = GLMakie.RGBf(1,0,1),
)

(though it seems that the commented parameters do nothing).

out
out800×600 105 KB

9

If ambient is 1 then there’s maximum light everywhere already. If backlight is on, the other parameters should become visible no matter the normal orientation.

10

backlight might get ignored by the backend if it’s not a Float32. It wasn’t set up with a conversion originally and I’m not sure if that ever got fixed

11

Yes, I do not see any the effects of other parameters even if I set ambient = GLMakie.RGBf(0.5,0.5,0.5) and backlight = true.

12

@jules, @ffreyer maybe you know how to prepare a proper mesh for Makie.mesh function. The use of Meshes.viz is very unconvinient. For example, I can not change the colorrange attribute no matter what I do. Or I can not change the colormap through set_theme!; it seems Meshes.viz redefines it.

13

Does this help?
https://docs.makie.org/stable/reference/plots/mesh/#using_geometrybasicsmesh_and_buffersampler_type
Btw, we found out that the lighting isn’t all that correct at the moment and @ffreyer has been fixing it here

If you want to give that a try, I think it also contains some better docs:
https://docs.makie.org/previews/PR3246/reference/scene/lighting/

14

@fedoroff take a look at the docstring of the viz function. It explains the options and rationale. If some keyword argument is missing that should be forwarded to the Makie.jl internal call, feel free to suggest the addition.

15

Not really. Unfortunately, now I do not have enough knowledge to understand what is going on there. I will try to read more on how GeometryBasics define meshes. Thank you.

The examples there look really amazing. I hope we will see it soon in the production.

16

The docstring of viz is very bref. For example, I want to change the colorrange attribute to enhance the colors. The docstring does not mention this argument. If I write something like
Meshes.viz!(ax, simple_mesh; color=Fcolors, colorscheme, colorrange=(0,0.1)), nothing happens. Does it mean that the Meshes owners have to manually add every keyword argument which is available in the original Makie function?

And, sorry, but I still can not find the exact place in Meshes sources where viz is defined. I searched through the sources using “viz” keyword and found nothing. Can you, please, give me a direct link where I can see what is the internal machinery of viz?

17

The color range is set automatically based on the extrema(colors). This is what most users would expect in practice, specially when they are using the viewer function that adjusts the color bar accordingly when a different variable is selected from the menu. Can you please clarify what you mean by enhacing the colors?

The viz recipe is a high-level utility function that simply forwards all defined options to Makie.jl internal functions. The current list of options was defined based on real needs from users, and we can certainly brainstorm more options if a new use case shows up. It is not clear to me yet what do you mean by “enhance colors”.

The ext folder contains all the viz methods as I explained above. The specific method you are looking for SimpleMesh is written in simplemesh.jl

18

When the data has a very large dynamic range (a lot of details at small values) or it is corrupted by, for example, some “broken pixels”, then, like in the example below, by changing the colorrange we can visualize some specific details which are not visible with the default color scale. In other words, we can oversaturate the colors to force some tiny details to be visible.

import GLMakie as mak

x = range(-10, 10, 101)
y = range(-10, 10, 101)
F = [0.1*exp(-(xi^2+yi^2)/4^2) + 1e-4*sin((xi^2+yi^2)/2^2)  for xi=x, yi=y]
F[25,25] = 1   # broken pixel

fig = mak.Figure(resolution=(1500,600))

ax1 = mak.Axis(fig[1,1]; title="Default colorrange")
hm1 = mak.heatmap!(ax1, x, y, F)
mak.Colorbar(fig[2,1], hm1; vertical=false, flipaxis=false)

ax2 = mak.Axis(fig[1,2]; title="colorrange=(0,0.1)")
hm2 = mak.heatmap!(ax2, x, y, F; colorrange=(0,0.1))
mak.Colorbar(fig[2,2], hm2; vertical=false, flipaxis=false)

ax3 = mak.Axis(fig[1,3]; title="colorrange=(0,0.001)")
hm3 = mak.heatmap!(ax3, x, y, F; colorrange=(0,0.001))
mak.Colorbar(fig[2,3], hm3; vertical=false, flipaxis=false)

mak.display(fig)

out
out1500×600 59.2 KB

In the file simplemesh.jl I see only the following functions:
Makie.plottype(::SimpleMesh) = Viz{<:Tuple{SimpleMesh}}
function Makie.plot!(plot::Viz{<:Tuple{SimpleMesh}})
function vizmesh1D!(plot)
function vizmesh2D!(plot)
function vizmesh3D!(plot)
function segmentsof(topo, vert)
function segmentsof(topo::GridTopology, vert)

I would expect that when I call Meshes.viz function, there is somewhere an entry point defined as function viz(some_arguments). And I can not find it. What does exactly happens when I call Meshes.viz?

19

Perhaps that could be achieved with a different scale other than the linear scale? Sometimes log scale can address some of these visualization issues. We could also say that the issue is in the data, and that you should preprocess it first before sending to viz. Check the DropExtrema, Scale and similar transforms in chapter 5:

https://juliaearth.github.io/geospatial-data-science-with-julia/05-transforms.html

This is the Makie.jl recipe system. The keyword arguments are defined in the plot object. If you read the code you will find lines accessing plot[:colorscheme] for example.

20

I believe that the data is primary and thus I would prefer not to preprocess it (besides, in some cases it can take sufficient amount of time), but instead use already existing tools for its visualization. That is why I would like to use Makie directly and not rely on third-part wrappers that cut the functionality of the original functions. I assume that viz covers all your needs. However, I need a little bit more.

Unfortunately, I am not familiar with the recipe system. It looks for me like some black magic. There is a function which you can call, but its body is defined nowhere. I have to read more about it.