Note to self: Find out about residuals of planefits using singular value decomposition. And while you're at it: find out, how least square fitting works with singular value decomposition anyhow. Right now - as with any transformation matrix bigger than 3x3 - i mostly think it's magic..
Here is my current planefitting code in python:
-
import numpy as np
-
from Geometry import Point, Line, Plane
-
-
def fitplane(XYZ):
-
[rows,npts] = XYZ.shape
-
-
if not rows == 3:
-
print XYZ.shape
-
raise ('data is not 3D')
-
return None
-
-
if npts <3:
-
raise ('too few points to fit plane')
-
return None
-
-
# Set up constraint equations of the form AB = 0,
-
# where B is a column vector of the plane coefficients
-
# in the form b(1)*X + b(2)*Y +b(3)*Z + b(4) = 0.
-
t = XYZ.T
-
p = (np.ones((npts,1)))
-
A = np.hstack([t,p])
-
-
if npts == 3: # Pad A with zeros
-
A = [A, np.zeros(1,4)]
-
-
[u, d, v] = np.linalg.svd(A) # Singular value decomposition.
-
#print v[3,:]
-
B = v[3,:]; # Solution is last column of v.
-
nn = np.linalg.norm(B[0:3])
-
B = B / nn
-
plane = Plane(Point(B[0],B[1],B[2]),D=B[3])
-
return plane
The geometry-module is the one i've posted on github. It enables me to do fast transformations and - the most important part - intersections between planes and lines. Go check that out. I would like to hear some feedback..
Back to the problem: I would like to also return a measure of error or residual. If anyone could help..
Update (25.10.2010): There was a small but still significant error in the algorithm. I assumed that B[0:2] returns the elements [B[0],B[1],B[2]], but it appears the slicing in python is in line with the syntax for a "for"-statement. So, the line with the norm had to be changed to
-
nn = np.linalg.norm(B[0:3])
On another note, i'd have to look up the np.linalg.solve()-function which should be able to solve my problem a little bit more efficiently.