我有一个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和Google结果，比如this webpage似乎是正确的？）但是这并没有给我正确的结果。
我希望打印的矩阵是：
[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]
我不确定问题是我的代码还是我对旋转主轴的假设，但如果有任何帮助，我将不胜感激。
最佳答案
Here是我为实现这种对齐而开发的代码的链接。
给定一组具有坐标（x，y，z）的散乱点，目标是将与最小特征值相关联的特征向量与三维笛卡尔坐标轴的x轴和与来自同一三维笛卡尔轴的y轴相关的中值特征值相关的特征向量相匹配。
为此，我遵循以下步骤：
仅通过执行以下操作，将质心位于（xmn，ymn，zmn）的点集转换为质心位于（0,0,0）的新点集：（x-xmn，y-ymn，z-zmn）。
计算与最小特征值（最小特征值）和X轴在笛卡尔轴中相关联的特征向量的XY投影之间的角度θ（z的旋转）。在获得所得到的系链之后，旋转给定的θ的最小本征值，使其包含在xy平面中。我们把这个结果向量叫做rotz
计算y轴和x轴之间的夹角φ，以便在y轴附近进行旋转。一旦获得PHI，旋转被施加到绕Y轴的ROZ。在最后一个旋转中，与介质特征向量相关的特征向量（MeimuMuthigeIn）是在笛卡尔轴的YZ投影中，所以我们只需要找到中间特征和笛卡尔轴的Y轴之间的夹角。
计算中间特征线和Y轴之间的角度α。在X轴AAAX周围应用旋转：完成了！
注意：将步骤1、2、3应用于点集后，必须重新计算3D_SVD（3D_单值分解）和由此得到的特征向量集，然后使用新的介质_特征来实现第4步。
我真的希望这能有帮助。
旋转通过这里定义的旋转矩阵来实现：Rotating a Vector in 3D Space
关于python - 如何将numpy中的3D密度贴图的主轴与笛卡尔坐标轴对齐？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/37006542/