python三维坐标轴单位长度_如何将numpy三维密度图的主轴与笛卡尔坐标轴对齐？...

我有一个nxnxnnumpy数组，它包含立方体网格上的密度值。我试着把密度图的惯性主轴和网格的笛卡尔坐标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]

这里我假设惯性张量的特征向量定义了旋转矩阵(基于this question和谷歌搜索结果，比如this webpage似乎是正确的？)但是这个结果不正确。在

我希望印刷矩阵是：

^{pr2}$

(这可能是错误的)但是不要从单位向量开始。我得到的是：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]

我不确定问题是我的代码还是我关于旋转主轴的假设，但是如果有任何帮助，我将不胜感激。在