Как выровнять главные оси трехмерной карты плотности в NumPy с декартовыми осями?

Question

У меня есть массив nxnxn numpy, который содержит значения плотности в кубической сетке. Я пытаюсь совместить главные оси инерции карты плотности с декартовыми осями x,y,z сетки. Пока у меня есть следующее:

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]

Здесь я предполагаю, что собственные векторы тензора инерции определяют матрицу вращения (которая основана на этом вопросе, и результаты Google, такие как эта веб-страница, кажутся правильными?) Однако это не дает мне правильный результат.

Я ожидаю, что напечатанная матрица будет:

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

(что может быть неправильно), но даже не начинайте с единичных векторов. Что я получаю это:

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]

Я не уверен, является ли проблема моим кодом или моими предположениями относительно вращения главных осей, но любая помощь будет принята с благодарностью.

Answer 1

Вот ссылка на код, который я разработал для такого выравнивания.

При заданном наборе точек рассеяния с координатами (x,y,z) цель состоит в том, чтобы сопоставить собственный вектор, связанный с минимальным собственным значением, с осью X трехмерной декартовой оси и собственным вектором, связанным с медианным собственным значением с осью Y с той же трехмерной декартовой оси.

Для этого я выполнил следующие шаги:

1. Переведите набор точек с центроидом в (xmn, ymn, zmn) в новый набор точек с центроидом в (0,0,0), выполнив: (x-xmn, y-ymn, z-zmn).

2. Вычислить угол THETA (вращение вокруг z) между проекцией xy собственного вектора, связанного с минимальным собственным значением (min_eigen), и осью x на декартовой оси. Получив полученную тету, поверните min_eigen заданной тэты, чтобы она находилась в плоскости xy. Давайте назовем этот результирующий вектор: rotz

3. Рассчитайте угол PHI между осью X и X, чтобы выполнить вращение вокруг оси Y. Как только фи получен, вращение применяется вокруг оси у. При этом последнем вращении собственный вектор, связанный с собственным вектором среды (medium_eigen), находится затем в проекции yz декартовой оси, поэтому нам просто нужно найти угол между medium_eigen и осью y декартовой оси.

4. Рассчитайте угол ALPHA между medium_eigen и осью y. Примените вращение вокруг оси Xaaa: ЭТО СДЕЛАНО!

ПРИМЕЧАНИЕ. После применения шагов 1,2,3 к вашему набору точек вы должны пересчитать 3D_SVD (разложение по значениям 3D_single) и из полученного набора собственных векторов, а затем выполнить 4-й шаг с новым medium_eigen.

Я действительно надеюсь, что это поможет.

Вращения осуществляются с помощью матрицы вращения, определенной здесь: Вращение вектора в трехмерном пространстве