Plotting an ellipsoid as 2D surface in 3D, given the principal axes

Hi there,

Let’s say I have some lengths [5.0, 2.0, 1.0] and I want to plot an ellipsoid whose axis have these lengths as a surface in 3D. Extra bonus if I can also provide two angles of rotation for azimuth / polar angles for the orientation of the ellipsoid.

How do I do this in Makie.jl?

1 Like

Isn’t there maybe a sphere primitive in geometrybasics that you could transform to a mesh, which you then transform to an ellipse with a matrix multiplication?

Would really appreciate some help, as GeometryBasics is a bit hard. I have:

using GLMakie, AbstractPlotting, StaticArrays
sphere = Sphere(Point3f0(0), 1)
spheremesh = AbstractPlotting.GeometryBasics.mesh(sphere)

fig = Figure(resolution = (800, 600))
ax3D = fig[1,1] = Axis3(fig)
mesh!(ax3D, spheremesh)

this works, but is pure black

I’m trying to color the surface but

c = [RGBAf0(rand(3)...) for i in 1:length(spheremesh)]
mesh!(ax3D, spheremesh; color = c)

errors claiming that lengths do not match, which I’m not sure it makes sense.

Nevertheless, I really don’t know how to transform this spheremesh and then plot it either, because:

mat = rand(3,3)
newmesh = [mat*point for point in spheremesh]

errors

I think meshes are iterated as triangles, so it’s slightly more complicated than just multiplying points. And can’t you set color = :red or something? The length if you pass an array should probably be equal to the number of vertices, not faces.

A few days back I was discussing with Andrew Bylard about the same thing… but our problem was to properly rotate the ellipsoid…I (we) didn’t succeed, nevertheless here is the code.

using GLMakie, Rotations
fig = Figure(resolution = (600,600))
box_center = Vec(0.,0.,0.)
box_widths = Vec(10.,10.,10.)
box = Rect(box_center - box_widths/2, box_widths)
r = [2., 3., 5.]
p_center = zeros(3)
ax = Axis3(fig, aspect = :data)
mesh!(ax, box, color = (:dodgerblue, 0.5), transparency = true)
sphere = Sphere(Point(p_center...), 1)

meshplot = mesh!(ax, sphere, color = :red)

scale!(meshplot, r[1], r[2], r[3])
obj =  ax.scene[end]

rotation = RotY(π/2)
q = UnitQuaternion(rotation)
qrot = Quaternionf0(q.w, q.x, q.y, q.z) # Quaternionf0(0.5, 0.5, 0.6, 0.5)
GLMakie.rotate!(obj, qrot) # 

fig[1,1] = ax
fig

1 Like

A version for “uncertainty balls”

using Distributions, StaticArrays, LinearAlgebra
using GLMakie

"""
    quaternion(m::SMatrix{3,3})

Compute the (rotation-) quarternion of a 3x3 rotation matrix. Useful to create
isodensity ellipses from spheres in GL visualizations.
"""
function quaternion(m)
    qw2 = 1 + m[1,1] + m[2,2] + m[3,3]
    if qw2 < 0
        qw2 = -qw2
    end

    qw = sqrt(qw2)/2
    qx = (m[3,2] - m[2,3])/(4qw)
    qy = (m[1,3] - m[3,1])/(4qw)
    qz = (m[2,1] - m[1,2])/(4qw)
    SVector{4}(qx,qy,qz,qw)
end

"""
    svd3(mcstates) -> mean, q, sv

Decompose variance into singular values `sv` and 
3d rotation quaternion `q`. 

Returns vectors of appropriate `GeometryTypes`.
"""
function svd3(vs)
    scal = Point3f0[]
    rot = Vec4f0[]
    for v in vs
        s = svd(v)
        push!(scal, sqrt.(s.S))
        push!(rot, quaternion(s.U))
    end
    rot, scal
end

qu = 0.95

xs = cos.(1:0.5:20)
ys = sin.(1:0.5:20)
zs = LinRange(0, 3, length(xs))
outer(x) = x*x'
rot, scal = svd3([0.005outer(randn(3,3) + 0.2I) for i in 1:length(xs)])

c = sqrt(quantile(Chisq(3), qu))

meshscatter(xs, ys, zs, color = zs,
    markersize = c*(scal),      
    rotations = rot,        
    markerspace=SceneSpace 
)

5 Likes

Thank you both, I’ll use your answers. I am not sure how I can take advantage of the answer of @mschauer , because I only need to plot 1 ellipsoid whose color changes along the ellipsoid. To be precise, I want to color it so that it is red at the top and green at the bottom, etc. (I have a map from coordinates to colors already)

I tried out something as well, although I always struggle with the GeometryBasics interface. Never know how I’m supposed to do things like transformations, accessing points, constructing meshes, etc.

Anyway:

import AbstractPlotting.GeometryBasics as GB

f = Figure()

sphere = Sphere(Point3f0(0), 1)
spheremesh = GB.mesh(sphere)

function transform_mesh(msh, mat4x4)
    pos_trans = Point3f0.(Ref(mat4x4) .* to_ndim.(Point4f0, spheremesh.position, 0))
    GB.Mesh(pos_trans, GB.faces(msh))
end

m = AbstractPlotting.rotationmatrix_y(pi/4) *
        AbstractPlotting.scalematrix(Vec3f0(1, 1, 2))

colors = [randn() for p in spheremesh.position]

ax = Axis3(f[1, 1], aspect = :data)
mesh!(transform_mesh(spheremesh, m), color = colors)

f

1 Like

something like this: (this code you can find it the gallery beautiful makie,

using AbstractPlotting: get_dim, surface_normals
using GLMakie, Rotations, GeometryBasics

function getMesh(x,y,z)
    positions = vec(map(CartesianIndices(z)) do i
    GeometryBasics.Point{3, Float32}(
        get_dim(x, i, 1, size(z)),
        get_dim(y, i, 2, size(z)),
        z[i])
    end)
    faces = decompose(GLTriangleFace, Rect2D(0f0, 0f0, 1f0, 1f0), size(z))
    normals = surface_normals(x, y, z)
    vertices = GeometryBasics.meta(positions; normals=normals) 
    meshObj = GeometryBasics.Mesh(vertices, faces)
    meshObj 
end


Θ = LinRange(0, 2π, 100) # 50
Φ = LinRange(0, π, 100)
r = 0.5
x = [r * cos(θ) * sin(ϕ) for θ in Θ, ϕ in Φ]
y = [r * sin(θ) * sin(ϕ) for θ in Θ, ϕ in Φ]
z = [r * cos(ϕ) for θ in Θ, ϕ in Φ]


meshSphere = getMesh(x,y, z)

fig = Figure(resolution = (600,600))

r = [0.5, 0.6, 0.95]

ax = Axis3(fig, aspect = :data)

meshplot = mesh!(ax, meshSphere, color= [v[3] for v in coordinates(meshSphere)], # color = z, # v[3]
    colormap = (:green,:red), shading = true)
scale!(meshplot, r[1], r[2], r[3])

fig[1,1] = ax
fig

I suppose the right way to do the transformations, rotations should be done as @mschauer just mentioned above.

1 Like

this solution is nice (@jules):
for your case @Datseris, just change

colors = [p[3] for p in spheremesh.position]
and add
colormap = (:green, :red) in mesh!(…). Although the surface looks too opaque. Probably something to do with normals.

Using Makie and Rotations.jl:

using Makie, Rotations

M(u,v) = [a/2*cos(u)*sin(v), b/2*sin(u)*sin(v), c/2*cos(v)]    # ellipsoid

RM(u,v) = RotXYZ(α,β,γ) * M(u,v) .+ C     # rotated ellipsoid

C = [1,1,2]                 # center of ellipsoid
α, β, γ = π/4, 0, π/3       # Euler angles
a, b, c = 5.0, 2.0, 1.0     # ellipsoid principal axes
u, v = range(0, 2π, length=72), range(0, π, length=72)
xs, ys, zs = [[p[i] for p in RM.(u, v')] for i in 1:3]
Makie.surface(xs,ys,zs, lightposition = Vec3f0(0,-10,0))

4 Likes

Thank you very much @jules , in the end your code was the one that was most useful to me and suitable to my application. With it was able to make the following educative illustration:

(only thing left for me to do is find out how the purple line in left plot can be plotted always on top of the grey lines)

This discussion made me realize that this is not something a user could find on their own. Should we have some kind of documentation page about this online? If yes, where? Do we have some kind of “official examples library” in the documentation, like PyPlot does?

3 Likes

Try transparency = true or overdraw = true for the purple line?

1 Like

also using tightlimits!(ax, Bottom()) with the barplots would probably look nicer. And you can try transparency = true on the lorenz attractor lines as well, they lose their white jagged edges that way

1 Like

Thanks, looks much crisper now:

4 Likes