3

我有一个 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]

在这里，我假设惯性张量的特征向量定义了旋转矩阵（基于这个问题和谷歌结果，比如这个网页似乎是正确的？）但这并没有给我正确的结果。

我希望打印的矩阵是：

[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]

我不确定问题是我的代码还是我对旋转主轴的假设，但我们将不胜感激。

4

1 回答 1

0

这是我为进行这种对齐而开发的代码的链接。

给定一组具有坐标 (x,y,z) 的散点，目标是将与最小特征值相关联的特征向量与 3D 笛卡尔轴的 X 轴相匹配，并将与中值特征值相关联的特征向量与 Y 轴相匹配来自相同的 3D 笛卡尔轴。

为此，我遵循了以下步骤：

  1. 仅通过执行以下操作将质心在 (xmn, ymn, zmn) 的点集转换为质心在 (0,0,0) 的新点集：(x-xmn, y-ymn, z-zmn)。

  2. 计算与最小特征值 ( min_eigen ) 相关的特征向量的 xy 投影与笛卡尔轴中的 x 轴之间的角度 THETA（绕 z 旋转） 。获得结果 tetha 后，将 min_eigen 旋转给定的 theta，使其包含在 xy 平面中。我们称这个结果向量为：rotz

  3. 计算rotz和 x 轴之间的角度 PHI 以执行围绕 y 轴的旋转。获得 phi 后，对rotz应用绕y 轴的旋转。通过最后一次旋转，与中间特征向量 (medium_eigen) 关联的特征向量位于笛卡尔轴的 yz 投影中，因此我们只需要找到 medium_eigen 与笛卡尔轴的 y 轴之间的角度。

  4. 计算 medium_eigen 和 y 轴之间的角度 ALPHA。应用围绕 x 轴的旋转 aaaand：它完成了！

注意：将步骤 1、2、3 应用于您的点集后，您必须重新计算 3D_SVD（3D_single 值分解）并根据结果特征向量集，然后使用新的 medium_eigen 实施第四步。

我真的希望这会有所帮助。

旋转是通过此处定义的旋转矩阵实现的：Rotating a Vector in 3D Space

于 2019-04-26T12:17:17.997 回答