Makie: Colouring 3D Contours by another Array

Specific Domains Visualization
1

I would like to build a 3D contour plot, with the contours constructed based on a given array and the colours on said contours defined by a second array. An example relevant to me is the Inviscid Taylor-Green Vortex:

using GLMakie

function InitConds(x, y, z)
    γ = 1.4; R = 287.15
    vx = sin(x) * cos(y) * cos(z)
    vy = -cos(x) * sin(y) * cos(z)
    vz = 0
    p = 1e5 + ((cos(2 * z) + 2) * (cos(2 * x) + cos(2 * y)) - 2) / 16
    ρ = 1
    T = p / (R * ρ)
    return ρ * sqrt(vx^2 + vy^2 + vz^2), ρ * (log(T^(γ / (γ - 1))) / log(p))
end

domain_size = 2π
npts = 64; dx = domain_size / npts
points = LinRange(dx / 2, domain_size - dx / 2, npts)

data = [InitConds(x, y, z) for x = points, y = points, z = points]
KinEnergy = getindex.(data, 1)
Entropy = getindex.(data, 2)

fig = Figure()
ax = Axis3(fig[1, 1], title = "Inviscid Taylor Green Vortex")
contour!(ax, points, points, points, KinEnergy, levels = 5)
display(fig)

At the moment, this plots contours of KinEnergy, with the contours coloured by the value of KinEnergy. I would like the contours to be coloured by Entropy- is this possible? I haven’t been able to find an example in the Makie docs or in BeautifulMakie.

As an aside, is there a more elegant way to pull out values from broadcasted functions that return multiple values than what I did with KinEnergy and Entropy?

2

There’s no builtin way to specify colors for the levels, as it’s really a volume plot and they don’t work with “discrete” colors. There might be a way to do it by precomputing a colormap where each level in the plot falls on a color that would be produced if the level’s value had been taken not from KinEnergy but from Entropy. Like contour!(..., colormap = [some_color for level in levels])

3

Ah that’s a shame. That method would prescribe a single colour for each level right? I think it would be possible by algorithmically locating the coordinates of the desired contours and plotting them as surfaces then colour those surfaces, with some time and effort.

4

You can get/color your contours according to the values of the Entropy, using MarchingCubes.jl, respectively PlotlyJS.jl.
PlotlyJS can plot volumes and contours, but not according to the values of an additional array.
You can use instead MarchingCubes to carve each contour into the volume, KinEnergy, and PlotlyJS to display it as a 3d Mesh, colored according to the values of Entropy at its points.
Below is defined an arbitrary function (sin(z)) to color each contour, just to see how it works.
Hence you have to get the values of entropy at the points of coordinates
(xv[i], yv[i], zv[i]), i=1:N, through interpolation.

using PlotlyJS, MarchingCubes
function mesh_data(fvals::Array{T,3};  v₀=0.0) where T <: AbstractFloat
    #instantiate the structure MC (Marching Cube)
    mc = MC(fvals, Int)
    m = march(mc, v₀) #v₀ is the isovalue
    #extract data for the 3d mesh to be plotly plotted
    xv, yv, zv = eachrow(reduce(hcat, mc.vertices))    
    I, J, K = eachrow(reduce(hcat, mc.triangles) .- 1) #0-based  indices for triangles
    return xv, yv, zv, I, J, K
end
begin
    domain_size = 2π
    npts = 64; dx = domain_size / npts
    points = LinRange(dx / 2, domain_size - dx / 2, npts)

    data = [InitConds(x, y, z) for x = points, y = points, z = points]
    KinEnergy = getindex.(data, 1)
    Entropy = getindex.(data, 2)
    pl_data= AbstractTrace[]
    for v0 in [0.1, 0.3, 0.5, 0.7, 0.9]
        xv, yv, zv, I, J, K = mesh_data(KinEnergy; v₀= v0)
        push!(pl_data, mesh3d(x=xv, y=yv, z=zv,
                     i=I, j=J, k=K,
                     intensity=sin.(zv),#here would go the Entropy evaluated, through interpolation, at the contour points
                     colorscale=colors.plasma,   
                     showscale=false,
                     opacity=0.6,
                     lighting=attr(ambient=0.65,
                                   diffuse=0.5,
                                   fresnel=0.25,        
                                   specular=0.25,
                                   roughness=0.23),
                     lightposition=attr(x=1000, y=1000, z=100)
             ))
    end            
    layout = Layout(
                 width=600,
                 height=600,
                 scene=attr( 
                    aspectmode="data",
                    camera_eye=attr(x=1.6, y=1.15, z=0.8),
                    ))
         
    pl = Plot(pl_data, layout) 
end

contours

1 Like
5

Two more lines must be added:

emin, emax= extrema(Entropy) #before `for v0 in...`
cmin=emin, cmax=emax, # after opacity=0.6

to ensure that points of the same entropy are mapped to the same color in the colorscale, no matter the contour they belong.

6

Ah yeah if you do marching cubes, you should be able to use Makie’s mesh as well :slight_smile:

7

That method using Plotly worked well. I’m trying to make it work in Makie, having issues making the mesh formats match what Makie expects. I’ll post a solution here if I can make it work.

8

Worked out how to do it in Makie, using the MarchingCubes, Meshes and MeshViz packages. Meshes to convert the output of MarchingCubes to a Julia Mesh, and MeshViz to extend the plot() signature to accept the Mesh.

using GLMakie, Interpolations, Meshes, MeshViz, MarchingCubes

function InitConds(x, y, z)
    γ = 1.4; R = 287.15
    vx = sin(x) * cos(y) * cos(z)
    vy = -cos(x) * sin(y) * cos(z)
    vz = 0
    p = 1e5 + ((cos(2 * z) + 2) * (cos(2 * x) + cos(2 * y)) - 2) / 16
    ρ = 1
    T = p / (R * ρ)
    return ρ * sqrt(vx^2 + vy^2 + vz^2), ρ * (log(T^(γ / (γ - 1))) / log(p))
end

# Don't unpack the locations and connections
function mesh_data(fvals::Array{T, 3}; v0 = 0.0) where T <: AbstractFloat
    # Generate the marching cube
    mc = MC(fvals, Int)
    m = march(mc, v0)
    return mc
end

begin
    # Configure the domain
    domain_size = 2π
    npts = 64; dx = domain_size / npts
    points = LinRange(dx / 2, domain_size - dx / 2, npts)

    # Generate the data
    data = [InitConds(x, y, z) for x = points, y = points, z = points]
    KinEnergy = getindex.(data, 1)
    Entropy = getindex.(data, 2)

    # Extract an isosurface of the KinEnergy at 3/4 the maximum value
    meshes = mesh_data(KinEnergy, v0 = 0.75 * maximum(KinEnergy))
    mesh_points = Meshes.Point3[Tuple(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 the field to colour by onto the vertex locations
    # Marching cubes uses cell indices starting at 0
    itp = interpolate((0:npts-1, 0:npts-1, 0:npts-1), Entropy, Gridded(Linear()))
    Entropy_at_verts = [itp(pt[1], pt[2], pt[3]) for pt = coordinates.(to_SimpleMesh.points)]

    fig = Figure()
    ax = Axis3(fig[1, 1], title = "Inviscid Taylor Green Vortex")
    plot!(ax, to_SimpleMesh, color = Entropy_at_verts)
    display(fig)
end

Results in this:

ContourColouredDemo
ContourColouredDemo800×600 195 KB

2 Likes
9

You can probably also do the mesh part with Makie’s dependency GeometryBasics, although I’m not a 3d guy so I don’t remember how exactly :slight_smile: But you probably don’t need Meshes and Meshviz just for that unless it runs some complex algorithm that GeometryBasics doesn’t have

10

I suspect it is, but the documentation of Makie’s mesh() function is very brief and I couldn’t work it how it fit in with this MarchingCubes mesh. The MarchingCubes package falls under the JuliaGeometry umbrella so they should work together in some form, but I am far from an expert on mesh generation for visualisation. The MWE from the MarchingCubes docs actually use Meshes and MeshViz.