# Lab 1 : Extraction de caractéristiques et réduction de dimension
###### tags: 'SYS800'

## 1 Extraction d'attributs (caractéristiques)

On va utiliser deux méthodes d'extraction : classique et automatique.

### 1.1 Extraction manuelle : méthode du codage rétinien

L'bjectif est de coder la méthode d'extraction manuelle sur Matlab.
Celle-ci fonctionne comme suit :

1. Centrer l'image de dimesion $[n\times m]$ dans une matrice carré de dimensions $[max(n,m) \times max(n,m)]$.
2. Il faut ensuite agrandir l'image pour conserver le maximum de détails. La dimension devient $[l*nz \times c*mz]$ avec $l$ le nombre de ligne de la rétine et $c$ le nombre de colonne.
3. Finalement, on extrait de chaque zone la caractéristique correspondante :
$Z(i) = \frac{\sum{Pixels\ de\ la\ zone\ i}}{Nombre\ de\ pixels\ de\ la\ zone\ i}$
4. On conserve chaque valeur dans un vecteur de dimensions $[l*c \times 1]$.

#### Fonction ```retinalcoding()```
Cette fonction prend en entrée trois paramètres :
- `image` : image d'entrée. Matrice de dimension $[n\times m]$.
- `nb_l` : nombre de ligne de la rétine.
- `nb_c` : nombre de colonne de la rétine.

et fournit une sortie :
- `output` : flattened vector/features vector. Vecteur caractéristique de dimensions $[nb_l*c \times 1]$.

La database fournit les images à traiter sous forme de vecteurs. L'image est "applatie".
Pour simplifier la rédaction et la compréhension de l'algorithme, le vecteur est redimensionné en une matrice de dimension 28x28.
L'algorithme permettant de réaliser cette opération est le suivant :

 for n = 1:LDataset
 % go through all image pixel
 
        % sqrtLD corresponds to image dimension (usually 28)
        Mod = mod(n, SqrtLD);
        %used to pass to a new line when dimension has been reached
        
        % case of dimension reached by index
        if Mod == 0
            img(ind, SqrtLD) = tmp(n);
            
            % jump to next line
            ind=ind+1;
            
        % case of index under the dimension value
        else
            img(ind, Mod) = tmp(n);
        end
    end  

Il est maintenant nécessaire de centrer l'image afin que l'ensemble des élements d'une même classe soit le plus similaires possible. Pour obtenir des résultats optimaux, il est important de limiter au maximum la variance entre les élements d'une même classe.
En effet, prennons l'exemple d'un "6". Si celui-ci est collé au bord gauche et qu'un autre est collé au bord droit, leur vecteur rétinien sera impacté par cette position. Centrer les élements permet donc de supprimer cette variance inutile.

L'algorithme réalisant cette operation est le suivant :

% Store size of the image 
[rows, cols] = find(img);

% Set the limit the digit
cropL = min(cols);
cropR = max(cols);
cropT = min(rows);
cropB = max(rows);

% Crop the image to fit exactly the digit
imgc = imcrop(img, [cropL cropT cropR-cropL cropB-cropT]);

% New height and width
height = cropB-cropT +1;
width = cropR-cropL +1;

diff = abs(height - width);
marge = round(diff/2);

% Increase the numbers of columns in order to obtain a square image
if gt(height, width)
    imgNew = zeros(height, height);
    imgNew(: , marge+1:width+marge) = imgc;

% Increase the numbers of lines in order to obtain a square image
elseif gt(width, height)
    imgNew = zeros(width, width);
    imgNew(marge+1:height+marge, :) = imgc;
else
    imgNew = imgc;
end

Maintenant que la base de donnees d'image a été formatée correctement, il faut développer la méthode statistique de codage retinien.

Le code de la fonction est le suivant :

function [output] = retinalcoding(image,nb_l,nb_c)
%% Input Output
% image : image matrix (square) to proceed
% nb_l : retina`s number of lines
% nb_c : retina`s number of columns
% output : final feature vector

%% Step/Etape 1:
% setting the dimensions of the image (image) to (nb_l*M x nb_c*M)

    [rows col] = size(image);
    % nearest to stay in the range 0 253
    imresized = imresize(image,[nb_l*rows nb_c*col], 'nearest');
    
    
%% Step/Etape 2:
% Dividing the image (imresized) into a grid and
% calculating the ratio shape/all in every blocks

    output = zeros(nb_c*nb_c,1);
    for i= 0:1:nb_l-1
        for j = 0:1:nb_c-1
            % Store the current part of the image
            i_temp = imcrop(imresized,[i*rows+1 j*col+1 rows-1 col-1]);
            
            % Compute mean of image part i_temp
            mean = sum(i_temp,'all') / (rows*col*254);
            
            % Store current feature (mean) into the vector
            output(i*nb_c+j+1) = mean;
            
        end
    end

Une base de données de vecteurs caracteristiques propre a chaque image est generée pour les dimensions de retines suivantes :
- 2x2
- 3x3
- 10x10
- 15x15
- 20x20

La comparaison des taux de chevauchements (cf. Preambule) pour chaque type de retine permet de définir la plus efficace.
Comme la dimension des vecteurs caractéristiques augmente fortement avec la dimension de la retine, le temps de calcul du taux de chevauchement est fortement impacté.
Afin de limiter cet effet contraignant tout en obtenant des résultats cohérents, le taux sera calculé en utilisant les 10000 premières images du dataset.
Le tableau suivant repertorie les résultats :

| Dimensions de la retine | 2x2 | 3x3 | 10x10 | 15x15 | 20x20 |
| ------------------------ | --- | --- | --- | --- | --- |
| Nombre d'attributs | 4 | 9 | 100| 225 |400|
| Taux de chevauchement | 0.4935 | 0.1883 | 0.0490| 0.0499 | 0.0620 |

On constate que la retine la plus efficace est celle de dimesion 10x10. Il faut tout de même remarquer que le temps de calcul de la rétine 3x3 est beaucoup plus court que celle de dimension 10x10 au vu du nombre d'attributs.
Le choix de la dimension de la rétine varie en fonction de l'objectif principal recherché : compromis compléxité/précision.
Dans le cadre de ce TP, la rétine de dimension 10x10 sera tout de même considéré comme le choix de la rétine aux meilleurs résultats pour la suite du TP.

#### A.3 ACP
L'objectif de la méthode ACP est de réduire la dimension de l'espace des caractéristiques. Nous l'appliquerons sur les caractéristiques de la base de données obtenue avec la rétine 10x10.
Le script ```acp.m``` est utilisé pour calculer les composantes principales (la matrice V) :

Il prend en entrée :
- Les données X, de dimension ```60 000 x 100```. Nos images y sont stockées horizontalement : une ligne représente une image de 100 pixels
- La valeur ```thres```, un nombre compris entre 0 et 1 qui représente l'énergie cumultaive, soit le pourcentage de variabilité conservée

Il a pour sortie :
- La matrice ```V```, dont les colonnes sont les composantes principales
- La valeur ```G```, qui est le nombre

function [V, G] = acp(X, thres)

% 1) Calcule et soustraction de la moyenne empirique

moyenne_empirique = mean(X);

B = X - ones(size(X,1),1)*moyenne_empirique;

% 2) Calcule de la covariance
C = cov(B);

% 3) Calcule des vecteurs propres et des valeurs propres
[V, D] = eig(C);

lambda = diag(D);

% Triage des vecteurs propres selon l'ordre des valeurs propres du plus
% grand au plus petit
[lambda, ind] = sort(lambda, 'descend');
V = V(ind,:);

% Calcule du contenue d'energy cumulative de chaque vecteurs propre
G = cumsum(lambda)/sum(lambda)*100;

% Seuillage juste audessus du seuil%
nb_dim = length(G) - length(G(G >= thres)) + 1;

% Reduction du nombre de dimension
V = V(:,1:nb_dim);

La fonction ```mean(X)``` calcule d'abord la moyenne de chaque "axe". On obtient 100 moyennes, qui correspondent aux valeurs moyennes de chaque pixel. On soustrait ensuite ces moyennes aux colonnes correspondantes de la matrice ```X```, ce qui permet de centrer les axes par rapport aux données, pour le point de vue géometrique.
On calcule ensuite la covariance à l'aide de la fonction ```cov(B)``` qui retourne le produit $B^TB$.
La matrice ```C```, de dimension 100x100, représente ainsi la covariance entre les "axes" : on les compare les un avec les autres pour ensuite maximiser cette covariance.
Pour cela on utilise les valeurs propres obtenus avec la fonction ```eig(C)```.
On trie ensuite ces valeurs propres dans l'ordre décroissant et on les range dans une matrice diagonale.
En fonction de l'énergie qui a été prédéfinie, on calcule le nombre de dimensions jugées pertinentes (```nb_dim```). On se débarasse du reste en conservant uniquement les ```nb_dim``` premières colonnes de la matrice diagonale ```V``` contenant les valeurs propres ordonnées.
Cette nouvelle matrice ```V``` contient ainsi un nombre réduit de dimensions, dont la covariance est maximale et donc la plus pertinente possible pour représenter nos données.
### A.4 Projection de la base de test
On peut maintenant utiliser la matrice ```V``` pour projeter nos données de test et ainsi réduire la dimensionnalité.
Pour ce faire, on utilise le script ```projection_acp``` :

function [T] = projection_acp(X, V)

% 1) Calcule et soustraction de la moyenne empirique
moyenne_empirique = mean(X);

% 2) Normalization des données
B = X - ones(size(X,1),1)*moyenne_empirique;

% 3) Projections des données
T = B*V;

Il répéte simplement la méthode pour centrer les données, puis les projète avec la matrice ```V```. Avec une énergie fixée à 95%, le nombre de caractéristiques passe de 100 (10x10 pixels) à 56.

### 1.2 Extraction automatique : Autoencodeur

Cette fois-ci, on va utiliser une approche réseau de neuronnes.
Le fonctionnement d'un auto-encodeur est le suivant :
1. **Encodeur :** compresse l'entrée (ici l'image) en un code (réprésentation latente) dans un espace de dimension réduit appelé espace latent.
2. **Code :** Couche contenant le code généré représentant l'entrée dans l'espace réduit.
3. **Décodeur :** reconstruit du mieux possible l'image d'entrée en se prenant comme entrée le code générée.

L'autoencodeur est ainsi une méthode non supervisée ayant pour objectif de reconstruire un ensemble d'objet en utilisant leur représentation dans l'espace latent, leur "code".

Les autoencodeurs qui seront ont été mis en place sous matlab afin de réutiliser facilement les fonctions de préprocessing déjà développées.

#### 1.2.1 Autoencodeur à une couche
##### Définition de l'autoencodeur

Le premier autoencodeur mis en place est composé d'une seule couche.
L'image suivante présente son format. Les dimensions indiquées varient en fonction des attributs définit.


Son code est le suivant :

%% Step One : Preprocesing

% Load training dataset 
mnisttrain = csvread('mnist_train.csv');
[nb_train lg_images_train] = size(mnisttrain);
% nb_train = 5000; 
% can be hardly fixed to accelerate the training for trial
Xtrain = preprocessing(mnisttrain(1:nb_train, 2:785),nb_train);
train_labels = mnisttrain(:, 1);

% Load test dataset 
mnisttest = csvread('mnist_test.csv');
[nb_test lg_images_test] = size(mnisttest);
Xtest = preprocessing(mnisttest(1:nb_test, 2:785),nb_test);
test_labels = mnisttest(:, 1);

%% Step Two : Autoencoder Training and testing
hiddenSize = 25;

autoenc = trainAutoencoder(Xtrain,hiddenSize, 'L2WeightRegularization',0.004,'DecoderTransferFunction','purelin','SparsityRegularization',4,'SparsityProportion',0.15);

xReconstructed = predict(autoenc,Xtest);


%% Step Three :  Plot tested pictures

% Reset image in the integer space
Xtest_int = cell(1,nb_test);
for i = 1:nb_train
    Xtest_int(1,i) = {uint8(255 * mat2gray(cell2mat(Xtest(1,i))))};
end

% Plot 20 first images
figure;
for i = 1:20
    subplot(4,5,i);
    imshow(Xtest_int{i});
end


%% Plot reconstructed pictures

% Reset image in the integer space
xReconstructed_int = cell(1,nb_test);
for i = 1:nb_train
    xReconstructed_int(1,i) = {uint8(255 * mat2gray(cell2mat(xReconstructed(1,i))))};
end

% Plot 20 first images
figure;
for i = 1:20
    subplot(4,5,i);
    imshow(xReconstructed_int{i});
end

La première partie du code (`%% Step One`) se concentre sur le formatage des données.

Dans la même idée que pour le traitement rétinien, les images sont centrées. Une étape de normalisation est mise en place afin de permettre une meilleur convergence du gradient.

En effet, sans la normalization, le gradient des poids liés à une caractéristique sera beaucoup plus important que le gradient d'une autre caractéristique pendant la formation, puisque les gradients dépendent de l'entrée.

Une normalisation min = 0, max = 1 a été choisie.
Cette étape est réalisée dans la fonction `preprocessing()` décrite ci-dessous :

function [database] = preprocessing(dataset,nb_images)
% Format test and train image for autoencodeur
% INPUT :
%
% dataset : Images flatten in row vectors
% nb_images : nb images to prepocess from the dataset
% OUTPUT :
%
% database : cell of images preprocessed


LDataset = length(dataset(1,:));
SqrtLD = sqrt(LDataset);

database = cell(1,nb_images);

img = zeros(SqrtLD, SqrtLD);

for i = 1:1:nb_images  
    ind = 1; 
    tmp = dataset(i, :);

    for n = 1:LDataset
        Mod = mod(n, SqrtLD);
        if Mod == 0
            img(ind, SqrtLD) = tmp(n);
            ind=ind+1;
        else
            img(ind, Mod) = tmp(n);
        end
    end
    ind = 1;
    
    % normalization to increase autoencodeur's performanceÒ
    ct_nm_im = im2double(CenterImage(img));
    database(1,i) = {ct_nm_im};
end

La seconde partie (`%%Step Two`) entraine l'autoencodeur. Les hyperparamètres sont d'abord définit de la même manière que le propose la documentation MatLab. Ceux-si seront par la suite modifier pour améliorer les performances de l'autoencodeur.

###### Hyperparamètres
Un méthode de régularisation permettant d'éviter un sur-apprentissage est ajouté. La technique retenue, qui est la plus utilisé est la norme L2. La formule mathématique est la suivante :
$$Loss = Error(Y - \widehat{Y}) + fact{\lambda}{2} \sum_1^n w_i^{2}$$
Le facteur $\lambda$ est définit à 0,004. Ce facteur définit l'importance de la régularisation, soit son impact sur la loss function.

La dernière partie du code (`%% Step Three`), permet d'afficher les images testées et les images reconstruites correspondantes, et de calculer l'erreur de de reconstruction.

##### Résultats et analyses

Pour définir la performance de l'autoencodeur, l'erreur quadratique moyenne entre les images tests et les images reconstruites est calculé. La formule est la suivante :
$$ MSE = \sum_{i=1}^{D}(x_i-y_i)^2 $$
avec D la nombre de pixel présents dans l'image.

Dans le cas de l'autoencodeur présenté dans cette partie, avec 25 features, le tracé de l'erreur selon les epochs est la suivante :
[](https://i.imgur.com/LkoqiS6.png)
*Figure i : Erreur de reconstruction selon les epochs | 25 features*

Les figures suivantes présentese les images tests et les images reconstruites.

*Figure X : Images test avec un autoencodeur à une couche | 25 features*

*Figure X : Images reconstruites avec un autoencodeur à une couche| 25 features*

Résultats avec 32 caractéristiques extraites :
*Figure i : Erreur de reconstruction selon les epochs | 32 features*

*Figure i : Erreur de reconstruction selon les epochs | 64 features*

*Figure i : Erreur de reconstruction selon les epochs | 128 features*

On constate que les résultats sont meilleurs avec une dimenison plus importante, cela est aussi visible en observant les images reconstituées.
*Figure X : Images tests (gauche) et reconstruites (droite) avec un autoencodeur à une couche *