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

# load data
_data = PlyData.read('plane_tango.ply').elements[0].data
plane_tango = np.stack((_data['x'], _data['y'], _data['z']), axis=1)
_data = PlyData.read('plane_riegl1.ply').elements[0].data
plane_riegl1 = np.stack((_data['x'], _data['y'], _data['z']), axis=1) 
_data = PlyData.read('plane_riegl2.ply').elements[0].data
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$ 

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

Compare the noise of the different planes to each other? How can the values be interpreted? Where do possible differences arise from?

Calculate the normal vectors and their intersection angles. Are they perpendicular to each other?

Author: Artem Leichter
Last modified: 2019-04-30