Calculating the covariance matrix of an object's points - Rotation of an object

Question

I am trying to calculate the axes of largest variances of an object's points. As you might know, the eigenvectors of the covariance matrix of the points should yield the result.

For a cube (which is made out of 8 evenly spread points) I would therefore expect to get the three axes of the cube. However, this is only the case when I did not rotate the cube, but when I do, my method yields axes that are not correct.

Here is some code I use to calculate the axes:

import bpy
import numpy
from numpy import linalg as LA
vertices = bpy.context.active_object.data.vertices
matrix_world = bpy.context.active_object.matrix_world
verticesArray = numpy.zeros((3,len(vertices)))
counter = 0
for v in vertices:
     world_coord = matrix_world * v.co
     verticesArray[0,counter]= world_coord[0]
     verticesArray[1,counter] = world_coord[1]
     verticesArray[2,counter] = world_coord[2]
     counter = counter + 1
convMat = numpy.cov(verticesArray)
eigenvalues, eigenvectors = LA.eig(convMat)       # the eigenvectors determine the axes I want

As you see, I firstly transform all vertices to world coordinates. I guess somewhere there has to be an error that does not account for the rotation or something... Do you have an idea what I am missing?

Answer 1

I suggest using the addon Math vis to display vectors, in addition to testing your code in the console.

Paste the following in the console to see that the vertices are placed at the right locations. I don't know enough maths, but maybe the problem is in the numpy operations...

import bpy
import numpy

#from mathutils import Vector

from numpy import linalg as LA
vertices = bpy.context.active_object.data.vertices
matrix_world = bpy.context.active_object.matrix_world
verticesArray = numpy.zeros((3,len(vertices)))
counter = 0
for index, v in enumerate(vertices):
    globals()["vertex"+str(index)] = matrix_world * v.co ##vector to global variable
    world_coord = matrix_world * v.co
    verticesArray[0,counter]= world_coord[0]
    verticesArray[1,counter] = world_coord[1]
    verticesArray[2,counter] = world_coord[2]
    counter = counter + 1

convMat = numpy.cov(verticesArray)
eigenvalues, eigenvectors = LA.eig(convMat)
for index, vec in enumerate(eigenvectors):
    print(index, vec)
    globals()["vector"+str(index)] = Vector((vec[0], vec[1], vec[2]))

del world_coord
del matrix_world
print(eigenvectors)