# Determination of plane parameters with PCA

Principal Component Analysis (PCA) calculates the eigenvector of point cloud as the normal vector for finding out the parameters of the optimal plane fitting the data. Program PCA as a function (input: point cloud as a $n \times 3$ array) and estimate corresponding planes of all three point clouds. This page is available for practice as an interactive jupyter notebook.

# imports
from plyfile import PlyData
import numpy as np

plane_tango = np.stack((_data['x'], _data['y'], _data['z']), axis=1)
plane_riegl1 = np.stack((_data['x'], _data['y'], _data['z']), axis=1)
plane_riegl2 = np.stack((_data['x'], _data['y'], _data['z']), axis=1)

print(plane_tango)
print(plane_riegl1)
print(plane_riegl2)

[[ 0.0693865   2.54306006  0.0111053 ]
[ 0.074412    2.54540992  0.0109674 ]
[ 0.0797195   2.54717994  0.010857  ]
...,
[ 0.321962    2.79351997 -0.0689245 ]
[ 0.330892    2.80245996 -0.0707128 ]
[ 0.339948    2.81185007 -0.0725622 ]]
[[ 0.36899999  2.77800012 -0.51099998]
[ 0.36899999  2.77600002 -0.63800001]
[ 0.37        2.78600001 -0.66500002]
...,
[-0.07        2.43600011 -0.248     ]
[-0.07        2.44099998 -0.37200001]
[-0.07        2.44099998 -0.54500002]]
[[ 1.39300001  2.71700001 -0.005     ]
[ 1.39300001  2.71600008 -0.017     ]
[ 1.38100004  2.70600009  0.034     ]
...,
[ 1.05400002  3.04299998 -0.505     ]
[ 1.05299997  3.03999996 -0.51099998]
[ 1.04299998  3.03299999 -0.38699999]]

def pca(points):
### BEGIN SOLUTION
#Reduce center of mass
com = np.sum(points, axis=0)/points.shape[0]

points = points - com

# determine plane

x = 0
y = 1
z = 2

A_xx = np.sum(points[:,x]*points[:,x], axis=0)
A_yy = np.sum(points[:,y]*points[:,y], axis=0)
A_zz = np.sum(points[:,z]*points[:,z], axis=0)
A_xy = np.sum(points[:,x]*points[:,y], axis=0)
A_xz = np.sum(points[:,x]*points[:,z], axis=0)
A_yz = np.sum(points[:,y]*points[:,z], axis=0)

ATA = [[A_xx, A_xy, A_xz],
[A_xy, A_yy, A_yz],
[A_xz, A_yz, A_zz]
]

w, v = np.linalg.eig(ATA)
lambda_ = w[w.argmin()]
n = v[:,w.argmin()]
n = n / np.linalg.norm(n)

d = -1* np.dot(n, com)
### END SOLUTION
return lambda_ , n, d  # plane paramters lambda, n: normal vector, d distance to


                  # Expected values
pca(plane_tango)  # 0.297639 [-0.69960505  0.71404618  0.02628372] -1.75849068165
pca(plane_riegl1) # 0.859337 [-0.67446429  0.73830175 -0.00290692] -1.8286588192
pca(plane_riegl2) # 0.747634 [-0.73383665 -0.67925513  0.00981276] 2.85578155518

[ 0.12740748  2.59833336 -0.29312617]
[ 0.1308195   2.59506202 -0.32786474]
[ 1.22546291  2.876863   -0.2413134 ]


## Calculate the noise 𝜎 for every point cloud

Calculate the noise using the iterative approach by calculating and summing single $\nu_i$ $\sigma=\sqrt{\frac{\sum \nu_i^2}{m-3}}$

def noise_it(points, n, d):
### BEGIN SOLUTION

### END SOLUTION
return sigma


## Compare the noise to the smallest eigenvalue 𝜆 from the PCA (cf. course script No.5, slide 23)

Calcultae the noise from smallest eigenvalue and write comparison to noise_it in the comment block.

def noise_lambda(lambda_, m):
### BEGIN SOLUTION

### END SOLUTION
return sigma