How to function fit a plane through a collection with points with minimal square root of error

Question

I'm working on a computer program that has to stabilize a set of points which should all appear on a plane in 3D space. Currently the program does not use the knowledge that all points should appear on a plane. So I want to write a program that projects these points on a plane. Now I know the math on how to project points onto a plane along a vector. Additionally rather then fitting for minimal distance I have noticed that the program tends to make errors in estimating distance but this error typically only occurs for a few data points meaning that I think it's best to function fit to the square root of the distance rather then the distance or the square of the distance (as I know is more common).

Note that while I'm programming in unity/C# I'm fine with a more or less mathematical answer (as long as the notation is not overly exotic)

Answer 1

I'd just write the error (in your case the sum of the square root of the distance of each point to the plane) as a function of the parameters that define that plane. After that I'd try to minimize the value of the function using any technique you like, e.g. the Nelder–Mead method or some kind of Gradient Descent with step size control or a Newton's Method, as you can easily calculate the derivatives of your function.

But I fear that the square root of the distances will not result in any useful fits. If you first want to play around with different kinds of error functions for fitting your plane (or just a line) I suggest trying it in Octave or Matlab where it is very easy to do. You can e.g. Just program a function that returns the error and use e.g. fminsearch (which uses the Nelder-Mead method)

Answer 2

If your points really are nearly on a plane $d=ax+by+cz$, then you get a fast result using the SVD. Construct a matrix $A$ with rows $$ w_k·[1,x_k,y_k,z_k] $$ where the weights $w_k$ will determine the importance of the error of that point in the overall equation, everything from $w_k=1$ to $w_k=(1+x_k^2-y_k^2-z_k^2)^{-\frac12}$ to $w_k=(1+x_k^2-y_k^2-z_k^2)^{-1}$ might be sensible (provided that $0$ is close to the center of the cluster of points).

Then $$ A=U\Sigma V^T $$ should have $2$ large singular values, one of large to middle size corresponding to the subspace spanned by the vectors $(a,d,0,0)$, $(b,0,d,0)$ and $(c,0,0,d)$. The remaining smallest or zero singular value has $v_4$ proportional (up to small errors) to $(-d,a,b,c)$. Where you can now read off the coefficients of the plane equation.