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?
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
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
)
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
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.
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))
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?
Try transparency = true or overdraw = true for the purple line?
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
Thanks, looks much crisper now:
