Orthogonal projection with numpy

Question

I have a list of three-dimensional points for which I calculate the plane using the numpy.linalg.lstsq method. But now I want to make an orthogonal projection for each point on this plane, but I can not find my error:

from numpy.linalg import lstsq def VecProduct(vek1, vek2): return (vek1[0]*vek2[0] + vek1[1]*vek2[1] + vek1[2]*vek2[2]) def CalcPlane(x, y, z): # x, y and z are given in lists n = len(x) sum_x = sum_y = sum_z = sum_xx = sum_yy = sum_xy = sum_xz = sum_yz = 0 for i in range(n): sum_x += x[i] sum_y += y[i] sum_z += z[i] sum_xx += x[i]*x[i] sum_yy += y[i]*y[i] sum_xy += x[i]*y[i] sum_xz += x[i]*z[i] sum_yz += y[i]*z[i] M = ([sum_xx, sum_xy, sum_x], [sum_xy, sum_yy, sum_y], [sum_x, sum_y, n]) b = (sum_xz, sum_yz, sum_z) a,b,c = lstsq(M, b)[0] ''' z = a*x + b*y + c a*x = z - b*y - c x = -(b/a)*y + (1/a)*z - c/a ''' r0 = [-c/a, 0, 0] u = [-b/a, 1, 0] v = [1/a, 0, 1] xn = [] yn = [] zn = [] # orthogonalize u and v with Gram-Schmidt to get u and w uu = VecProduct(u, u) vu = VecProduct(v, u) fak0 = vu/uu erg0 = [val*fak0 for val in u] w = [v[0]-erg0[0], v[1]-erg0[1], v[2]-erg0[2]] ww = VecProduct(w, w) # P_new = ((x*u)/(u*u))*u + ((x*w)/(w*w))*w for i in range(len(x)): xu = VecProduct([x[i], y[i], z[i]], u) xw = VecProduct([x[i], y[i], z[i]], w) fak1 = xu/uu fak2 = xw/ww erg1 = [val*fak1 for val in u] erg2 = [val*fak2 for val in w] erg = [erg1[0]+erg2[0], erg1[1]+erg2[1], erg1[2]+erg2[2]] erg[0] += r0[0] xn.append(erg[0]) yn.append(erg[1]) zn.append(erg[2]) return (xn,yn,zn) 

This returns me a list of points that are all in the plane, but when I show them, they are not in the positions that should be. I believe that there is already a certain built-in method to solve this problem, but I could not find any = (

Answer 1

You are using np.lstsq very poorly since you feed it with a precomputed 3x3 matrix instead of letting it do the job. I would do it like this:

 import numpy as np def calc_plane(x, y, z): a = np.column_stack((x, y, np.ones_like(x))) return np.linalg.lstsq(a, z)[0] >>> x = np.random.rand(1000) >>> y = np.random.rand(1000) >>> z = 4*x + 5*y + 7 + np.random.rand(1000)*.1 >>> calc_plane(x, y, z) array([ 3.99795126, 5.00233364, 7.05007326]) 

In fact, it is more convenient to use a formula for your plane, which does not depend on the coefficient, when z not equal to zero, i.e. uses a*x + b*y + c*z = 1 . You can also calculate a , b and c :

 def calc_plane_bis(x, y, z): a = np.column_stack((x, y, z)) return np.linalg.lstsq(a, np.ones_like(x))[0] >>> calc_plane_bis(x, y, z) array([-0.56732299, -0.70949543, 0.14185393]) 

To project points onto a plane using my alternative equation, the vector (a, b, c) perpendicular to the plane. It is easy to verify that the point (a, b, c) / (a**2+b**2+c**2) is on the plane, so the projection can be done by referencing all the points on this point on the plane, projecting points to a normal vector, subtracting the projection from the points, then referring to them back to the beginning. You can do it as follows:

 def project_points(x, y, z, a, b, c): """ Projects the points with coordinates x, y, z onto the plane defined by a*x + b*y + c*z = 1 """ vector_norm = a*a + b*b + c*c normal_vector = np.array([a, b, c]) / np.sqrt(vector_norm) point_in_plane = np.array([a, b, c]) / vector_norm points = np.column_stack((x, y, z)) points_from_point_in_plane = points - point_in_plane proj_onto_normal_vector = np.dot(points_from_point_in_plane, normal_vector) proj_onto_plane = (points_from_point_in_plane - proj_onto_normal_vector[:, None]*normal_vector) return point_in_plane + proj_onto_plane 

So now you can do something like:

 >>> project_points(x, y, z, *calc_plane_bis(x, y, z)) array([[ 0.13138012, 0.76009389, 11.37555123], [ 0.71096929, 0.68711773, 13.32843506], [ 0.14889398, 0.74404116, 11.36534936], ..., [ 0.85975642, 0.4827624 , 12.90197969], [ 0.48364383, 0.2963717 , 10.46636903], [ 0.81596472, 0.45273681, 12.57679188]]) 

Answer 2

You can just do everything in matrices - this is one option.

If you add your points to numerical vectors to the matrix X , and y is a vector, then the vector of beta parameters for solving the least squares:

 import numpy as np beta = np.linalg.inv(XTdot(X)).dot(XTdot(y)) 

but there is an easier way if we want to make projections: the QR decomposition gives us an orthonormalized projection matrix, since QT and Q itself is a matrix of orthonormalized basis vectors. So, we can first generate a QR , then get beta , then use QT to project the points.

QR:

 Q, R = np.linalg.qr(X) 

beta

 # use R to solve for beta # R is upper triangular, so can use triangular solver: beta = scipy.solve_triangular(R, QTdot(y)) 

So now we have beta , and we can project points using QT very simply:

 X_proj = QTdot(X) 

Here it is!

If you need more information and graphic drawings and stuff, I made a whole bunch of notes doing something like this: https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

(Edit: note that if you want to add an offset term, so the best match should not go through the origin, you can simply add an extra column with all-1s to X , which acts as a term / offset attribute)