Construct Ellipsoid Given Inertia Tensor

[daniel-a-diaz]

So my research is in Additive Manufacturing and I work with 3D X-ray CT scans that capture gas pores left over from processing. To extract the pores and calculate physical descriptors I use scikit-image's measure module to label the segmented CT scan and calculate region properties of the labeled pores. With measure.regionprops() I can use inertia_tensor to calculate things like aspect ratio and orientation with respect to the xy plane. What I'm trying to do now is figure out how to construct equivalent ellipsoids for each extracted pore using the inertia tensor. With its eigenvalues I can get the axis lengths and aspect ratio, and the eigenvectors give me the orientation, but I'm not sure how to translate that into the axis parameters for vedo's Ellipsoid. Below is an example script to calculate these quantities and render a visualization of the eigenvectors. How do I get to an Ellipsoid from this information?

from skimage import measure
from skimage.measure._regionprops import _inertia_eigvals_to_axes_lengths_3D as axis_lengths_3D
from vedo import Arrow, Ellipsoid, Mesh, Volume, show
import numpy as np

# Create Ellipsoid
ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(2,2,0), axis3=(3,0,3))
mesh = Mesh(ellipsoid).normalize()
vol = mesh.binarize(spacing=(0.03,0.03,0.03))
np_vol = vol.tonumpy()

# Get Region Properties
labels = measure.label(np_vol, connectivity=2)
properties = measure.regionprops(labels)

centroid = properties[0].centroid
inertia_tensor = properties[0].inertia_tensor
inertia_tensor_eigvals = properties[0].inertia_tensor_eigvals

# Calculate Aspect Ratio
axis_lengths = np.array(axis_lengths_3D(inertia_tensor_eigvals))
aspect_ratio = axis_lengths[0] / axis_lengths[-1]

# Calculate Orientation
eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor)
eigenvalues = np.clip(eigenvalues, 0, None, out=eigenvalues)
sorted_indices = np.argsort(eigenvalues)
eigenvectors = np.take_along_axis(eigenvectors, np.stack([sorted_indices]*3), axis=1)

eigenvector_0 = eigenvectors[:,0]
eigenvector_1 = eigenvectors[:,1]
eigenvector_2 = eigenvectors[:,2]

theta = np.arctan(np.abs(eigenvector_0[2]/np.sqrt(eigenvector_0[0]**2+eigenvector_0[1]**2)))*180/np.pi

# Visualize Eigenvectors
vol = Volume(np_vol)
vol.alpha((0, 0.05))

vector_0 = Arrow(start_pt=centroid, end_pt=centroid+40*eigenvector_0, c='red')
vector_1 = Arrow(start_pt=centroid, end_pt=centroid+30*eigenvector_1, c='blue')
vector_2 = Arrow(start_pt=centroid, end_pt=centroid+20*eigenvector_2, c='green')

print(f"Aspect Ratio = {aspect_ratio:.3f}, Orientation = {theta:.3f} degrees")
show((vol, vector_0, vector_1, vector_2), axes=1, bg2='lightblue', interactive=True, viewup='z', new=True)

Aspect Ratio = 4.253, Orientation = 30.025 degrees

[marcomusy]

Hi Daniel, what about this solution:

from vedo import *

# easier to inspect things without perspective:
settings.use_parallel_projection = True 

t = LinearTransform()
ev = [4,3,2] # some fake eigenvalues and eigenvectors:
t.scale(ev).rotate_z(45).rotate_y(20).rotate_z(20)

# position of the ellipsoid and arrows
p = [11,12,13]

elli = Sphere().apply_transform(t)
elli.color('green').alpha(0.2)

a1, a2, a3 = t.matrix3x3.T
arr1 = Arrow([0,0,0], a1).c('red5')
arr2 = Arrow([0,0,0], a2).c('green5')
arr3 = Arrow([0,0,0], a3).c('blue5')

# we can group the 4 objects into a single one for convenience
tens = Assembly(arr1, arr2, arr3, elli).pos(p)
show(tens, axes=1)

You can create the linear transform directly with your 3x3 matrix with t = LinearTransform(my_matrix).
(always make some tests to double check that you get what you expect to see!).

[daniel-a-diaz]

So it's almost there. I've tried using the inertia tensor and the eigenvector tensor as the matrix in the LinearTransform, and while it looks like it gets the correct proportions of everything, the rotations are off. I've also tried transposing them in case the convention is different. And it creates a gigantic Ellipsoid so I had to scale it down considerably to get it in the same plot as the original Ellipsoid. Below is the code where I added the Ellipsoid created with LinearTransform.

from skimage import measure
from skimage.measure._regionprops import _inertia_eigvals_to_axes_lengths_3D as axis_lengths_3D
from vedo import *
import numpy as np

settings.use_parallel_projection = True 

# Create Ellipsoid
ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(2,2,0), axis3=(3,0,3))
#mesh = Mesh(ellipsoid).normalize()
vol = ellipsoid.binarize(spacing=(0.05,0.05,0.05))
np_vol = vol.tonumpy()

# Get Region Properties
labels = measure.label(np_vol, connectivity=2)
properties = measure.regionprops(labels)

centroid = properties[0].centroid
inertia_tensor = properties[0].inertia_tensor
inertia_tensor_eigvals = properties[0].inertia_tensor_eigvals

# Calculate Aspect Ratio
axis_lengths = np.array(axis_lengths_3D(inertia_tensor_eigvals))
aspect_ratio = axis_lengths[0] / axis_lengths[-1]

# Calculate Orientation
eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor)
eigenvalues = np.clip(eigenvalues, 0, None, out=eigenvalues)
sorted_indices = np.argsort(eigenvalues)
eigenvectors = np.take_along_axis(eigenvectors, np.stack([sorted_indices]*3), axis=1)

eigenvector_0 = eigenvectors[:,0]
eigenvector_1 = eigenvectors[:,1]
eigenvector_2 = eigenvectors[:,2]

theta = np.arctan(np.abs(eigenvector_0[2]/np.sqrt(eigenvector_0[0]**2+eigenvector_0[1]**2)))*180/np.pi

# Visualize Eigenvectors
vol = Volume(np_vol)
vol.alpha((0, 0.05))

vector_0 = Arrow(start_pt=centroid, end_pt=centroid+40*eigenvector_0, c='red')
vector_1 = Arrow(start_pt=centroid, end_pt=centroid+30*eigenvector_1, c='blue')
vector_2 = Arrow(start_pt=centroid, end_pt=centroid+20*eigenvector_2, c='green')

# Construct Ellipsoid With Inertia Tensor
transform = LinearTransform(inertia_tensor).scale(0.08)
ell = Sphere().apply_transform(transform).pos(centroid).color('green').alpha(0.5)

assembly = Assembly(vol, vector_0, vector_1, vector_2, ell)

print(f"Aspect Ratio = {aspect_ratio:.3f}, Orientation = {theta:.3f} degrees")
show(assembly, axes=1, bg2='lightblue', interactive=True, viewup='z', new=True)

[marcomusy]

It lloks to me that there must be some mismatch in xz convention + maybe sorting (?) the eigenvalues

from skimage import measure
from vedo import *

settings.use_parallel_projection = True 

# Create Ellipsoid
ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(0,2,0), axis3=(0,0,3))
vol = ellipsoid.binarize(spacing=(0.1,0.1,0.1)).print()
np_vol = vol.tonumpy()
# np_vol = np.transpose(np_vol, axes=[2, 1, 0])

# Get Region Properties
labels = measure.label(np_vol, connectivity=2)
properties = measure.regionprops(labels)

centroid = properties[0].centroid
inertia_tensor = properties[0].inertia_tensor * 0.1
inertia_tensor_eigvals = properties[0].inertia_tensor_eigvals

# Calculate Orientation
eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor)

# eigenvalues = np.clip(eigenvalues, 0, None, out=eigenvalues)
# sorted_indices = np.argsort(eigenvalues)
# eigenvectors = np.take_along_axis(eigenvectors, np.stack([sorted_indices]*3), axis=1)
# eigenvector_0 = eigenvectors[:,0]*eigenvalues[sorted_indices[0]]
# eigenvector_1 = eigenvectors[:,1]*eigenvalues[sorted_indices[1]]
# eigenvector_2 = eigenvectors[:,2]*eigenvalues[sorted_indices[2]]

eigenvector_0 = eigenvectors[:,0]*eigenvalues[0]
eigenvector_1 = eigenvectors[:,1]*eigenvalues[1]
eigenvector_2 = eigenvectors[:,2]*eigenvalues[2]
print("eigenvectors\n", eigenvectors)
print("eigenvalues ", eigenvalues)

# Visualize Eigenvectors
vol = Volume(np_vol).alpha((0, 0.05))

vector_0 = Arrow(start_pt=centroid, end_pt=centroid+eigenvector_0, c='red')
vector_1 = Arrow(start_pt=centroid, end_pt=centroid+eigenvector_1, c='blue')
vector_2 = Arrow(start_pt=centroid, end_pt=centroid+eigenvector_2, c='green')

# Construct Ellipsoid With Inertia Tensor
transform = LinearTransform(inertia_tensor).scale(1).translate(centroid)

ell = Sphere().apply_transform(transform).color('green').alpha(0.5)
assembly = Assembly(vol, vector_0, vector_1, vector_2, ell)

print("inertia_tensor\n", inertia_tensor)
print(transform)
show(assembly, 
    ellipsoid.scale(10).pos(centroid).wireframe(),
    axes=1, bg2='lightblue', viewup='z')

i will only have time to look into this tuesday evening.. I also found a bug in the Ellipsoid class (but this should not affect the above),
thanks for reporting!

[marcomusy]

I had some time to play around and indeed ti seems that this inertia_tensor is computed under some different orientation/convention.
I tried with np.transpose(np_vol, axes=[2, 1, 0]) but could not guess the right combination..
A possible way around this is to rotate around y:

from skimage import measure
from vedo import *

settings.use_parallel_projection = True 

# Create Ellipsoid
ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(0,2,0), axis3=(0,0,3))
ellipsoid.rotate_z(45).rotate_x(20).pos(30,40,50)

vol = ellipsoid.binarize(spacing=(0.1,0.1,0.1)).print()
np_vol = vol.tonumpy()
# np_vol = np.transpose(np_vol, axes=[2, 1, 0])

labels = measure.label(np_vol, connectivity=2)
properties = measure.regionprops(labels)

centroid = properties[0].centroid
inertia_tensor = properties[0].inertia_tensor * 0.1
eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor)
e0 = eigenvectors[:,0]*eigenvalues[0]
e1 = eigenvectors[:,1]*eigenvalues[1]
e2 = eigenvectors[:,2]*eigenvalues[2]
# print("eigenvectors\n", eigenvectors)
# print("eigenvalues ", eigenvalues)

# Construct Ellipsoid With Inertia Tensor
transform = LinearTransform(inertia_tensor)
transform.rotate(90, e1)
transform.shift(centroid)

tens = Sphere(res=24).apply_transform(transform)
tens.color('w').alpha(0.25).lw(1)

line0 = Line([0,0,0], [1,0,0]).apply_transform(transform).color('r').lw(4)
line1 = Line([0,0,0], [0,1,0]).apply_transform(transform).color('g').lw(4)
line2 = Line([0,0,0], [0,0,1]).apply_transform(transform).color('b').lw(4)

evector0 = Arrow(end_pt=e0, c='r').rotate(90, e1).shift(centroid)
evector1 = Arrow(end_pt=e1, c='g').rotate(90, e1).shift(centroid)
evector2 = Arrow(end_pt=e2, c='b').rotate(90, e1).shift(centroid)

triads = Assembly(line0, line1, line2, evector0, evector1, evector2)

print(transform)
print("inertia_tensor\n", inertia_tensor)
show(
    # ellipsoid.scale(10).pos(centroid).wireframe(), # should roughly match the tensor
    tens,
    triads,
    axes=1, bg2='lightblue', viewup='z',
)

[daniel-a-diaz]

Ok, I was able to take what you had and consistently match it to the same orientation of the original numpy ellipsoid by sorting the eigenvectors first, and then rotating 90 degrees about the eigenvector axis with the middle eigenvalue. The one thing I can't figure out though is how to match the scale of the original numpy ellipsoid. I scaled it by 0.08 to match this specific case, but the scale needed to match to the original is dependent on its size.

from skimage import measure
from skimage.measure._regionprops import _inertia_eigvals_to_axes_lengths_3D as axis_lengths_3D
from vedo import *
import numpy as np

settings.use_parallel_projection = True 

# Create Ellipsoid
ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(2,2,0), axis3=(3,0,3))
vol = ellipsoid.binarize(spacing=(0.05,0.05,0.05))
np_vol = vol.tonumpy()

# Get Region Properties
labels = measure.label(np_vol, connectivity=2)
properties = measure.regionprops(labels)

centroid = properties[0].centroid
inertia_tensor = properties[0].inertia_tensor

# Calculate and Sort Eigenvalues and Eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor)
eigenvalues = np.clip(eigenvalues, 0, None, out=eigenvalues)
sorted_indices = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sorted_indices]
eigenvectors = np.take_along_axis(eigenvectors, np.stack([sorted_indices]*3), axis=1)

# Calculate Aspect Ratio
axis_lengths = np.array(axis_lengths_3D(eigenvalues))
aspect_ratio = axis_lengths[0] / axis_lengths[-1]

# Calculate Orientation
eigenvector0 = eigenvectors[:,0]*eigenvalues[0]
eigenvector1 = eigenvectors[:,1]*eigenvalues[1]
eigenvector2 = eigenvectors[:,2]*eigenvalues[2]

theta = np.arctan(np.abs(eigenvector0[2]/np.sqrt(eigenvector0[0]**2+eigenvector0[1]**2)))*180/np.pi

# Construct Ellipsoid With Inertia Tensor
scale = 0.08
transform = LinearTransform(inertia_tensor).rotate(90, eigenvector1).shift(centroid).scale(scale)

ellipsoid = Sphere(res=24).apply_transform(transform)
ellipsoid.color('w').alpha(0.25).lw(1)

evector0 = Arrow(end_pt=eigenvector0, c='r').rotate(90, eigenvector1).shift(centroid).scale(scale)
evector1 = Arrow(end_pt=eigenvector1, c='g').rotate(90, eigenvector1).shift(centroid).scale(scale)
evector2 = Arrow(end_pt=eigenvector2, c='b').rotate(90, eigenvector1).shift(centroid).scale(scale)

evectors = Assembly(evector0, evector1, evector2)

orig_vol = Volume(np_vol).alpha((0, 0.05))

print(f"Aspect Ratio = {aspect_ratio:.3f}, Orientation = {theta:.3f} degrees")
show(
    ellipsoid,
    evectors,
    orig_vol,
    axes=1, bg2='lightblue', viewup='z',
)

[daniel-a-diaz]

@marcomusy thanks so much for all the help. That LinearTransform(inertia_tensor) solution was exactly what I needed. And thanks for the great tool. It makes dealing with 3D data so much easier.

[marcomusy]

Thant happens because the numpy array has no information about positioning in space, but vedo has it!
so when you create it from a numpy array you need to specify the voxel size (spacing) and positioning (origin) otherwise they are assumed to be [1,1,1], and [0,0,0].
Note that spacing can be different in xyz in vedo but skimage.image cannot cope with that (i guess).

Please do

pip install -U git+https://github.com/marcomusy/vedo.git

then

from vedo import *

settings.use_parallel_projection = True 

# Create Ellipsoid
# ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(2,2,0), axis3=(3,0,3)) # why??
ellipsoid = Ellipsoid(axis1=(1,0,0), axis2=(0,2,0), axis3=(0,0,3))
ellipsoid.rotate_z(45).rotate_x(20).pos([300,400,500])

s = 0.05
vol = ellipsoid.binarize(spacing=(s,s,s)) # keep equal spacing
vol.print()
origin = vol.origin()
spacing = vol.spacing()

np_vol = vol.tonumpy()

redo_vol = Volume(np_vol)
redo_vol.spacing(spacing).origin(origin)

show(vol, redo_vol, N=2, axes=1, bg2='lightgreen', viewup='z').close()

# ...
# inertia_tensor = properties[0].inertia_tensor * s

[daniel-a-diaz]

Where I'm having issues with size is in the equivalent ellipsoid created from the inertia tensor LinearTransform. I also need to make into a numpy array for future use, and it needs to match the size of the original numpy array the inertia tensor was based on. So for instance the shape of the original numpy array is (53,73,53), so after I perform the LinearTransform(inertia_tensor).rotate(90, eigenvector1).shift(centroid) on Sphere() and make that into a numpy array, it needs to end up (53,73,53).

I think I found a solution though. After I perform all the transforms in vedo to get to the equivalent ellipsoid numpy array, I then use skimage.transform.resize(ellipsoid, (53,73,53)) to match the (53,73,53) shape of the original numpy array.