I fitted a plane into a set of points in R3 using minimization of L2 error (taken from this post: https://gist.github.com/amroamroamro/1db8d69b4b65e8bc66a6).
data = cluster_data[:ALLOWED_PEAKS, :] grid_steps = 20 # regular grid covering the domain of the data mn = np.min(data, axis=0) mx = np.max(data, axis=0) X, Y = np.meshgrid(np.linspace(mn, mx, grid_steps), np.linspace(mn, mx, grid_steps)) XX = X.flatten() YY = Y.flatten() # best-fit linear plane (1st-order) A = np.c_[data[:, 0], data[:, 1], np.ones(data.shape)] C, _, _, _ = scipy.linalg.lstsq(A, data[:, 2]) # coefficients # evaluate it on grid Z = C * X + C * Y + C gci = int((grid_steps / 2) - 1 ) # grid_center_index centroid = np.array([X[gci, gci], Y[gci, gci], Z[gci, gci]]).reshape(1,-1) # take center point of grid steps surface_norm = np.array([-C, -C, 1]).reshape(1,-1) surface_norm_unit = normalize(surface_norm, norm='l2') centroid_norm = centroid + surface_norm_unit
The plane looks like a great fit by visual inspection.
Now I want to find the normal vector of the plane. When looking at the code, the plane was determined by using an equation of the form
my normal vector should be
[ C, C, -1] or
[ -C, -C, 1].
When plotting this normal vector, I receive a line that does not appear to be perpendicular on the plane.
Especially when looking from on side and making the plane appear to be a line (very thin, I hope you can see it)
Am I missing something fundamental about normal vectors? Is my assumption incorrect, that on the second image it should be orthogonal to the thin line? Here is my code for plotting:
# plot points and fitted surface using Matplotlib fig = plt.figure(figsize=(10, 10)) ax = fig.gca(projection='3d') ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2) ax.scatter(data[:, 0], data[:, 1], data[:, 2], c='r', s=50) ax.quiver(*[*centroid.reshape(-1,).tolist(),*surface_norm_unit.reshape(-1,).tolist()]) ax.scatter(centroid[0,0], centroid[0,1], centroid[0,2], c='b', s=100) ax.scatter(centroid_norm[0,0], centroid_norm[0,1], centroid_norm[0,2], c='b', s=100) plt.xlabel('X') plt.ylabel('Y') ax.set_zlabel('Z') ax.axis('equal') ax.axis('tight')
When multiplying the calculated normal vector with a vector on the plane, the result is
2.77555756e-17, so nearly 0, which should be an indicator of both vectors being orthogonal.