Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def getEigenvaluePCA(a, b, c, dim, decim=2):
"""
Function to estimate the eigenvalues of the PCA
for every points into a given sphere of radius dim/2
------
INPUT:
@a, @b, @c: numpy arrays with the coordinates
@dim : float - diameter of interest of the neighborhood ball
@decim : integer - decimation factor to lower the time processing.
-----
OUTPUT:
@comp : array of array with the 3 eigenvalues
EXAMPLE:
[[0.68693784 0.31106007 0.0020021 ]
[0.75319626 0.24592732 0.00087642]
[0.78309201 0.21595601 0.00095197]
...]
"""
from sklearn.decomposition import PCA
from scipy import spatial
print('Data are decimed by a factor of ', decim)
# Remove the mean of each vector
a -= np.mean(a, axis=0)
b -= np.mean(b, axis=0)
c -= np.mean(c, axis=0)
# Concat the coordinates arrays into a single matrix.
Y = np.c_[a[::decim], b[::decim], c[::decim]]
print("shape of Y " ,np.shape(Y))
# Estimate the distance of one point with every neighbour around
dist = spatial.distance.squareform(spatial.distance.pdist(Y))
n_samples = Y.shape[0] # number of points in Y
print(' Treating distance of', dim, 'm')
comp = [] # Initialize the eigenvalues components array
# Loop over every points to have a look on its neighbour.
# 1. We consider a sphere of radius dim/2 and we only consider
# the neighbours inside this sphere
# 2. We create a subsample with these points
# 3. Using PCA scikit-learn function, we estimate the 3 eigenvalues of these points
# 4. We concatenate these 3 eigenvalues in a single array which will be return at the end.
for kk in range(0,n_samples):
pts = np.where((dist[kk,:]<=dim/2)) # Mask of values corresp. to the criteria
if np.size(pts)<3: # Check if the number of neighbour is sufficient
print("Revise dims: ")
print(pts)
print(dist)
else:
Ytmp = Y[pts,:] # Apply the criteria on the dataset
Ytmp = Ytmp[0,:,:] # Reduce the depth of the array
# (not necessary depending on the data structure used)
# Equivalent to Ytmp = np.squeeze(Ytmp)
pca = PCA(n_components=3) # Instantiate the PCA object
pca.fit(Ytmp)
eigenvalues = pca.explained_variance_ratio_
#sval = pca.singular_values_
#print('------')
#print( eigenvalues)
#print( sval**2 / np.sum(sval**2) )
comp.append(eigenvalues) # Append the 3 eigenvalues to "comp" array
comp = np.array(comp)
return comp