Minimize distances to plane

Question

I have a set of points in 3D space and I know that all of these points belong to a plane. However, there is some noise present in these points, so I cannot just extract a plane directly from it. I would like to find the formula of the plane(ax+by+c*z+d=0) that best fits these points. In other words, the sum of the (squared) distances from the points to the plane should be minimized.

I am doing all this using python and numpy, but I can't seem to figure out how to exactly implement this.

Does this answer your question? [Fit plane to a set of points in 3D: scipy.optimize.minimize vs scipy.linalg.lstsq](https://stackoverflow.com/questions/35070178/fit-plane-to-a-set-of-points-in-3d-scipy-optimize-minimize-vs-scipy-linalg-lsts) — Zaraki Kenpachi, Mar 11 '20 at 11:03
I did not immediatly figure out how to get it to work in my case(it worked for some planes, but it didn't for others), however, I searched some more knowing about those things and found something similar which did the trick for me. — DrRonne, Mar 11 '20 at 16:31
There is an interesting answer here: https://math.stackexchange.com/questions/99299/best-fitting-plane-given-a-set-of-points — GiulioSanto, Mar 11 '20 at 17:07

score 0 · Accepted Answer · answered Mar 11 '20 at 16:29

I did some more research on what python, scipy and numpy had to offer to solve this problem and I came across this piece of code which happened to do exactly what I needed:

def calcBestNormal(points, origin):
p = np.subtract(points, origin)

# Inital guess of the plane
p0 = np.array([0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129])

def f_min(X, p):
    plane_xyz = p[0:3]
    distance = (plane_xyz * X.T).sum(axis=1) + p[3]
    return distance / np.linalg.norm(plane_xyz)

def residuals(params, signal, X):
    return f_min(X, params)

sol = leastsq(residuals, p0, args=(None, p.T))[0]
normal = np.divide(sol[0:3], np.linalg.norm(sol[0:3]))
return normal

score 0 · Answer 2 · answered Mar 16 '20 at 00:19

If you calculate the singular value decomposition of the point distribution, the vector corresponding to the smallest eigenvalue will be normal to the plane adjusted to the distribution.

Assuming points is a point matrix containing a point per row:

centroid = np.mean(points, axis=0)               # point on the plane
normal = np.linalg.svd(points - centroid)[2][2]  # plane normal (a,b,c)
d = -centroid.dot(normal)                        # distance to the origin (d)
plane = np.append(normal, d)                     # plane coefficients (a,b,c,d)

Minimize distances to plane

2 Answers2