title: How to align principal axes of 3D density map in numpy with Cartesian axes?
body: <p>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:</p>
  
  <pre><code>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 &lt;= 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]
  </code></pre>
  
  <p>Here I'm assuming that the eigenvectors of the inertia tensor define the rotation matrix (which based on <a href="https://math.stackexchange.com/questions/145023/rotation-matrix-from-an-inertia-tensor">this question</a> and Google results such as <a href="http://www.euclideanspace.com/physics/dynamics/inertia/rotation/index.htm" rel="nofollow noreferrer">this webpage</a> seems correct?) However this doesn't give me the correct result.</p>
  
  <p>I expect the printed matrix to be:</p>
  
  <pre><code>[1 0 0]
  [0 1 0]
  [0 0 1]
  </code></pre>
  
  <p>(which could be wrong) but don't even get unit vectors to start with. What I get is:</p>
  
  <pre><code>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]
  </code></pre>
  
  <p>I'm not sure if the problem is my code or my assumptions about rotating principal axes, but any help would be appreciated.</p>