Update: the problem has been solved. The solution can be found further down in this thread.
Original message below the line
For nearly a week now I’ve been stuck on a relatively simple problem: calculating a best fit plane to a number of vertices and then projecting the vertices on it.
What I got so far
First I calculate the center of mass of all vertices. This is simply the average of all locations of the vertices.
After that I calculate a 3x3 matrix (a covariance matrix):
where:
li = X - averageX
mi = Y - averageY
ni = Z - averageZ
Using Cardano’s method I find the 3 eigenvalues that belong to this matrix (lambda1 = lambda2 < lambda3).
A new matrix is calculated: M = T-lambda3*I
where I is an identity matrix. So basically M is just the first matrix where lambda is subtracted from each of the diagonal values.
The eigenvector that belongs to lambda3 is perpendicular to the two linearly independent rows in the matrix M. So it can be calculated by the dotproduct between the two rows.
The final step is projecting the vertices on the plane defined by the eigenvector (which should be the normal of that plane).
The problem
The above works fine, except from one thing. The eigenvector that is calculated isn’t the normal of the plane, but is a vector along that plane. Now everywhere I look says that this method should give you the normal, so I’ve probably made a mistake somewhere, but I really can’t find it.
Any help would be appreciated. Below you will find an image to illustrate the problem and my current code.
As you can see, the vector calculated by the script is perfectly located on the best fit plane, but is not the actual normal.
import bpy
import Blender
from Blender import *
# calculates the eigenvalues of a matrix
def eigenValues(mat):
a = -1*(mat[0][0]+mat[1][1]+mat[2][2])
b = -1*(mat[0][1]*mat[1][0]+mat[0][2]*mat[2][0]+mat[1][2]*mat[2][1]-mat[0][0]*mat[1][1]-mat[0][0]*mat[2][2]-mat[1][1]*mat[2][2])
c = -1*(mat.determinant())
eigenvalues = cubicCardano(a, b, c)
eigenvalues.sort()
return eigenvalues
# finds the roots of a cubic function
def cubicCardano(a,b,c):
p = b-((a**2)/3.0)
q = c+((2*a**3-9*a*b)/(27.0))
u = abs(-q/2.0 + abs((q**2)/4.0+(p**3)/27.0)**0.5)**(1/3.0)
roots = [u, u*complex(-0.5, (3**0.5)/2.0), u*complex(-0.5, -1*((3**0.5)/2.0))]
solutions = []
for j in roots:
x = j - p/(3.0*j) - a/3.0
try:
x = x.real
except:
pass
solutions.append(x)
return solutions
# make certain we're in objectmode, so we have up to date information
editmode = Window.EditMode()
Window.EditMode(0)
# getting the vertex locations
scn = bpy.data.scenes.active
me = scn.objects.active.getData(mesh = True)
locs = [v.co for v in me.verts]
# calculating the center of masss
com = Mathutils.Vector(0.0, 0.0, 0.0)
for loc in locs:
com += loc
com /= len(locs)
x, y, z = com
# creating the covariance matrix
mat = Mathutils.Matrix([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
for loc in locs:
mat[0][0] += (loc[0]-x)**2
mat[0][1] += (loc[0]-x)*(loc[1]-y)
mat[0][2] += (loc[0]-x)*(loc[2]-z)
mat[1][0] += (loc[1]-y)*(loc[0]-x)
mat[1][1] += (loc[1]-y)**2
mat[1][2] += (loc[1]-y)*(loc[2]-z)
mat[2][0] += (loc[2]-z)*(loc[0]-x)
mat[2][1] += (loc[2]-z)*(loc[1]-y)
mat[2][2] += (loc[2]-z)**2
# caculating the normal to the plane
eigenvalues = eigenValues(mat)
eigen = eigenvalues[2]
mat[0][0] -= eigen
mat[1][1] -= eigen
mat[2][2] -= eigen
normal = Mathutils.CrossVecs(mat[0], mat[1])/Mathutils.CrossVecs(mat[0], mat[1]).length
# project vertices onto the plane
for v in me.verts:
v.co = v.co - Mathutils.DotVecs((v.co-com), normal)*normal
# clean up
Window.Redraw()
Window.EditMode(editmode)
Any help is highly appreciated.
Thanks for your time.
