How to align principal axes of 3D density map in numpy with Cartesian axes?

Question

I have an n x n x n numpy array that contains density values on a cubic grid. I'm trying to align the principal axes of inertia of the density map with the cartesian x,y,z axes of the grid. I have the following so far:

import numpy as np
from scipy import ndimage

def center_rho(rho):
    """Move density map so its center of mass aligns with the center of the grid"""
    rhocom = np.array(ndimage.measurements.center_of_mass(rho))
    gridcenter = np.array(rho.shape)/2.
    shift = gridcenter-rhocom
    rho = ndimage.interpolation.shift(rho,shift,order=1,mode='wrap')
    return rho

def inertia_tensor(rho,side):
    """Calculate the moment of inertia tensor for the given density map."""
    halfside = side/2.
    n = rho.shape[0]
    x_ = np.linspace(-halfside,halfside,n)
    x,y,z = np.meshgrid(x_,x_,x_,indexing='ij')
    Ixx = np.sum(rho*(y**2 + z**2))
    Iyy = np.sum(rho*(x**2 + z**2))
    Izz = np.sum(rho*(x**2 + y**2))
    Ixy = -np.sum(rho*x*y)
    Iyz = -np.sum(rho*y*z)
    Ixz = -np.sum(rho*x*z)
    I = np.array([[Ixx, Ixy, Ixz],
                  [Ixy, Iyy, Iyz],
                  [Ixz, Iyz, Izz]])
    return I

def principal_axes(I):
    """Calculate the principal inertia axes and order them in ascending order."""
    w,v = np.linalg.eigh(I)
    return w,v

#number of grid points along side
n = 10
#note n <= 3 produces unit eigenvectors, not sure why
#in practice, n typically between 10 and 50
np.random.seed(1)
rho = np.random.random(size=(n,n,n))
side = 1. #physical width of box, set to 1.0 for simplicity

rho = center_rho(rho)
I = inertia_tensor(rho,side)
PAw, PAv = principal_axes(I)

#print magnitude and direction of principal axes
print "Eigenvalues/eigenvectors before rotation:"
for i in range(3):
    print PAw[i], PAv[:,i]

#sanity check that I = R * D * R.T
#where R is the rotation matrix and D is the diagonalized matrix of eigenvalues
D = np.eye(3)*PAw
print np.allclose(np.dot(PAv,np.dot(D,PAv.T)),I)

#rotate rho to align principal axes with cartesian axes
newrho = ndimage.interpolation.affine_transform(rho,PAv.T,order=1,mode='wrap')

#recalculate principal axes
newI = inertia_tensor(newrho,side)
newPAw, newPAv = principal_axes(newI)

#print magnitude and direction of new principal axes
print "Eigenvalues/eigenvectors before rotation:"
for i in range(3):
    print newPAw[i], newPAv[:,i]

Here I'm assuming that the eigenvectors of the inertia tensor define the rotation matrix (which based on this question and Google results such as this webpage seems correct?) However this doesn't give me the correct result.

I expect the printed matrix to be:

[1 0 0]
[0 1 0]
[0 0 1]

(which could be wrong) but don't even get unit vectors to start with. What I get is:

Eigenvalues/eigenvectors before rotation:
102.405523732 [-0.05954221 -0.8616362   0.5040216 ]
103.177395578 [-0.30020273  0.49699978  0.81416801]
104.175688943 [-0.95201526 -0.10283129 -0.288258  ]
True
Eigenvalues/eigenvectors after rotation:
104.414931478 [ 0.38786    -0.90425086  0.17859172]
104.731536038 [-0.74968553 -0.19676735  0.63186566]
106.151322662 [-0.53622405 -0.37896304 -0.75422197]

I'm not sure if the problem is my code or my assumptions about rotating principal axes, but any help would be appreciated.

I know this is old, but I'd like to ping this to see if someone can help out. — tomerg, Nov 15 '17 at 16:47
It turns out this does work, you just have to run it in a loop a few times. The interpolation causes issues. Once its applied 3 times or so it does give the identity matrix to within an acceptable tolerance. — tomerg, May 24 '19 at 16:34

pedro_galher · Answer 1 · 2019-04-26T16:33:32.493

Here is the link to the code I developed to do such alignment.

Given a set of scatter points with coordinates (x,y,z), the objective is to match the eigenvector associated to the minimum eigenvalue with the X-axis of a 3D cartesian axis and the eigenvector associated to the median eigenvalue with the Y axis from the same 3D cartesian axis.

For this purpose, I followed the following steps:

Translate the set of points with centroid in (xmn, ymn, zmn) to a new set of points with centroid in (0,0,0) only by doing: (x-xmn, y-ymn, z-zmn).
Calculate the angle THETA (rotation around z) between the xy-projection of the eigenvector associated to the minimum eigenvalue (min_eigen) and the x-axis in a cartesian axis. After obtention of the resulting tetha, rotate the min_eigen the given theta so that it is contained in the xy-plane. Let's call this resulting vector: rotz
Calculate the angle PHI between rotz and x-axis in order to perform a rotation around the y-axis. Once the phi is obtained, a rotation is applied to rotz aound the y axis. With this last rotation, the eigenvector associated to the medium eigenvector (medium_eigen) is then in the yz proyection of the cartesian axis, so we will just need to find the angle between medium_eigen and the y-axis of the cartesian axis.
Calculate the angle ALPHA between the medium_eigen and y-axis. Apply the rotation around the x-axis aaaand: IT'S DONE!

NOTE: After applying steps 1,2,3 to your set of points, you have to recalculate the 3D_SVD (3D_single value decomposition) and from the resulting set of eigenvectors, then implement the 4th step with the new medium_eigen.

I really hope this helps.

The rotations are implemented by means of the rotation matrix defined here: Rotating a Vector in 3D Space

Thanks for the effort. I'm actually looking for how to rotate an N x N x N array though, not a set of N points (which is an N x 3 array). The issue primarily is interpolation I think. I've now got some code to solve this however. — tomerg, May 24 '19 at 16:33

How to align principal axes of 3D density map in numpy with Cartesian axes?

1 Answers1

Linked