Delete 3dlidar_classif.py

Labs/Lab_3_LiDAR_ML_Segmentation/3dlidar_classif.py

-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Jan 15 22:19:54 2021
-
-@author: alucas
-"""
-
-
-
-import numpy as np
-from numpy import savetxt
-
-from scipy import spatial
-from scipy import spatial
-from scipy.stats import gaussian_kde
-
-import pandas as pd
-
-import matplotlib.pyplot as plt
-import matplotlib.tri as tri
-from matplotlib.colors import ListedColormap
-from mpl_toolkits.mplot3d import Axes3D
-
-import seaborn as sns
-
-from sklearn.decomposition import PCA
-from sklearn.model_selection import train_test_split
-from sklearn.tree import DecisionTreeClassifier
-from sklearn.metrics import confusion_matrix
-from sklearn.metrics import pairwise_distances
-
-
-
-
-def readData(filename):
-    """
-        Function to read the landscape data
-        INPUT:
-        ------
-            @filename: string
-
-        OUTPUT:
-        -------
-            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]  
-        y = floor[:,1]
-        z = floor[:,2]
-        f.close()
-        
-        return x, y, z
-        
-    else:
-        print(f'File {filename} is not accessible')
-    
-    
-    
-    
-def plot_figs(fig_num, elev, azim, dx, dy, dz, density=1):
-    """
-    Draw (x, y z) coordinates into a 3 dimensional scatter plot.
-    INPUT:
-    ------
-        @fig_num: int   -  Figure number
-        @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 
-    
-    """
-   
-    # Subsample input data 
-    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()
-    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 getEigenvaluePCA(a, b, c, dim, decim=1):
-    """
-    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.
-    
-    """
-    print('Data are decimed by a factor of ',decim)
-    ddComp = []
-    Y = np.c_[a[::decim], b[::decim], c[::decim]]
-    #data = Y.copy()
-    Y = np.float32(Y)
-    #dist = np.float32(spatial.distance.pdist(Y))
-    dist = spatial.distance.squareform(spatial.distance.pdist(Y))    
-    a -= np.mean(a, axis=0)
-    b -= np.mean(b, axis=0)
-    c -= np.mean(c, axis=0)
-    
-    ##Y = np.c_[a[::decim], b[::decim], c[::decim]] <<<--- pourquoi une nouvelle fois le decim' ? on ne devrait pas
-
-    n_samples = Y.shape[0]
-    feat = []   
-    
-    print('  Treating distance of', dim, 'm')
-    feat = []
-    comp = []
-
-    # Loop over every points
-    for kk in range(0,n_samples):
-        pts = np.where((dist[kk,:]<=dim/2))
-
-        if np.size(pts)<2:
-            print("Revise dims")
-            print(pts)
-            print(dist)
-
-        Ytmp = Y[pts,:]
-        Ytmp = Ytmp[0,:,:]
-
-        pca = PCA(n_components=3)
-        pca.fit(Ytmp)
-        eigenvalues = pca.explained_variance_ratio_
-        comp.append(np.array(eigenvalues))
-    comp = np.array(comp)
-    return comp,Y
-
-
-
-
-
-def estimateTernaryCoord(comp):       
-    """
-    Function which estimate the ternary coordinates 
-    from the list of eigenvalues.
-    ----
-    INPUT:
-        @comp: np array - list of eigenvalue for every 
-                points with a given dimension
-    OUTPUT:
-    -------
-        @x_ter, @y_ter, @z_ter
-    
-    """
-    # Conversion towards a,b,c space  
-    ct = 3 * comp[:,2]
-
-    xp1p2 = (1 - (comp[:,1] - comp[:,0]))/2
-    yp1p2 = 1 - xp1p2
-
-    dS2 = np.sqrt((.5 -  xp1p2)**2 + (.5 - yp1p2)**2)
-    
-    at = (1 - ct)*dS2 /(np.sqrt((0.5 -  1)**2 + (0.5 - 0)**2))
-    bt = 1 - (ct + at) 
-    v = (at + bt + ct)  
-
-    # Conversion towards ternary graph 
-    x_ter = .5 * ( 2.*bt+ct ) / v        
-    y_ter = 0.5*np.sqrt(3) * ct / v
-
-    # Compute density for clarity
-    xy = np.vstack([x_ter, y_ter])
-    z_ter = gaussian_kde(xy)(xy)
-    
-    return x_ter, y_ter, z_ter 
-
-
-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()
-    
-    
-    
-    
-    
-workingdir='./LiDARDunes/training/'
-filename=workingdir+'floor.xyz'
-dx1, dy1, dz1 = readData(filename)
-
-dims = [0.1, 0.3]
-floor = pd.DataFrame()
-
-for dim in dims:
-    eigv,Y1 = getEigenvaluePCA(dx1, dy1, dz1, dim, decim=5)
-    x_ter, y_ter, z_ter = estimateTernaryCoord(eigv)
-    data_dim = pd.DataFrame({'x_ter_'+str(dim):x_ter, 
-                             'y_ter_'+str(dim):y_ter, 
-                             'z_ter_'+str(dim):z_ter})
-    
-    floor = pd.concat([floor, data_dim], axis=1)
-    title = "Components at "+str(dim)+" m"
-    ##plotTernary(x_ter, y_ter, z_ter, title)
-floor["class"] = 1
-
-del dx1, dy1, dz1 
-
-workingdir='./LiDARDunes/training/'
-filename=workingdir+'vegetation.xyz'
-dx, dy, dz = readData(filename)
-
-# All together
-dims = [0.1, 0.3]
-veget = pd.DataFrame()
-
-for dim in dims:
-    eigv,Y2 = getEigenvaluePCA(dx, dy, dz, dim, decim=5)
-    x_ter, y_ter, z_ter = estimateTernaryCoord(eigv)
-    data_dim = pd.DataFrame({'x_ter_'+str(dim):x_ter, 
-                             'y_ter_'+str(dim):y_ter, 
-                             'z_ter_'+str(dim):z_ter})
-    
-    veget = pd.concat([veget, data_dim], axis=1)
-    title = "Components at "+str(dim)+" m"
-    #plotTernary(x_ter, y_ter, z_ter, title)
-veget["class"] = 2
-del dx, dy, dz,data_dim,eigv,x_ter,y_ter,z_ter 
-
-
-
-data = pd.DataFrame()
-data = pd.concat([veget, floor])
-del veget,floor
-#X = data.drop("class", axis=1)
-X = data.drop({"z_ter_0.1","z_ter_0.3","class"}, axis=1)
-y = data["class"].copy()
-X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)
-
-
-from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
-#from sklearn.metrics import plot_confusion_matrix
-ClaName="LDA"
-lda = LinearDiscriminantAnalysis(store_covariance=True)
-lda.fit(X_train, y_train)
-y_pred_dt = lda.predict(X_test)
-cnf_matrix=confusion_matrix(y_test, y_pred_dt)
-#plot_confMat(cnf_matrix,ClaName)
-y_pred_lda = lda.predict(X_test)
-y_pred_all = lda.predict(X)
-y_pred_train = lda.predict(X_train)
-y_pred = lda.fit(X, y.ravel()).predict(X)
-print("Classifier: Decision Tree")
-print("Accuracy on the train data: {:.2f}".format(lda.score(X_train, y_pred_train)*100)+"%")
-print("Accuracy on test data: {:.2f}".format(lda.score(X_test, y_pred_lda)*100)+"%")
-print("Accuracy on the whole data: {:.2f}".format(lda.score(X, y_pred_all)*100)+"%")
-print(" ")
-
-
-
-
-
-y_scene = lda.predict(X)
-Y = np.concatenate((Y2, Y1))
-
-del Y1,Y2
-#dy = np.concatenate((dy1, dy))
-#dz = np.concatenate((dz1, dz))
-    
-plot_3dcladd(Y[:,0],Y[:,1],Y[:,2],y_scene,density=1)
-
-
-####
-
-workingdir='./LiDARDunes/data/'
-filename=workingdir+'scene2.xyz'
-dx, dy, dz = readData(filename)
-
-
-#Y = np.c_[dx[::5],dy[::5],dz[::5]]
-#dist = pairwise_distances(Y)
-
-# # All together
-dims = [0.1, 0.3]
-scene = pd.DataFrame()
-
-for dim in dims:
-     eigv,YY = getEigenvaluePCA(dx, dy, dz, dim, decim=5)
-     x_ter, y_ter, z_ter = estimateTernaryCoord(eigv)
-     data_dim = pd.DataFrame({'x_ter_'+str(dim):x_ter, 
-                               'y_ter_'+str(dim):y_ter, 
-                               'z_ter_'+str(dim):z_ter})
-    
-     scene = pd.concat([scene, data_dim], axis=1)
-#     title = "Components at "+str(dim)+" m"
-#     #plotTernary(x_ter, y_ter, z_ter, title)
-
-
-
-# y_scene = lda.predict(scene)
-
-# plot_3dcladd(YY[:,0],YY[:,1],YY[:,2],y_scene,density=1)