HW3-2 Linear Algebra - Naguib

1903 days ago by bigdata2016

In this problem, I have few 3D point cloud. the measurement of these points is erroneous. one of the interesting features we would like to measure for geometrical representation and matching is "surface normal".

 

simply, if these points are coming from a small patch in space, we may assume that this patch is almost planner. we are interested in modeling this plane, that is, finding the normalized vector perpendicular to this plane: "surface normal".

 

I generated 25 3D points on XY plane and added some noise to them.

 

I demonstrate 2 approaches to find the surface normal:

 

1. using SVD: surface normal is the Eigen vector with smallest Eigen value

2. using linear algebra: any 2 vectors formed by any 3 points will only exist in their plane (XY plane). thus, surface normal can be found by cross product. to minimize error, we can take average of multiple sampled estimation. there might be better ways than just average, but I think it's suitable enough for an example
import matplotlib.pyplot as plt from matplotlib import cm import numpy as np from mpl_toolkits.mplot3d import Axes3D from matplotlib.patches import FancyArrowPatch from mpl_toolkits.mplot3d import proj3d class Arrow3D(FancyArrowPatch): def __init__(self, xs, ys, zs, *args, **kwargs): FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs) self._verts3d = xs, ys, zs def draw(self, renderer): xs3d, ys3d, zs3d = self._verts3d xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M) self.set_positions((xs[0],ys[0]),(xs[1],ys[1])) FancyArrowPatch.draw(self, renderer) ################################################################################################################################ A=matrix(RDF, [[0,0 ,0],[5,0 ,0],[10,0 ,0],[-5,0 ,0],[-10,0 ,0], [0,5 ,0],[5,5 ,0],[10,5 ,0],[-5,5 ,0],[-10,5 ,0], [0,10 ,0],[5,10 ,0],[10,10 ,0],[-5,10 ,0],[-10,10 ,0], [0,-5 ,0],[5,-5 ,0],[10,-5 ,0],[-5,-5 ,0],[-10,-5 ,0], [0,-10,0],[5,-10,0],[10,-10,0],[-5,-10,0],[-10,-10,0]]) for i in range(25): A[i,0]=A[i,0]+ZZ.random_element(100)/24.0-2.0 A[i,1]=A[i,1]+ZZ.random_element(100)/24.0-2.0 A[i,2]=A[i,2]+ZZ.random_element(100)/10.0-5.0 pretty_print("25 noisy point vectors of a 3D plane:") #show(A) show(A) pretty_print("How to find the surface normal (red arrow). i.e., fit a 3D plane?") fig = plt.figure() ax = fig.add_subplot(111, projection='3d') x_surf=[-15, 15] y_surf=[-15, 15] x_surf, y_surf = np.meshgrid(x_surf, y_surf) z_surf=x_surf*0 ax.plot_surface(x_surf, y_surf, z_surf, cmap=cm.hot, alpha=0.2) ax.scatter(A.transpose()[0], A.transpose()[1], A.transpose()[2], c='r', marker='o') arrow = Arrow3D([0, 0], [0, 0], [0, 20], mutation_scale=20, lw=3, arrowstyle="-|>", color="r") ax.add_artist(arrow) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax.set_xlim(xmin=-15.0, xmax=15) ax.set_ylim(xmin=-15.0, xmax=15) ax.set_zlim(xmin=-5.0, xmax=25) #plt.show() plt.savefig('Scatter.png') plt.close() 
       
S,V,_D = A.SVD() D=copy(_D) if D[2,2]<0: D=D*-1 D[2]=D[2]/norm(D[2]) pretty_print("SVD approach: the vector with smallest eigen value:") show(D[2]) V=matrix(RDF,25^3,3) Vavg=matrix(RDF,1,3) Vavg[0,0]=0 Vavg[0,1]=0 Vavg[0,2]=0 pretty_print("Geometrical approach: sampling every 3 points, form 2 vectors") pretty_print("surface normal = cross product of 2 vectors") indx=0 for i in range(25): for j in range(25): if i!=j : for k in range(25): if k!=j and k!=i : V[indx]=(A[i]-A[k]).cross_product(A[j]-A[k]) if V[indx,2]<0 : V[indx]=V[indx]*-1 Vavg[0]=Vavg[0]+V[indx]/norm(V[indx]) indx=indx+1 Vavg=Vavg/norm(Vavg) pretty_print("normalized average vector:") show(Vavg) 
       
h=matrix(RDF,indx,3) for i in range(indx): h[i]=V[i]/norm(V[i]) pretty_print("Histogram of x component of generated normals") plt.hist(h.transpose()[0], bins=50, range=(-1.,1.)) plt.savefig('Histogram_x.png') plt.close() 
       
pretty_print("Histogram of y component of generated normals") plt.hist(h.transpose()[1], bins=50, range=(-1.,1.)) plt.savefig('Histogram_y.png') plt.close() 
       
pretty_print("Histogram of z component of generated normals") plt.hist(h.transpose()[2], bins=50, range=(-1.,1.)) plt.savefig('Histogram_z.png') plt.close()