Fitting a plane by Orthogonal Regression in Python

Question

I want to fit a plane to a set of points (x, y, z) in Python. I found various answers how to perform the fitting if the error is measured with respect to the z-axis but I want to consider errors in orthogonal direction. I found the following question (Best fit plane by minimizing orthogonal distances) which addresses the same question - but it's not clear to me how to implement this in Python (likely with NumPy/SciPy). Further details regarding the mathematical derivation can also be found here: http://www.ncorr.com/download/publications/eberlyleastsquares.pdf (section 2).

Answer 1

The first link you gave does describe the algorithm for orthogonal distance fitting, but rather tersely. Here, in case it helps, is a more prolix description:

I suppose you have points (in your case 3d, but the dimension makes no odds to the algotithm) P[i], i=1..N You want to find a (hyper-) plane that is of mininmal orthogonal distance from your points.

A hyper-plane can be described by a unit vector n and a scalar d. The set of points on the plane is

{ P | n.P + d = 0 }

and the (orthogonal) distance of a point P from the plane is

n.P + d

So we want to find n and d to minimise

Q(n,d) = Sum{ i | (n.P[i]+d)*(n.P[i]+d) } /N

(The division by N isn't essential, and makes no difference to the values of n and d that are found, but to my mind makes the algebra neater)

The first thing to notice is that if we knew n, the d that minimises Q will be

d = -n.Pbar where
Pbar = Sum{ i | P[i]}/N, the mean of the P[]

We may as well use this value of d, so that, after a little algebra the problem reduces to minimising Q^:

Q^(n) = Sum{ i | (n.P[i]-n.Pbar)*(n.P[i]-n.Pbar) } /N
      = n' * C * n
where
C = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar) } /N

The form of Q^ tells us that the value of n to minimise Q^ will be an eigenvector of C correseponding to a minimal eigenvalue.

So (sorry I can't give code but my python is contemptible):

a/ compute

Pbar = Sum{ i | P[i]}/N, the mean of the points

b/ compute

C = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar) } /N, the covariance matrix of the points

c/ diagonalise C, and pick out a minimal eigenvalue and the corresponding eigenvector n

d/ compute

d = -Pbar.n

Then n, d define the hyperplane you want.

Answer 2

I've also had to deal with this situation and at first the mathematical notation can be overwhelming, but in the end the solution is fairly simple.

Once you get the intuition that the vector (A,B,C) that defines the best fitting plane Ax+By+Cz+D=0 is the one that explains the minimum variance of your set of coordinates, then the solution is straightforward.

First thing to do is center your coordinates (this way D will be 0 in your plane equation)

coords -= coords.mean(axis=0)

Then you have 2 options to get the vector you are interested in: (1) use the PCA implementation from sklearn or scipy to get the vector that explains minimal variance

pca = PCA(n_components=3)
pca.fit(coords)
# The last component/vector is the one with minimal variance, see PCA documentation
normal_vector = pca.components_[-1]

(2) re-implement the procedure described in the Geometric Tool reference you've linked.

@njit
def get_best_fitting_plane_vector(coords):

    # Calculate the covariance matrix of the coordinates
    covariance_matrix = np.cov(coords, rowvar=False) # Variables = columns

    # Calculate the eigenvalues & eigenvectors of the covariance matrix
    e_val, e_vect = np.linalg.eig(covariance_matrix)

    # The normal vector to the plane is the eigenvector associated to the minimum eigenvalue
    min_eval = np.argmin(e_val)
    normal_vector = e_vect[:, min_eval]

    return normal_vector

In terms of speed, the re-implemented procedure is faster than using PCA, and can be a lot faster if you use numba (just decorate the function with @njit).

Answer 3

Based on your second refernce []

Say you have n samples (x,y,z)

I'll call the 3 terms M*A=V, and define the column arrays

X=[ x_0, x_1 .. x_n ]'
Y=[ y_0, y_1 .. y_n ]'
Z=[ z_0, z_1 .. z_n ]'

Define the (n by 3) matrix XY1=[X,Y,1n]:

      [[x_0,y_0,1],
XY1=   [x_1,y_1,1],
       ...
       [x_n,y_n,1]]

The matrix M can be obtained as

M = XY1' * XY1

Where apostrophe (') is the transposition operator and (*) the matrix product.

And the array V is

V = XY1'*Z

The least squares solution can be obtained through the moore-penrose pseoudoinverse: [(M'*M)^-1 * M']

~A  = [(M'*M)^-1 * M'] * V

Sample code:

import numpy as np
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt

#Input your values
A=3
B=2
C=1

#reserve memory
xy1=np.ones([n,3])

#Make random data, n  ( x,y ) tuples.
n=30 #samples
xy1[:,:2]=np.random.rand(n,2)

#plane: A*x+B*y+C = z , the z coord is calculated from random x,y
z=xy1.dot (np.array([[A,B,C],]).transpose() )
    
#addnoise
xy1[:,:2]+=np.random.normal(scale=0.05,size=[n,2])
z+=np.random.normal(scale=0.05,size=[n,1])

#calculate M and V    
M=xy1.transpose().dot(xy1)
V=xy1.transpose().dot(z)

#pseudoinverse:
Mp=np.linalg.inv(M.transpose().dot(M)).dot(M.transpose())

#Least-squares Solution
ABC= Mp.dot(V) 

Output

In [24]: ABC
Out[24]: 
array([[3.11395111],
       [2.02909874],
       [1.01340411]])