(Although there are a number of questions regarding how to best fit a plane to some 3D data on SO, I couldn't find an answer for this issue.)
Given N (x, y, z) points, I need the best fit plane
a*x + b*y + c*z + d = 0
defined through the a, b, c, d coefficients that minimize the mean of the orthogonal distances from the points to the plane. The point-plane orthogonal distance (for a given (x0, y0, z0) point) is defined as:
d = |a*x0 + b*y0 + c*z0 + d|/sqrt(a^2 + b^2 + c^2)
I set up two methods (code below):
- Singular-Value Decomposition (source)
- Basin-Hopping minimization of the mean orthogonal distances
As I understand it, the SVD method should produce the exact best fit plane by minimizing the orthogonal distances analytically. What I find instead is that the BH method gives better results than the supposedly exact SVD method, even for a low number of BH runs.
By "better" I mean that the final mean orthogonal distance value is smaller with the BH method, than with the SVD method.
What am I missing here?
import numpy as np
import scipy.optimize as optimize
def perp_error(params, xyz):
"""
Mean of the absolute values for the perpendicular distance of the
'xyz' points, to the plane defined by the coefficients 'a,b,c,d' in
'params'.
"""
a, b, c, d = params
x, y, z = xyz
length = np.sqrt(a**2 + b**2 + c**2)
return (np.abs(a * x + b * y + c * z + d) / length).mean()
def minPerpDist(x, y, z, N_min):
"""
Basin-Hopping method, minimize mean absolute values of the
orthogonal distances.
"""
def unit_length(params):
"""
Constrain the vector perpendicular to the plane to be of unit length.
"""
a, b, c, d = params
return a**2 + b**2 + c**2 - 1
# Random initial guess for the a,b,c,d plane coefficients.
initial_guess = np.random.uniform(-10., 10., 4)
# Constrain the vector perpendicular to the plane to be of unit length.
cons = ({'type': 'eq', 'fun': unit_length})
min_kwargs = {"constraints": cons, "args": [x, y, z]}
# Use Basin-Hopping to obtain the best fit coefficients.
sol = optimize.basinhopping(
perp_error, initial_guess, minimizer_kwargs=min_kwargs, niter=N_min)
abcd = list(sol.x)
return abcd
def SVD(X):
"""
Singular value decomposition method.
Source: https://gist.github.com/lambdalisue/7201028
"""
# Find the average of points (centroid) along the columns
C = np.average(X, axis=0)
# Create CX vector (centroid to point) matrix
CX = X - C
# Singular value decomposition
U, S, V = np.linalg.svd(CX)
# The last row of V matrix indicate the eigenvectors of
# smallest eigenvalues (singular values).
N = V[-1]
# Extract a, b, c, d coefficients.
x0, y0, z0 = C
a, b, c = N
d = -(a * x0 + b * y0 + c * z0)
return a, b, c, d
# Generate a random plane.
seed = np.random.randint(100000)
print("Seed: {}".format(seed))
np.random.seed(seed)
a, b, c, d = np.random.uniform(-10., 10., 4)
print("Orig abc(d=1): {:.3f} {:.3f} {:.3f}\n".format(a / d, b / d, c / d))
# Generate random (x, y, z) points.
N = 200
x, y = np.random.uniform(-5., 5., (2, N))
z = -(a * x + b * y + d) / c
# Add scatter in z.
z = z + np.random.uniform(-.2, .2, N)
# Solve using SVD method.
a, b, c, d = SVD(np.array([x, y, z]).T)
print("SVD abc(d=1): {:.3f} {:.3f} {:.3f}".format(a / d, b / d, c / d))
# Orthogonal mean distance
print("Perp err: {:.5f}\n".format(perp_error((a, b, c, d), (x, y, z))))
# Solve using Basin-Hopping.
abcd = minPerpDist(x, y, z, 500)
a, b, c, d = abcd
print("BH abc(d=1): {:.3f} {:.3f} {:.3f}".format(a / d, b / d, c / d))
print("Perp err: {:.5f}".format(perp_error(abcd, (x, y, z))))