https://math.stackexchange.com/a/99317 and https://stackoverflow.com/q/15959411 give two different ways to use SVD to fit a plane to a set of points. The former first subtracts out the centroid of the points, constraining the plane to contain the centroid. The latter does not do this, which gives an extra degree of freedom to the SVD. The two solutions sometimes give slightly different solutions, as the following code illustrates:
import numpy as np
def fitPlaneSVDNoConstraint(XYZ):
[rows,cols] = XYZ.shape
# 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.
p = (np.ones((rows,1)))
AB = np.hstack([XYZ,p])
[u, d, v] = np.linalg.svd(AB,0)
B = v[3,:]; # Solution is last column of v.
nn = np.linalg.norm(B[0:3])
B = B / nn
return B[0:3]
def fitPlaneSVDNoConstraintConstrainToCentroid(XYZ):
[rows,cols] = XYZ.shape
# 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.
AB = XYZ - np.mean(XYZ, axis=0)
[u, d, v] = np.linalg.svd(AB,0)
B = v[2,:]; # Solution is last column of v.
nn = np.linalg.norm(B[0:3])
B = B / nn
return B[0:3]
if name=="main":
XYZ = np.array([
[0,0,1],
[0,1,2],
[0,2,3],
[1,0,1],
[1,1,2],
[1,2,3],
[2,0,1],
[2,1,2],
[2,2,3]
])
print( "SVD, constrained to centroid: ",fitPlaneSVDNoConstraint(XYZ))
print( "SVD, no constraint : ",fitPlaneSVDNoConstraintConstrainToCentroid(XYZ))
XYZ = np.array([
[ -1.38446357, -1.88062266, -1.54779387],
[ 0.17198868, -1.80153589, -0.21210896],
[ 0.44548632, 1.42599871, 1.79542404],
[ 0.08578674, 0.14314378, 0.05247274],
[ 9.03135451, -0.90611959, -11.72476911]
])
print( "SVD, constrained to centroid: ",fitPlaneSVDNoConstraint(XYZ))
print( "SVD, no constraint : ",fitPlaneSVDNoConstraintConstrainToCentroid(XYZ))
This gives:
SVD, constrained to centroid: [ 0. 0.70710678 -0.70710678]
SVD, no constraint : [-0. 0.70710678 -0.70710678]
SVD, constrained to centroid: [ 0.57230214 -0.67672136 0.46316138]
SVD, no constraint : [-0.57665504 0.66508349 -0.47448174]
The first case gives the same result, but the second case gives two slightly difference results. The sign difference doesn't matter since the direction of the normal vector to the plane is arbitrary, but I'm having trouble figuring out why exactly there's a difference. Why is there a difference, and which method would be preferable?
