# Makie: Isosurface colored by values from another array

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:

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:

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?

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.

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.

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).

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:

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.

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.

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.

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).

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.

`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

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`.

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

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/

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

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.

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`?

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

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)
``````

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`?

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.

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.