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

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
verticesArray[1,counter] = world_coord
verticesArray[2,counter] = world_coord
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?

• the object vertices are probably stored in the local coordinate system of the mesh, so you'll either need to create or specify a transformation matrix to retreive the world coordinates them. i don't have a lot of python experience, but I did run into this issue in AutoCAD quite often – John Burrill Aug 1 '14 at 17:36
• @JohnBurrill: I believed that is what I am doing by multiplying the coordinates with "matrix_world", but there seems to be sth wrong with that. Maybe someone knows what and which matrix I need to use... – user1809923 Aug 4 '14 at 7:00

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
verticesArray[1,counter] = world_coord
verticesArray[2,counter] = world_coord
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, vec, vec))

del world_coord
del matrix_world
print(eigenvectors)