Makie: Isosurface colored by values from another array

You can read the code and extract the underlying Makie command for additional customization. The viz function is part of a framework, i.e., a set of guidelines for working with geospatial data. If you don’t adhere to the framework, then it is certainly a bit painful to use the function in isolation.

Good luck with the additional customization.

Finally, I managed to plot the mesh using the native Makie functions. The key point was to convert Meshes.SimpleMesh to the coordinates and faces recognizable by Makie.mesh function. Though this conversion is present in the source code of Meshes, it is buried inside a sophisticated (for me) code. I found it in its original form here.
The final solution is the following:

import GLMakie as mak
import Interpolations
import MarchingCubes
import Meshes

# 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 for a given isovalue:
isovalue = 1
mc = MarchingCubes.MC(G, Int; x, y, z)
MarchingCubes.march(mc, isovalue)

# Prepare the mesh using Meshes:
mesh_points = @. Meshes.Point3(mc.vertices)
mesh_triangles = @. Meshes.connect(Tuple(mc.triangles), Meshes.Triangle)
mesh = Meshes.SimpleMesh(mesh_points, mesh_triangles)

# Extract coordinates and faces recognizable by GLMakie.mesh:
mesh_vertices = Meshes.vertices(mesh)
xyz = @. collect(Meshes.coordinates(mesh_vertices))
xyz = transpose(reduce(hcat, xyz))
mesh_elements = Meshes.elements(Meshes.topology(mesh))
faces = @. collect(Meshes.indices(mesh_elements))
faces = transpose(reduce(hcat, faces))

# Convert F to colors by interpolating it onto the coordinates of the mesh vertices:
itp = Interpolations.linear_interpolation((x, y, z), F)
Fcolors = [itp(xyz[i,1], xyz[i,2], xyz[i,3]) for i=1:size(xyz,1)]

# Visualize:
mak.set_theme!(lightposition=mak.Vec3f(1,1,1), ambient=mak.RGBf(1,1,1))
fig = mak.Figure()
ax = mak.Axis3(fig[1,1]; aspect=:data)
mak.mesh!(ax, xyz, faces; color=Fcolors)
mak.display(fig)
# mak.save("out.png", fig)


Thank you everyone for your help!

You shouldn’t be using Meshes.jl if all you need is a list of vertices and a list of connectivity tuples as your low-level representation. The main benefit of the package is in the abstraction layer that allows you to perform additional operations on the surface, representation of n-gons, etc.

1 Like

Sorry, I am not a specialist in meshes.
Can you show how I can extract the coordinates and faces for the Makie.mesh function directly from the MarchingCubes data?

You need to simply check the output of MarchingCubes.jl, it should be a list of vertices and connecitivities already. They show examples in the README and you can always use the Julia REPL to inspect the object. Take advantage of the dynamic nature of the language to explore the fields of the structs.

1 Like

Thank you for your very valuable comment.

An updated solution without Meshes:

import GeometryBasics
import GLMakie as mak
import Interpolations
import MarchingCubes

# 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 for a given isovalue:
isovalue = 1
mc = MarchingCubes.MC(G, Int; x, y, z)
MarchingCubes.march(mc, isovalue)
mesh = MarchingCubes.makemesh(GeometryBasics, mc)

# Convert F to colors by interpolating it onto the coordinates of the mesh vertices:
itp = Interpolations.linear_interpolation((x, y, z), F)
Fcolors = [itp(pt[1], pt[2], pt[3]) for pt=mc.vertices]

# Visualize:
mak.set_theme!(lightposition=mak.Vec3f(1,1,1), ambient=mak.RGBf(1,1,1))
fig = mak.Figure()
ax = mak.Axis3(fig[1,1]; aspect=:data)
mak.mesh!(ax, mesh; color=Fcolors)
mak.display(fig)
# mak.save("out.png", fig)
1 Like