Skip to content
Snippets Groups Projects
Commit 1d6b7bd0 authored by Antoine Lucas's avatar Antoine Lucas ⛷️
Browse files

Replace libLab3_Lidar.py

parent d1452021
Branches main
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Labs/Lab_3_LiDAR_ML_Segmentation/libLab3_Lidar.py
+ 237
− 183
View file @ 1d6b7bd0
......@@ -4,13 +4,15 @@
Created on Mon Jan 18 15:09:57 2021
@authors: Grégory Sainton & Antoine Lucas
@purpose: Lib linked to the Lab3 notebook EDS_3D_LiDAR_Claiff_ML_*.ipynb
@purpose: Lib linked to the Lab3 notebook EDS_2020_2021_3D_LiDAR_Claiff_ML_*.ipynb
It contains useful function to run the lab, to plot several sets of data
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
......@@ -25,20 +27,20 @@ def readData(filename):
-------
x, y, z: numpy arrays
"""
"""
import numpy as np
import os
if os.path.isfile(filename):
f = open(filename)
floor = np.genfromtxt(filename)
x = floor[:,0]
x = floor[:,0]
y = floor[:,1]
z = floor[:,2]
f.close()
return x, y, z
else:
print(f'File {filename} is not accessible')
......@@ -51,28 +53,28 @@ def plot_figs(fig_num, elev, azim, dx, dy, dz, density=2):
INPUT:
------
@fig_num: int - Figure number
@elev : float - Elevation
@elev : float - Elevation
@azim : float - Azimuth
@dx, @dy, dz : numpy arrays - data to plot
@density: int - optionnal - factor to decimate the data on plot
OUTPUT:
-------
None
None
"""
# Subsample input data with a factor "density"
X = dx[::density]
Y = dy[::density]
Z = dz[::density]
# Plot the data in 3D
fig = plt.figure(figsize=(10, 8))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)
ax.scatter(dx[::density], dy[::density], dz[::density], c=dz[::density], marker='+', alpha=.4)
# OPTIONNAL - BUT NICE TO IMPROVE YOUR POINT OF VIEW
# Create a blanck cubic bounding box to simulate equal aspect ratio in 3D plots
max_range = np.array([X.max()-X.min(), Y.max()-Y.min(), Z.max()-Z.min()]).max()
......@@ -86,125 +88,150 @@ def plot_figs(fig_num, elev, azim, dx, dy, dz, density=2):
def getEigenvaluePCA(a, b, c, dim, decim=1,Verbose=False):
"""
Function to estimate the eigenvalues of the PCA
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.
"""
from sklearn.decomposition import PCA
from sklearn import preprocessing
from scipy import spatial
comp = None
Y = None
return comp, Y
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
def estimateTernaryCoord(comp):
"""
Function which estimate the ternary coordinates
"""
Function which estimate the ternary coordinates
from the list of eigenvalues.
----
INPUT:
@comp: np array - list of eigenvalue for every
@comp: np array - list of eigenvalue for every
points with a given dimension
OUTPUT:
-------
@X, @Y, @pdf
"""
# Conversion towards a,b,c space
# Conversion towards a,b,c space
#
from scipy.stats import gaussian_kde
P1 = comp[:,0]
P2 = comp[:,1]
P3 = comp[:,2]
# 2021-01/25 - Bug fixed in the calculation of xp1p2
# Consider a point of coordinates (p1, p2) -> Red point on the graphic above.
# What is the projection of this point on the line (p1p2) ?
# Estimate distances along p1-p2 axis
# Estimate a,b and c coordinates from the results just above
# Norm factor v: the sum of the eigenvalue (p1, p2 and p3 is equal to 1)
# Conversion towards ternary graph
# Consider a point of coordinates (p1, p2) -> Red point on the graphic above.
# What is the projection of this point on the line (p1p2) ?
# What is the projection of this point on the line (p1p2) ?
# We know that p1+p2+p3 = 1 -> p3 = 1-p1-p2
# To find the projection x_p1p2 : x_p1p2 = p1 + dx
# To find the projection x_p1p2 : x_p1p2 = p1 + dx
# with dx = p3.cos(theta)
# with theta = angle of the slope (-pi/4)
# It comes x_p1p1 = p1 + (1-p1-p2) * cos(theta)
xp1p2 = P1 + (1-P2-P1)*np.cos(np.pi/4)
# yp1p2 is derive from the fact that the slope is -1.
yp1p2 = 1 - xp1p2
# Estimate distances along p1-p2 axis
# Estimate a,b and c coordinates from the results just above
SmallDelta = np.sqrt((.5 - xp1p2)**2 + (.5 - yp1p2)**2)
BigDelta = np.sqrt((0.5 - 1)**2 + (0.5 - 0)**2)
# Estimate a,b and c coordinates from the results just above
c = 3 * P3
a = (1 - c)*SmallDelta/BigDelta
b = 1 - (c + a)
# Norm factor: the sum of the eigenvalue (p1, p2 and p3 is equal to 1)
# Conversion towards ternary graph
X = None
Y = None
v = (a + b + c)
# Finally, let's compute density probability function for clarity
# Conversion towards ternary graph
X = .5 * ( 2.*b+c ) / v
Y = 0.5*np.sqrt(3) * c / v
# Compute density probability function for clarity
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html
from scipy.stats import gaussian_kde
xy = np.vstack([X, Y])
pdf = gaussian_kde(xy)(xy)
return X, Y, pdf
xy = np.vstack([X, Y])
pdf = gaussian_kde(xy)(xy)
return X, Y, pdf
def density_scatter( x , y, ax = None, sort = True, bins = 20, **kwargs ):
"""
Function to plot a density plot from x,y data
----
"""
from matplotlib.colors import Normalize
from scipy.interpolate import interpn
import matplotlib.tri as tri
if ax is None :
fig , ax = plt.subplots()
data , x_e, y_e = np.histogram2d( x, y, bins = bins, density = True )
z = interpn( ( 0.5*(x_e[1:] + x_e[:-1]) , 0.5*(y_e[1:]+y_e[:-1]) ) , data , np.vstack([x,y]).T , method = "splinef2d", bounds_error = False)
#To be sure to plot all data
z[np.where(np.isnan(z))] = 0.0
# Sort the points by density, so that the densest points are plotted last
if sort :
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
ax.scatter( x, y, c=z, s= 0.5 )
norm = Normalize(vmin = np.min(z), vmax = np.max(z))
# creating the TRI grid
corners = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)*0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
refiner = tri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=4)
# plotting the mesh
plt.triplot(trimesh,'k-', linewidth=0.1)
plt.axis('equal')
plt.axis('off')
plt.show()
def plotTernary(x, y, z, title):
def plotTernary(x, y, z, title):
"""
Function to plot a ternary graph from list of ternary coordinates.
----
......@@ -214,15 +241,15 @@ def plotTernary(x, y, z, title):
----
OUTPUT:
None
"""
import matplotlib.tri as tri
# Plot part
plt.figure(figsize=(7, 5))
plt.clf()
plt.scatter(x,y,c=z, s= 0.5)
plt.scatter(x,y,c=z, s= 0.5)
# creating the TRI grid
corners = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)*0.5]])
......@@ -233,42 +260,27 @@ def plotTernary(x, y, z, title):
# plotting the mesh
plt.triplot(trimesh,'k-', linewidth=0.1)
plt.title(title)
plt.axis('equal')
plt.axis('off')
plt.show()
def plot_3dcladd(dx,dy,dz,y,density,fileout):
X = dx[::density]
Y = dy[::density]
Z = dz[::density]
elev = 36
azim = -144
# Plot the data in 3D
fig = plt.figure(figsize=(10, 8))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)
ax.scatter(dx[::density], dy[::density], dz[::density], c=y, marker='.', alpha=.5,cmap='YlGn')
# OPTIONNAL - BUT NICE TO IMPROVE YOUR POINT OF VIEW
# Create a blanck cubic bounding box to simulate equal aspect ratio in 3D plots
max_range = np.array([X.max()-X.min(), Y.max()-Y.min(), Z.max()-Z.min()]).max()
Xb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][0].flatten() + 0.5*(X.max()+X.min())
Yb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][1].flatten() + 0.5*(Y.max()+Y.min())
Zb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][2].flatten() + 0.5*(Z.max()+Z.min())
for xb, yb, zb in zip(Xb, Yb, Zb):
ax.plot([xb], [yb], [zb], 'w')
plt.savefig(fileout)
def plot_contours_compare(X, y, X_train , y_train, X_test, y_test, classifiers, plot_input = False):
"""
Function to plot in the same plot several classifiers
Function to plot in the same plot several classifiers
contour to compare them.
The librairies of the classifiers are supposed to be loaded out of the functions
since we don't know a priori the type of classifiers in input.
......@@ -283,30 +295,30 @@ def plot_contours_compare(X, y, X_train , y_train, X_test, y_test, classifiers,
@classifiers: list - list of instance of classifiers
ie : classifiers = [LinearDiscriminantAnalysis(store_covariance=True),
DecisionTreeClassifier()]
plot_input: boolean - optionnal - Option to plot the input
plot_input: boolean - optionnal - Option to plot the input
data alone in a separate plot
----
OUTPUT
None
EXAMPLE:
EXAMPLE:
-------
classifiers = [
LinearDiscriminantAnalysis(store_covariance=True),
SVC(kernel="linear", C=0.025),
QuadraticDiscriminantAnalysis(),
DecisionTreeClassifier(),
RandomForestClassifier()]
QuadraticDiscriminantAnalysis(),
DecisionTreeClassifier(),
RandomForestClassifier()]
plot_contours_compare(X, y, X_train , y_train, X_test, y_test, classifiers, plot_input = False)
"""
# Code adapted from
# Code adapted from
# https://scikit-learn.org/stable/auto_examples/classification/plot_classifier_comparison.html
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib import colors
cmap = colors.LinearSegmentedColormap(
......@@ -316,24 +328,24 @@ def plot_contours_compare(X, y, X_train , y_train, X_test, y_test, classifiers,
'blue': [(0, 0.7, 0.7), (1, 1, 1)]})
plt.cm.register_cmap(cmap=cmap)
nb_plot=len(classifiers)+1 if plot_input else len(classifiers)
# Here is defined the list of classifier's name
names = ["LDA","Linear SVM", "QDA", "Decision Tree", "Random Forest"]
figure = plt.figure(figsize=(27, 9))
i = 1
h = .02 # step size in the mesh
x_min, x_max = X.iloc[:, 0].min() - .5, X.iloc[:, 0].max() + .5
y_min, y_max = X.iloc[:, 1].min() - .5, X.iloc[:, 1].max() + .5
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
ds_cnt = 0
# just plot the dataset first
cm = plt.cm.RdBu
#cm_bright = ListedColormap(['#FF0000', '#0000FF'])
......@@ -370,11 +382,11 @@ def plot_contours_compare(X, y, X_train , y_train, X_test, y_test, classifiers,
# Put the result into a color plot
Z = Z.reshape(xx.shape)
#ax.contourf(xx, yy, Z, cmap=cm, alpha=.8)
ax.pcolormesh(xx, yy, Z, cmap='red_blue_classes',
norm=colors.Normalize(0., 1.), zorder=0,shading='auto')
ax.contour(xx, yy, Z, [0.5], linewidths=2., colors='white')
#ax.contourf(xx, yy, Z, cmap=cm, alpha=.8)
#'red_blue_classes'
......@@ -398,64 +410,106 @@ def plot_contours_compare(X, y, X_train , y_train, X_test, y_test, classifiers,
plt.show()
###
def plot_contours(X, y, classifier, resolution=0.02):
"""
Function to plot single classifier contours
----
INPUT:
@X: Numpy array with data
@y: Numpy array with labels
----
OUTPUT:
None
----
EXAMPLE:
plot_contours(X_train.to_numpy(), y_train.to_numpy(), lda)
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib import colors
# setup marker generator and color map
markers = ('s', 'x', 'o', '^', 'v')
colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
cmap = ListedColormap(colors[:len(np.unique(y))])
C = X[:, 0:2]
# plot the decision surface
x1_min, x1_max = C[:, 0].min() - 1, C[:, 0].max() + 1
x2_min, x2_max = C[:, 1].min() - 1, C[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
np.arange(x2_min, x2_max, resolution))
Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
Z = Z.reshape(xx1.shape)
plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
plt.xlim(xx1.min(), xx1.max())
plt.ylim(xx2.min(), xx2.max())
for idx, cl in enumerate(np.unique(y)):
plt.scatter(x=C[y == cl, 0], y=C[y == cl, 1],
alpha=0.8, c=cmap(idx),
marker=markers[idx], label=cl)
def plot_3dcladd(dx,dy,dz,y,density=1):
X = dx[::density]
Y = dy[::density]
Z = dz[::density]
elev = 36
azim = -144
# Plot the data in 3D
fig = plt.figure(figsize=(10, 8))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)
ax.scatter(dx[::density], dy[::density], dz[::density], c=y, marker='.', alpha=.5,cmap='YlGn')
# OPTIONNAL - BUT NICE TO IMPROVE YOUR POINT OF VIEW
# Create a blanck cubic bounding box to simulate equal aspect ratio in 3D plots
max_range = np.array([X.max()-X.min(), Y.max()-Y.min(), Z.max()-Z.min()]).max()
Xb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][0].flatten() + 0.5*(X.max()+X.min())
Yb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][1].flatten() + 0.5*(Y.max()+Y.min())
Zb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][2].flatten() + 0.5*(Z.max()+Z.min())
for xb, yb, zb in zip(Xb, Yb, Zb):
ax.plot([xb], [yb], [zb], 'w')
plt.show()
def plot_confMat2(cnf_matrix,ClaName,ClassId):
def plot_confMat(cnf_matrix,ClaName):
"""
Usefull function to plot a confusion matrix
----
INPUT:
cnf_matrix: confusion matrix calculated with scikit
ClaName: string : name of the classifier
"""
import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
class_names=ClassId# [0,1,2,3] # name of classes
class_names=[0,1,2,3] # name of classes
fig, ax = plt.subplots()
tick_marks = np.arange(len(class_names))
plt.xticks(tick_marks, class_names)
plt.yticks(tick_marks, class_names)
# create heatmap
sns.heatmap(pd.DataFrame(cnf_matrix), annot=True, cmap="YlGnBu" ,fmt='g')
#ax.xaxis.set_label_position("top")
ax.xaxis.set_label_position("top")
plt.tight_layout();
plt.title("Confusion matrix using "+ClaName, y=1.1);
plt.ylabel('Actual label');
plt.xlabel('Predicted label');
tick_marks = np.arange(len(class_names))+.5
plt.xticks(tick_marks, class_names)
plt.yticks(tick_marks, class_names)
plt.margins(0.1)
plt.show()
def lazToXYZ(lasFile):
"""
Usefull function to read LAS/LAZ and extract the X,Y,Z position of points
----
INPUT:
lasFile: LAS file name
----
OUTPUT:
X,Y,Z vectors
"""
import laspy
import os
if os.path.isfile(lasFile):
las = laspy.read(lasFile)
x = las.x
y = las.y
z = las.z
return x, y, z
plt.show()
else:
print(f'File {lasFile} is not accessible')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment