Date

Last year, we have published a paper on feature extraction with Gaussian Mixture Model (see pdf). This work has been continued these past months by Adrien Lagrange during its internship at the lab, under co-supervision of Manuel Grizonnet (CNES).

We have extended the linear algebra to save computational resources during the main step of the optimization. A paper has been submitted and details of the computation can be seen on HAL. As usual, the python code is freely available from the github of Adrien.

Furthermore, the algorithm have been implemented in an Orfeo Toolbox module. It is now available from the external modules webpage.

In this post, I am going to illustrate how it works and behaves on the famous AVIRIS data set. See images falseColor and ref. The image has 220 spectral bands and a size of 145 \(\times\) 145 pixels. The original data can be download:

Original false colors images.
Original false colors images.
Reference data
Reference data

How it works ?

otbcli_TrainGMMSelectionApp -help
This is the TrainGMMSelectionApp application, version 5.4.0
Train a GMM classifier from multiple pairs of images and training vector data.

Complete documentation: http://www.orfeo-toolbox.org/Applications/TrainGMMSelectionApp.html

Parameters: 
        -progress      <boolean>        Report progress 
MISSING -io.il         <string list>    Input Image List  (mandatory)
MISSING -io.vd         <string list>    Input Vector Data List  (mandatory)
        -io.imstat     <string>         Input XML image statistics file  (optional, off by default)
        -io.confmatout <string>         Output confusion matrix  (optional, off by default)
MISSING -io.out        <string>         Output model  (mandatory)
        -sample.mt     <int32>          Maximum training sample size per class  (mandatory, default value is 1000)
        -sample.mv     <int32>          Maximum validation sample size per class  (mandatory, default value is 1000)
        -sample.bm     <int32>          Bound sample number by minimum  (mandatory, default value is 1)
        -sample.edg    <boolean>        On edge pixel inclusion  (optional, off by default)
        -sample.vtr    <float>          Training and validation sample ratio  (mandatory, default value is 0.5)
        -sample.vfn    <string>         Name of the discrimination field  (mandatory, default value is Class)
MISSING -gmm.varnb     <int32>          Number of variables to select  (mandatory)
        -gmm.method    <string>         Method used for selection [forward/sffs] (mandatory, default value is forward)
        -gmm.crit      <string>         Criterion function for selection [jm/divkl/accuracy/kappa/f1mean] (mandatory, default value is jm)
        -gmm.ncv       <int32>          Number of folds for cross-validation  (mandatory, default value is 5)
        -gmm.best      <int32>          If 1, choose optimal set of features based on criterion function after selection. If 0, all selected features are used  (mandatory, default value is 1)
        -gmm.seed      <int32>          Rand seed for cross-validation  (mandatory, default value is 0)
        -rand          <int32>          set user defined seed  (optional, off by default)
        -inxml         <string>         Load otb application from xml file  (optional, off by default)

Examples: 
otbcli_TrainGMMSelectionApp -io.il QB_1_ortho.tif -io.vd VectorData_QB1.shp -io.imstat EstimateImageStatisticsQB1.xml -sample.mv 100 -sample.mt 100 -sample.vtr 0.5 -sample.edg false -sample.vfn Class -gmm.varnb 20 -gmm.method forward -gmm.crit jm -gmm.ncv 5 -gmm.best 1 -gmm.seed 0 -io.out svmModelQB1.txt -io.confmatout svmConfusionMatrixQB1.csv

So we need an image, a reference file and to choose what is the method we are going to use. I refer to the aforementioned paper for details ;)

Here, we are going to apply the feature selection algorithm with the following parameters:

  • Criteria: the Kappa coefficient and the Jeffreys-Matusita distance.
  • Method: forward.
  • Number of variables: 30
  • Number of samples per class: 50

Let’s go:

for crit in jm kappa
do
otbcli_TrainGMMSelectionApp -io.il ./images/aviris.tif \
                -io.vd ./images/ref.shp \
                -io.confmatout ./images/res_simu_npfs/confu_${crit}.txt \
                -io.out ./images/res_simu_npfs/model_${crit}.txt \
                            -sample.mt 50 \
                -sample.bm 0 \
                -gmm.method forward \
                -gmm.crit ${crit} \
                -gmm.varnb 30 \
                -rand 0    
done 

Then we can extract the number of selected features, and the global accuracy, for each criteria, from the txt files:

import linecache
import scipy as sp

for name_ in ["kappa","jm"]:
    # Selected bands
    filename = "images/res_simu_npfs/model_{}_selection.txt".format(name_)
    NumberOfSelectedBands = int(linecache.getline(filename,2).split()[1])
    SelectedBands = linecache.getline(filename,2+NumberOfSelectedBands).split()
    print("Number of selected band for {} criteria: {}".format(name_,NumberOfSelectedBands))
    print("Selected Bands:")
    print(SelectedBands)
    print("")
Number of selected band for kappa criteria: 8
Selected Bands:
['27', '35', '74', '133', '40', '203', '10']

Number of selected band for jm criteria: 30
Selected Bands:
['180', '131', '71', '12', '38', '34', '167', '61', '119', '47', '89', '198', '26', '6', '4', '139', '116', '144', '96', '32', '163', '197', '146', '46', '28', '130', '135', '30', '59']

import linecache
import scipy as sp
import matplotlib.pyplot as plt

for name_ in ["kappa","jm"]:
    # Confusion Matrix
    filename = "images/res_simu_npfs/confu_{}.txt".format(name_)
    confusionMatrix = sp.genfromtxt(filename,delimiter=',',skip_header=2)
    OA = sp.diag(confusionMatrix).sum()/confusionMatrix.sum()
    print("Overall Accuracy for criteria {}: {}".format(name_,OA))
Overall Accuracy for criteria kappa: 0.81195335277
Overall Accuracy for criteria jm: 0.607871720117

From the classification accuracy, we can see that the JM criteria performs worse than the Kappa criteria (similar conclusion would have been done if F1 mean or OA criterion have been used). JM is widely used in remote sensing, but it might perform badly in some cases where the Gaussianity is not fulfilled, as illustrated in these short example.

It is also possible to get the learning curves, i.e., the Kappa/JM value in function of the number of selected bands.

import linecache
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib
matplotlib.style.use('ggplot')

for name_ in ["kappa","jm"]:
    # Selected bands
    filename = "images/res_simu_npfs/model_{}_selection.txt".format(name_)
    crit = []
    for lineNumber in range(35,65):
        crit.append(linecache.getline(filename,lineNumber).split()[1])
    plt.figure()
    plt.title("Learning curves for {}".format(name_))
    plt.plot(range(1,31),crit,'o')
    plt.grid(b=True)
    plt.savefig("images/learning_curve_{}.png".format(name_))
    print("[[./images/learning_curve_{}.png]]".format(name_))

A global maximum is reached with the Kappa for 8 extracted features. If we look to the table tab, this is the highest number of features for which the number of training samples is higher than the number of parameters to estimate.

For the JM criteria, the value still increased (slowly) for 30 features. /Note: rather than selecting the maximum value of the criteria, we could stop when the increased of the criteria is too small. It can be done in two steps with the module./

Number of variables (nv) to be estimated in function of the number of features (nf): nv = nf*(nf+1)/2 + nf (proportion is omitted).
# of variables # of parameters to estimate
1 2
2 5
3 9
4 14
5 20
6 27
7 35
8 44
9 54
10 65

Classify the whole data set

Now it is possible to classify the whole image:

otbcli_PredictGMMApp -help
This is the PredictGMMApp application, version 5.4.0
Performs a classification of the input image according to a GMM model file.

Complete documentation: http://www.orfeo-toolbox.org/Applications/PredictGMMApp.html

Parameters: 
        -progress  <boolean>        Report progress 
MISSING -in        <string>         Input Image  (mandatory)
        -mask      <string>         Input Mask  (optional, off by default)
MISSING -model     <string>         Model file  (mandatory)
        -modeltype <string>         Type of GMM model [basic/selection] (mandatory, default value is basic)
        -varnb     <int32>          New number of variables for prediction. (Ignore number of variables chosen during training)  (optional, off by default)
        -imstat    <string>         Statistics file  (optional, off by default)
MISSING -out       <string> [pixel] Output Image  [pixel=uint8/uint16/int16/uint32/int32/float/double] (default value is uint8) (mandatory)
        -confmap   <string> [pixel] Confidence map  [pixel=uint8/uint16/int16/uint32/int32/float/double] (default value is double) (optional, off by default)
        -ram       <int32>          Available RAM (Mb)  (optional, off by default, default value is 128)
        -inxml     <string>         Load otb application from xml file  (optional, off by default)

Examples: 
otbcli_PredictGMMApp -in QB_1_ortho.tif -imstat EstimateImageStatisticsQB1.xml -model model.txt -modeltype basic -out clLabeledImageQB1.tif

By default, the number of selected for the prediction corresponds to the number of variables for which the criteria is maximum. But it is possible to change it with the option -varnb. To classify the image using the model obtain with Kappa

otbcli_PredictGMMApp -in images/aviris.tif \
             -model images/res_simu_npfs/model_kappa.txt \
             -modeltype selection \
             -out images/res_simu_npfs/thematic_aviris.tif

The resulting images is given in fig:thematic. Off course, spatial information should be included in the classification process. But the subject of this post is about feature extraction: with only 8 features from 220, and 50 training samples per class we get an OA of 81 % … so that’s all folk !

Aviris thematic map
Aviris thematic map

And with SVM ?

If we do the classification with the classifier SVM ? Using the default implementation of OTB (LIBSVM), it can be done with the following command:

otbcli_ComputeImagesStatistics -il ./images/aviris.tif -out stat.xml
otbcli_TrainImagesClassifier -io.il ./images/aviris.tif \
                 -io.vd ./images/ref.shp -io.out .images/res_simu_npfs/model_svm.txt \
                 -io.imstat stat.xml \
                 -classifier libsvm \
                 -classifier.libsvm.k rbf -classifier.libsvm.opt 1 \
                 -rand 0 -sample.mt 50 -sample.bm 0 \
                 -io.confmatout .images/res_simu_npfs/confu_svm.txt

Then the OA is again found by

import scipy as sp
filename = "images/res_simu_npfs/confu_svm.txt"
confusionMatrix = sp.genfromtxt(filename,delimiter=',',skip_header=2)
OA = sp.diag(confusionMatrix).sum()/confusionMatrix.sum()
print("Overall Accuracy for SVM: {}".format(OA))
Overall Accuracy for SVM: 0.711370262391

What ? We are better than the super SVM ! Ok, it’s my website … But if you check the JSTARS paper, when the number of samples is very low, we are usually better. But when the number of training samples becomes large, then SVM (and RF) tends to be better in terms of classification accuracy.


Comments

comments powered by Disqus