# HW3-2 Linear Algebra-락가바

## 1917 days ago by bigdata2016

Since i'm not sure what kind of problems we should try to solve, I will try to use problems relevant to my research field as sample problems. In  , and . Since  we have .      □

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

solution is attached and sage code can be found here:

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,ys),(xs,ys)) 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(), A.transpose(), A.transpose(), 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()
 $\newcommand{\Bold}{\mathbf{#1}}\verb|25|\phantom{\verb!x!}\verb|noisy|\phantom{\verb!x!}\verb|point|\phantom{\verb!x!}\verb|vectors|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|a|\phantom{\verb!x!}\verb|3D|\phantom{\verb!x!}\verb|plane:|$ \newcommand{\Bold}{\mathbf{#1}}\left(\begin{array}{rrr} -1.875 & -1.375 & 1.9 \\ 6.83333333333 & 0.125 & -4.7 \\ 8.0 & 1.0 & 0.2 \\ -5.83333333333 & -1.16666666667 & -3.4 \\ -8.70833333333 & 1.54166666667 & 2.6 \\ 0.333333333333 & 4.04166666667 & -1.5 \\ 4.375 & 5.95833333333 & 2.4 \\ 9.125 & 3.91666666667 & 1.6 \\ -5.91666666667 & 3.45833333333 & -4.2 \\ -7.95833333333 & 3.95833333333 & 2.9 \\ 1.79166666667 & 8.91666666667 & 0.5 \\ 6.08333333333 & 9.70833333333 & -2.1 \\ 12.0 & 8.125 & -2.7 \\ -3.875 & 9.625 & -2.7 \\ -11.1666666667 & 11.6666666667 & 1.6 \\ -1.95833333333 & -5.0 & -5.0 \\ 4.91666666667 & -6.375 & 0.1 \\ 10.0416666667 & -5.16666666667 & -4.7 \\ -6.29166666667 & -4.66666666667 & 1.2 \\ -11.625 & -5.45833333333 & -2.9 \\ 1.16666666667 & -10.75 & -4.1 \\ 3.75 & -11.4583333333 & 4.1 \\ 11.0 & -9.16666666667 & 4.5 \\ -7.0 & -10.7916666667 & -0.6 \\ -8.125 & -8.25 & -4.7 \end{array}\right) $\newcommand{\Bold}{\mathbf{#1}}\verb|How|\phantom{\verb!x!}\verb|to|\phantom{\verb!x!}\verb|find|\phantom{\verb!x!}\verb|the|\phantom{\verb!x!}\verb|surface|\phantom{\verb!x!}\verb|normal|\phantom{\verb!x!}\verb|(red|\phantom{\verb!x!}\verb|arrow).|\phantom{\verb!x!}\verb|i.e.,|\phantom{\verb!x!}\verb|fit|\phantom{\verb!x!}\verb|a|\phantom{\verb!x!}\verb|3D|\phantom{\verb!x!}\verb|plane?|$ S,V,_D = A.SVD() D=copy(_D) if D[2,2]<0: D=D*-1 D=D/norm(D) pretty_print("SVD approach: the vector with smallest eigen value:") show(D) 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=Vavg+V[indx]/norm(V[indx]) indx=indx+1 Vavg=Vavg/norm(Vavg) pretty_print("normalized average vector:") show(Vavg)
 $\newcommand{\Bold}{\mathbf{#1}}\verb|SVD|\phantom{\verb!x!}\verb|approach:|\phantom{\verb!x!}\verb|the|\phantom{\verb!x!}\verb|vector|\phantom{\verb!x!}\verb|with|\phantom{\verb!x!}\verb|smallest|\phantom{\verb!x!}\verb|eigen|\phantom{\verb!x!}\verb|value:|$ \newcommand{\Bold}{\mathbf{#1}}\left(-0.0350423967582,\,0.0207107565043,\,0.999171204046\right) $\newcommand{\Bold}{\mathbf{#1}}\verb|Geometrical|\phantom{\verb!x!}\verb|approach:|\phantom{\verb!x!}\verb|sampling|\phantom{\verb!x!}\verb|every|\phantom{\verb!x!}\verb|3|\phantom{\verb!x!}\verb|points,|\phantom{\verb!x!}\verb|form|\phantom{\verb!x!}\verb|2|\phantom{\verb!x!}\verb|vectors|$ $\newcommand{\Bold}{\mathbf{#1}}\verb|surface|\phantom{\verb!x!}\verb|normal|\phantom{\verb!x!}\verb|=|\phantom{\verb!x!}\verb|cross|\phantom{\verb!x!}\verb|product|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|2|\phantom{\verb!x!}\verb|vectors|$ __main__:38: DeprecationWarning: The default norm will be changing from p='frob' to p=2. Use p='frob' explicitly to continue calculating the Frobenius norm. See http://trac.sagemath.org/13643 for details. $\newcommand{\Bold}{\mathbf{#1}}\verb|normalized|\phantom{\verb!x!}\verb|average|\phantom{\verb!x!}\verb|vector:|$ \newcommand{\Bold}{\mathbf{#1}}\left(\begin{array}{rrr} -0.0161161065356 & -0.0344666286846 & 0.999275899148 \end{array}\right) __main__:38: DeprecationWarning: The default norm will be changing from p='frob' to p=2. Use p='frob' explicitly to continue calculating the Frobenius norm. See http://trac.sagemath.org/13643 for details.
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(), bins=50, range=(-1.,1.)) plt.savefig('Histogram_x.png') plt.close()
 $\newcommand{\Bold}{\mathbf{#1}}\verb|Histogram|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|x|\phantom{\verb!x!}\verb|component|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|generated|\phantom{\verb!x!}\verb|normals|$ pretty_print("Histogram of y component of generated normals") plt.hist(h.transpose(), bins=50, range=(-1.,1.)) plt.savefig('Histogram_y.png') plt.close()
 $\newcommand{\Bold}{\mathbf{#1}}\verb|Histogram|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|y|\phantom{\verb!x!}\verb|component|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|generated|\phantom{\verb!x!}\verb|normals|$ pretty_print("Histogram of z component of generated normals") plt.hist(h.transpose(), bins=50, range=(-1.,1.)) plt.savefig('Histogram_z.png') plt.close()
 $\newcommand{\Bold}{\mathbf{#1}}\verb|Histogram|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|z|\phantom{\verb!x!}\verb|component|\phantom{\verb!x!}\verb|of|\phantom{\verb!x!}\verb|generated|\phantom{\verb!x!}\verb|normals|$ 