Drawing an ellipsoid to visualize a tensor

[KirillGerke]

Dear All,

I am struggling with some problem in writing my code in Julia and, thus, ask for your assistance…
What i want to do: visualize a tensor (3x3 matrix) as ellipsoid
Problem: 1) not sure about which graphical framework to use, e.g., PyPlot or Makie (the latter looks very cool, but so far i did not get it)
2) how to draw an ellipsoid for my tensor (see the code below)

Here is my code so far:

function draw_ellipsoid(K::Array{Float64,2},discr=50)
    F = svd(K)
    U = F.U
    V = F.V
    a = F.S[1]
    b = F.S[2]
    c = F.S[3]
    C = [0 0 0] #center of the ellipsoid

    #Here is the line in my matlab code i can't rewrite into Julia in any meaningful way... This ellipsoid function is the one i used to utilize in Matlab (https://www.mathworks.com/help/matlab/ref/ellipsoid.html?searchHighlight=ellipsoid&s_tid=doc_srchtitle)...
    [X,Y,Z] = ellipsoid(C[2],C[1],C[3],a,b,c,discr);

    XX = zeros(discr+1,discr+1);
    YY = zeros(discr+1,discr+1);
    ZZ = zeros(discr+1,discr+1);

    for k = 1:length(X), j = 1:length(X)
        point = [X(k,j) Y(k,j) Z(k,j)]'
        P = V * point;
        XX(k,j) = P[1];
        YY(k,j) = P[2];
        ZZ(k,j) = P[3];
    end
    #Here is another question - should i do it with PyPlot or may be Makie (Makie i do not understand at all, bit saw nice figs)
    surf(XX,YY,ZZ)
end

K =[  0.070087    -0.0386396   0.00209088;
 -0.0441344    0.105735    0.00012994;
 -0.00553273  -0.00292807  0.0896236]

The results should look something like this:

I would greatly appreciate any help!

[ElOceanografo]

This seems to work:

using Plots
pyplot()
using LinearAlgebra

function ellipsoid(center, coefs, ngrid=25)
    # Radii corresponding to the coefficients:
    rx, ry, rz = 1 ./ sqrt.(coefs)

    # Set of all spherical angles:
    u = range(0, 2pi, length=ngrid)
    v = range(0, pi, length=ngrid)

    # Cartesian coordinates that correspond to the spherical angles:
    # (this is the equation of an ellipsoid):
    x = [rx * x * y for (x, y) in  Iterators.product(cos.(u), sin.(v))]
    y = [ry * x * y for (x, y) in Iterators.product(sin.(u), sin.(v))]
    z = [rz * x * y for (x, y) in Iterators.product(ones(length(u)), cos.(v))]
    return x .+ center[1], y .+ center[2], z .+ center[3]
end
# Plot:

K =[  0.070087    -0.0386396   0.00209088;
 -0.0441344    0.105735    0.00012994;
 -0.00553273  -0.00292807  0.0896236]

F = svd(K)
U = F.U
V = F.V
C = [0 0 0]

x, y, z = ellipsoid(C, F.S)
surface(x, y, z, aspect_ratio=:equal)

The code above is a more-or-less direct translation of the Python in this StackOverflow answer, so there’s likely a more Julian way to do it if you want to get creative…

[KirillGerke]

Thank you, Sam!

While your solution provides a possibility to plot an ellipsoid, where is something wrong in the ellipsoid visualization generated by the code. Your ellipsoid is like a spheroid, which would mean that the tensor K contains only diagonal values. But in the K example in the code we have quite strong off-diagonal terms, thus, the ellispoid should be oriented somewhat angled to the orthogonal direction. Also, while K has the max magnitude of about 0.1, the ellipsoid produced goes up to 4…
At first glance i could not spot the problem, will try to look deeper…

[c42f]

Here’s a script which might help if you’d like to try Makie. Makie takes a while to load, but by gosh it is slickly interactive!

Before going on, I’d like to note that visualizing this tensor K of yours may be tricky because it is not symmetric. You tackle that in your example by taking svd(K) and considering only the singular values, but it’s worth thinking about whether this is what you want for your application. (Are the matrices U and V also important for you? They will define two different rotated frames)

With all those questions in mind, in the following I’ve cooked up a visualization which measures the length of K*v for each vector v on the unit sphere, and also color coded the output vector K*v in the color. Maybe this is not at all what you want though! (If you could describe your application briefly perhaps we could come up with a more appropriate visualization.)

using Makie
using LinearAlgebra

# Utility to split an array-of-length-3-arrays into x,y,z components
split_coords(ps) = getindex.(ps,1), getindex.(ps,2), getindex.(ps,3)

# Utility to convert vectors to colors
to_rgb(v) = RGBf0(v[1], v[2], v[3])

# Convert spherical coordinates to Cartesian coordinates.
# Note that θ,ϕ are named according to the "physicist convention"
# Note that using Vec3f0 here is efficient, but you could just use a
# normal array and everything would still work.
spherical_to_cartesian(r,θ,ϕ) = Vec3f0(r*sin(θ)*cos(ϕ), r*sin(θ)*sin(ϕ), r*cos(θ))

θ = range(0,pi,length=20)
ϕ = permutedims(range(0,2pi,length=30))

K =[ 0.070087    -0.0386396   0.00209088;
    -0.0441344    0.105735    0.00012994;
    -0.00553273  -0.00292807  0.0896236]

# A parametric unit sphere
v = spherical_to_cartesian.(1.0, θ, ϕ)

# One way of "seeing how K acts in each direction" is by taking the norm of
# K * v and using this as the radius in terms of the original angular
# coordinates.
r = norm.(Ref(K) .* v)
ps = spherical_to_cartesian.(r, θ, ϕ)

# Another way is by seeing how it distorts the unit sphere by computing K*v for
# every point v on the sphere:
# ps = Ref(K) .* v

# Ultimately both of the above techniques only present part of the information
# about K in the position because K contains more degrees of freedom than we
# can represent on an ellipse.

x,y,z = split_coords(ps)
scene = surface(x,y,z, color=to_rgb.(v))

x,y,z = split_coords(1.003.*ps) # Scaling hack to avoid some z-fighting
wireframe!(scene, x,y,z)

# Makie.save("tensor.png", scene)

[KirillGerke]

Thank you a lot, Chris!
Your solution with Makie fully covers my original question.

My answer to your questions/suggestions is threefold:

1. You are totally right - this tensor is non-symmetrical. I explain it a bit later on, while describing my applications below. Originally i was using S and V, and i do not remember exactly why U was not considered. I understand it in the following way - we have eigenvalues of the tensor which point along some “major” direction. Now we plot the “rotation” of this vector using our original coordinate system.
   May be it will get more clear from the application…
2. My original application is to visualize tensorial physical properties of different materials. My first application (here is the proof of the first publication on the topic currently in press: https://drive.google.com/open?id=1_fCv_4DxmkHElNVxTUuAdV-v1d674r7S) is to visualize fluid permeability tensors. It is important to visually observe the difference in orientation and magnitude of different tensors, but simply looking at 3x3 matrix it is really hard to compare them or obtain some general idea of its properties.
   Permeability tensors are theoretically symmetrical. But in practice due to inaccuracies of numerical methods and some problems with boundary conditions resulting tensors are not symmetrical. I think we can fix that in future, current common solution is to average off-diagonal terms, e.g., xy and yx = (xy+yx)/2.
   I think there are plenty of physical properties for which tensors can be non-symmetrical… But we still want to visualize them of course.
   Hope this clarifies the aims, if not, please, let me know and i shall expand.
3. Makie seems to be very able workframe. Now i am wondering if it is able to produce visualizations of numerous tensor ellipsoids like the one you produced here, e.g., something like was shown here: Array of boxes using Makie and volume - #2 by lazarusA
   Each subcube of a domain, say 10x10x10 subcubes, contains its own tensor visualization.
   I guess it can be solved by looping through 3d array of tensors and plotting ellipsoids at different locations corresponding to each position within 3d array. The trick would be to scale the size of ellipsoids to the largest one to prevent any overlaps…

Many thanks for your help again! I shall definitely use your solution in my research!

[c42f]

Oh right, these are permeability tensors, thanks for explaining! In that case I think you could symmetrize the tensor for visualization purposes (if it mattered) and use the eigendecomposition if you wanted. Though as you can see, working directly with the tensor as a linear map can be very easy and arguably less mysterious.

Here’s an alternative visualization technique which doesn’t rely on the tensor being symmetrical. The idea is to plot the vector field K*v on the unit sphere using an arrow plot. I’ve added |K*v| as the color of the sphere.

I’ve also iterated over several origins to place a small 1x1x3 grid of these visualizations so you can see that’s quite easy.

using Makie
using LinearAlgebra
using StaticArrays

# Utility to split an array-of-length-3-arrays into x,y,z components
split_coords(ps) = getindex.(ps,1), getindex.(ps,2), getindex.(ps,3)

# Spherical Fibonacci numbers for Quasi-random distribution on the sphere.
# See "Spherical Fibonacci Point Sets for Illumination Integrals" by Marques et al.,
# doi:10.1111/cgf.12190
function spherical_fib(j, N)
    θ = acos(1 - 2*j/N)
    ϕ = 2*j*π/MathConstants.φ
    (θ,ϕ)
end

# Convert spherical coordinates to Cartesian coordinates.
# Note that θ,ϕ are named according to the "physicist convention"
# Note that using Vec3f0 here is efficient, but you could just use a
# normal array and everything would still work.
spherical_to_cartesian(r,θ,ϕ) = Vec3f0(r*sin(θ)*cos(ϕ), r*sin(θ)*sin(ϕ), r*cos(θ))

K = SA[ 0.070087    -0.0386396   0.00209088;
       -0.0441344    0.105735    0.00012994;
       -0.00553273  -0.00292807  0.0896236]
# If desired, symmetrize K
# K = 0.5*(K + K')

scene = Scene()
for origin in Point3f0[(-3,0,0), (0,0,0), (3,0,0)]
    # Plot a parametric unit sphere
    θ = range(0,pi,length=100)
    ϕ = permutedims(range(0,2pi,length=100))
    ps = spherical_to_cartesian.(1.0, θ, ϕ)
    surface!(scene, split_coords(Ref(origin) .+ ps)..., color=norm.(Ref(K) .* ps))

    N = 500
    vs = map(j->spherical_to_cartesian(1, spherical_fib(j,N)...), 1:N)
    Kvs = Ref(K) .* vs
    arrows!(scene, Ref(origin) .+ vs, 0.1 .* Kvs, arrowsize=0.05, linewidth=3, lengthscale=30, arrowcolor=norm.(Kvs))
end
scene

# Makie.save("tensor.png", scene)

[c42f]

To make things more interesting, here’s an application on random tensors K which are not symmetric. It’s particularly interesting for those which have complex eigenvalues (you can see one of these in particular has quite a rotational component).

@sdanisch you might be entertained by this picture (Though probably will notice something wrong I’ve done with the Makie API )

using Makie
using LinearAlgebra
using StaticArrays

# Utility to split an array-of-length-3-arrays into x,y,z components
split_coords(ps) = getindex.(ps,1), getindex.(ps,2), getindex.(ps,3)

# Spherical Fibonacci numbers for Quasi-random distribution on the sphere.
# See "Spherical Fibonacci Point Sets for Illumination Integrals" by Marques et al.,
# doi:10.1111/cgf.12190
function spherical_fib(j, N)
    θ = acos(1 - 2*j/N)
    ϕ = 2*j*π/MathConstants.φ
    (θ,ϕ)
end

# Convert spherical coordinates to Cartesian coordinates.
# Note that θ,ϕ are named according to the "physicist convention"
# Note that using Vec3f0 here is efficient, but you could just use a
# normal array and everything would still work.
spherical_to_cartesian(r,θ,ϕ) = Vec3f0(r*sin(θ)*cos(ϕ), r*sin(θ)*sin(ϕ), r*cos(θ))

scene = Scene()
for x in -3:3:3, y = -3:3:3
    origin = Point3f0(x,y,0)
    K = @SArray randn(3,3)
    # Plot a parametric unit sphere
    θ = range(0,pi,length=100)
    ϕ = permutedims(range(0,2pi,length=100))
    ps = spherical_to_cartesian.(1.0, θ, ϕ)
    surface!(scene, split_coords(Ref(origin) .+ ps)..., color=norm.(Ref(K) .* ps))

    N = 500
    vs = map(j->spherical_to_cartesian(1, spherical_fib(j,N)...), 1:N)
    Kvs = Ref(K) .* vs
    arrows!(scene, Ref(origin) .+ vs, 0.1 .* Kvs, arrowsize=0.05, linewidth=3, lengthscale=3, arrowcolor=norm.(Kvs))
end
scene

# Makie.save("tensor.png", scene)

[KirillGerke]

Chris,

To this i can say only - wow! It really looks cool!
I always produce sucky figs and schemes with PowerPoint style at best… Makie can help to reach next level)
I guess i have to read more about Makie. While it is definitely a cool package, i have to admit that its documentation and examples did not provide enough information to start working with it.
Hope it is OK to ask for more in future)

Again, thanks a lot, mate! I really like it!

[c42f]

Yeah, I am not an expert in the Makie API actually, I’m an occasional contributor . Generally I browse the image gallery in the docs and click on a few which seem like they might be similar to what I want.

[gilgannonj]

Hey @c42f this visualisation is really cool and I would love to use it for my permeability tensors! I tried to run it but I can’t get it working. I was wondering what version of Julia and the other libs are that you used?

I am running Julia v1.6 and Makie v0.17.8.

I am new to Julia (Python is normally my jam) so I am still trying to get my head around the issues but it seems that when I run it I have problems with Vec3f0() and Point3f0() not being declared. I can’t find documentation on it but I assume that these are standard library functions so I am not sure how I am supposed to declare them.

[AndiMD]

Hi @gilgannonj
I tried to update the example for you, using GLMakie v0.6.9 and Julia 1.8-rc3 (use mouse to rotate, mouse wheel to zoom)

using GLMakie
GLMakie.activate!()
using LinearAlgebra
using StaticArrays

# Utility to split an array-of-length-3-arrays into x,y,z components
split_coords(ps) = getindex.(ps,1), getindex.(ps,2), getindex.(ps,3)

# Spherical Fibonacci numbers for Quasi-random distribution on the sphere.
# See "Spherical Fibonacci Point Sets for Illumination Integrals" by Marques et al.,
# doi:10.1111/cgf.12190
function spherical_fib(j, N)
    θ = acos(1 - 2*j/N)
    ϕ = 2*j*π/MathConstants.φ
    (θ,ϕ)
end

# Convert spherical coordinates to Cartesian coordinates.
# Note that θ,ϕ are named according to the "physicist convention"
# Note that using Vec3f0 here is efficient, but you could just use a
# normal array and everything would still work.
spherical_to_cartesian(r,θ,ϕ) = Vec3f(r*sin(θ)*cos(ϕ), r*sin(θ)*sin(ϕ), r*cos(θ))

scene = Scene()
cam = Camera3D(scene; projectiontype =:Perspective)
for x in -3:3:3, y = -3:3:3
    origin = Point3f(x,y,0)
    K = @SArray randn(3,3)
    # Plot a parametric unit sphere
    θ = range(0,pi,length=100)
    ϕ = permutedims(range(0,2pi,length=100))
    ps = spherical_to_cartesian.(1.0, θ, ϕ)
    surface!(scene, split_coords(Ref(origin) .+ ps)..., color=norm.(Ref(K) .* ps))

    N = 500
    vs = map(j->spherical_to_cartesian(1, spherical_fib(j,N)...), 1:N)
    Kvs = Ref(K) .* vs
    arrows!(scene, Ref(origin) .+ vs, 0.1 .* Kvs, arrowsize=0.05, linewidth=0.02, lengthscale=1, arrowcolor=norm.(Kvs))
end
scene

# Makie.save("tensor.png", scene)