## Chapter 19Classification

### 19.1 Introduction

In statistical classification, each object is represented by d features (a measurement vector), and the goal of classification becomes finding compact and disjoint regions (decision regions[42]) for classes in a d-dimensional feature space. Such decision regions are defined by decision rules that are known or can be trained. The simplest configuration of a classification consists of a decision rule and multiple membership functions; each membership function represents a class. Figure 19.1 illustrates this general framework.

This framework closely follows that of Duda and Hart[42]. The classification process can be described as follows:

1. A measurement vector is input to each membership function.
2. Membership functions feed the membership scores to the decision rule.
3. A decision rule compares the membership scores and returns a class label.

This simple configuration can be used to formulated various classification tasks by using different membership functions and incorporating task specific requirements and prior knowledge into the decision rule. For example, instead of using probability density functions as membership functions, through distance functions and a minimum value decision rule (which assigns a class from the distance function that returns the smallest value) users can achieve a least squared error classifier. As another example, users can add a rejection scheme to the decision rule so that even in a situation where the membership scores suggest a “winner”, a measurement vector can be flagged as ill defined. Such a rejection scheme can avoid risks of assigning a class label without a proper win margin.

Note that to use these concepts into your own programs, you might have to link to other OTB libraries (edit the TARGET_LINK_LIBRARIES to add them), for example OTBLearning, OTBMarkov.

#### 19.1.1 k-d Tree Based k-Means Clustering

The source code for this example can be found in the file
Examples/Classification/KdTreeBasedKMeansClustering.cxx.

K-means clustering is a popular clustering algorithm because it is simple and usually converges to a reasonable solution. The k-means algorithm works as follows:

1. Obtains the initial k means input from the user.
2. Assigns each measurement vector in a sample container to its closest mean among the k number of means (i.e., update the membership of each measurement vectors to the nearest of the k clusters).
3. Calculates each cluster’s mean from the newly assigned measurement vectors (updates the centroid (mean) of k clusters).
4. Repeats step 2 and step 3 until it meets the termination criteria.

The most common termination criteria is that if there is no measurement vector that changes its cluster membership from the previous iteration, then the algorithm stops.

The itk::Statistics::KdTreeBasedKmeansEstimator is a variation of this logic. The k-means clustering algorithm is computationally very expensive because it has to recalculate the mean at each iteration. To update the mean values, we have to calculate the distance between k means and each and every measurement vector. To reduce the computational burden, the KdTreeBasedKmeansEstimator uses a special data structure: the k-d tree ( itk::Statistics::KdTree) with additional information. The additional information includes the number and the vector sum of measurement vectors under each node under the tree architecture.

With such additional information and the k-d tree data structure, we can reduce the computational cost of the distance calculation and means. Instead of calculating each measurement vectors and k means, we can simply compare each node of the k-d tree and the k means. This idea of utilizing a k-d tree can be found in multiple articles [6] [107] [78]. Our implementation of this scheme follows the article by the Kanungo et al [78].

We use the itk::Statistics::ListSample as the input sample, the itk::Vector as the measurement vector. The following code snippet includes their header files.

#include "itkVector.h"
#include "itkListSample.h"

Since this k-means algorithm requires a itk::Statistics::KdTree object as an input, we include the KdTree class header file. As mentioned above, we need a k-d tree with the vector sum and the number of measurement vectors. Therefore we use the itk::Statistics::WeightedCentroidKdTreeGenerator instead of the itk::Statistics::KdTreeGenerator that generate a k-d tree without such additional information.

#include "itkKdTree.h"
#include "itkWeightedCentroidKdTreeGenerator.h"

The KdTreeBasedKmeansEstimator class is the implementation of the k-means algorithm. It does not create k clusters. Instead, it returns the mean estimates for the k clusters.

#include "itkKdTreeBasedKmeansEstimator.h"

To generate the clusters, we must create k instances of itk::Statistics::EuclideanDistance function as the membership functions for each cluster and plug that—along with a sample—into an itk::Statistics::SampleClassifier object to get a itk::Statistics::MembershipSample that stores pairs of measurement vectors and their associated class labels (k labels).

#include "itkMinimumDecisionRule.h"
#include "itkEuclideanDistance.h"
#include "itkSampleClassifier.h"

We will fill the sample with random variables from two normal distribution using the itk::Statistics::NormalVariateGenerator.

#include "itkNormalVariateGenerator.h"

Since the NormalVariateGenerator class only supports 1-D, we define our measurement vector type as one component vector. We then, create a ListSample object for data inputs. Each measurement vector is of length 1. We set this using the SetMeasurementVectorSize() method.

typedef itk::Vector<double,
1>
MeasurementVectorType;
typedef itk::Statistics::ListSample<MeasurementVectorType> SampleType;
SampleType::Pointer sample = SampleType::New();
sample->SetMeasurementVectorSize(1);

The following code snippet creates a NormalVariateGenerator object. Since the random variable generator returns values according to the standard normal distribution (The mean is zero, and the standard deviation is one), before pushing random values into the sample, we change the mean and standard deviation. We want two normal (Gaussian) distribution data. We have two for loops. Each for loop uses different mean and standard deviation. Before we fill the sample with the second distribution data, we call Initialize(random seed) method, to recreate the pool of random variables in the normalGenerator.

To see the probability density plots from the two distribution, refer to the Figure 19.3.

typedef itk::Statistics::NormalVariateGenerator NormalGeneratorType;
NormalGeneratorType::Pointer normalGenerator = NormalGeneratorType::New();

normalGenerator->Initialize(101);

MeasurementVectorType mv;
double                mean = 100;
double                standardDeviation = 30;
for (unsigned int i = 0; i < 100; ++i)
{
mv[0] = (normalGenerator->GetVariate()  standardDeviation) + mean;
sample->PushBack(mv);
}

normalGenerator->Initialize(3024);
mean = 200;
standardDeviation = 30;
for (unsigned int i = 0; i < 100; ++i)
{
mv[0] = (normalGenerator->GetVariate()  standardDeviation) + mean;
sample->PushBack(mv);
}

We create a k-d tree.

typedef itk::Statistics::WeightedCentroidKdTreeGenerator<SampleType>
TreeGeneratorType;
TreeGeneratorType::Pointer treeGenerator = TreeGeneratorType::New();

treeGenerator->SetSample(sample);
treeGenerator->SetBucketSize(16);
treeGenerator->Update();

Once we have the k-d tree, it is a simple procedure to produce k mean estimates.

We create the KdTreeBasedKmeansEstimator. Then, we provide the initial mean values using the SetParameters(). Since we are dealing with two normal distribution in a 1-D space, the size of the mean value array is two. The first element is the first mean value, and the second is the second mean value. If we used two normal distributions in a 2-D space, the size of array would be four, and the first two elements would be the two components of the first normal distribution’s mean vector. We plug-in the k-d tree using the SetKdTree().

The remaining two methods specify the termination condition. The estimation process stops when the number of iterations reaches the maximum iteration value set by the SetMaximumIteration(), or the distances between the newly calculated mean (centroid) values and previous ones are within the threshold set by the SetCentroidPositionChangesThreshold(). The final step is to call the StartOptimization() method.

The for loop will print out the mean estimates from the estimation process.

typedef TreeGeneratorType::KdTreeType                         TreeType;
typedef itk::Statistics::KdTreeBasedKmeansEstimator<TreeType> EstimatorType;
EstimatorType::Pointer estimator = EstimatorType::New();

EstimatorType::ParametersType initialMeans(2);
initialMeans[0] = 0.0;
initialMeans[1] = 0.0;

estimator->SetParameters(initialMeans);
estimator->SetKdTree(treeGenerator->GetOutput());
estimator->SetMaximumIteration(200);
estimator->SetCentroidPositionChangesThreshold(0.0);
estimator->StartOptimization();

EstimatorType::ParametersType estimatedMeans = estimator->GetParameters();

for (unsigned int i = 0; i < 2; ++i)
{
std::cout << "cluster[" << i << "] " << std::endl;
std::cout << "    estimated mean : " << estimatedMeans[i] << std::endl;
}

If we are only interested in finding the mean estimates, we might stop. However, to illustrate how a classifier can be formed using the statistical classification framework. We go a little bit further in this example.

Since the k-means algorithm is an minimum distance classifier using the estimated k means and the measurement vectors. We use the EuclideanDistance class as membership functions. Our choice for the decision rule is the itk::Statistics::MinimumDecisionRule that returns the index of the membership functions that have the smallest value for a measurement vector.

After creating a SampleClassifier object and a MinimumDecisionRule object, we plug-in the decisionRule and the sample to the classifier. Then, we must specify the number of classes that will be considered using the SetNumberOfClasses() method.

The remainder of the following code snippet shows how to use user-specified class labels. The classification result will be stored in a MembershipSample object, and for each measurement vector, its class label will be one of the two class labels, 100 and 200 (unsigned int).

typedef itk::Statistics::EuclideanDistance
<MeasurementVectorType> MembershipFunctionType;
typedef itk::MinimumDecisionRule DecisionRuleType;
DecisionRuleType::Pointer decisionRule = DecisionRuleType::New();

typedef itk::Statistics::SampleClassifier<SampleType> ClassifierType;
ClassifierType::Pointer classifier = ClassifierType::New();

classifier->SetDecisionRule((itk::DecisionRuleBase::Pointer) decisionRule);
classifier->SetSample(sample);
classifier->SetNumberOfClasses(2);

std::vector<unsigned int> classLabels;
classLabels.resize(2);
classLabels[0] = 100;
classLabels[1] = 200;

classifier->SetMembershipFunctionClassLabels(classLabels);

The classifier is almost ready to do the classification process except that it needs two membership functions that represents two clusters respectively.

In this example, the two clusters are modeled by two Euclidean distance functions. The distance function (model) has only one parameter, its mean (centroid) set by the SetOrigin() method. To plug-in two distance functions, we call the AddMembershipFunction() method. Then invocation of the Update() method will perform the classification.

std::vector<MembershipFunctionType::Pointer> membershipFunctions;

MembershipFunctionType::OriginType origin(
sample->GetMeasurementVectorSize());
int index = 0;
for (unsigned int i = 0; i < 2; ++i)
{
membershipFunctions.push_back(MembershipFunctionType::New());
for (unsigned int j = 0; j < sample->GetMeasurementVectorSize(); ++j)
{
origin[j] = estimatedMeans[index++];
}
membershipFunctions[i]->SetOrigin(origin);
}

classifier->Update();

The following code snippet prints out the measurement vectors and their class labels in the sample.

ClassifierType::OutputType membershipSample =
classifier->GetOutput();
ClassifierType::OutputType::ConstIterator iter = membershipSample->Begin();

while (iter != membershipSample->End())
{
std::cout << "measurement vector = " << iter.GetMeasurementVector()
<< "class label = " << iter.GetClassLabel()
<< std::endl;
++iter;
}

#### 19.1.2 K-Means Classification

##### Simple version

The source code for this example can be found in the file
Examples/Classification/ScalarImageKmeansClassifier.cxx.

This example shows how to use the KMeans model for classifying the pixel of a scalar image.

The itk::Statistics::ScalarImageKmeansImageFilter is used for taking a scalar image and applying the K-Means algorithm in order to define classes that represents statistical distributions of intensity values in the pixels. The classes are then used in this filter for generating a labeled image where every pixel is assigned to one of the classes.

#include "otbImage.h"
#include "otbImageFileWriter.h"
#include "itkScalarImageKmeansImageFilter.h"

First we define the pixel type and dimension of the image that we intend to classify. With this image type we can also declare the otb::ImageFileReader needed for reading the input image, create one and set its input filename.

typedef signed short PixelType;
const unsigned int Dimension = 2;

typedef otb::Image<PixelType, Dimension> ImageType;

With the ImageType we instantiate the type of the itk::ScalarImageKmeansImageFilter that will compute the K-Means model and then classify the image pixels.

typedef itk::ScalarImageKmeansImageFilter<ImageType> KMeansFilterType;

KMeansFilterType::Pointer kmeansFilter = KMeansFilterType::New();

const unsigned int numberOfInitialClasses = atoi(argv[4]);

In general the classification will produce as output an image whose pixel values are integers associated to the labels of the classes. Since typically these integers will be generated in order (0, 1, 2, ...N), the output image will tend to look very dark when displayed with naive viewers. It is therefore convenient to have the option of spreading the label values over the dynamic range of the output image pixel type. When this is done, the dynamic range of the pixels is divided by the number of classes in order to define the increment between labels. For example, an output image of 8 bits will have a dynamic range of [0:255], and when it is used for holding four classes, the non-contiguous labels will be (0, 64, 128, 192). The selection of the mode to use is done with the method SetUseContiguousLabels().

const unsigned int useNonContiguousLabels = atoi(argv[3]);

kmeansFilter->SetUseNonContiguousLabels(useNonContiguousLabels);

For each one of the classes we must provide a tentative initial value for the mean of the class. Given that this is a scalar image, each one of the means is simply a scalar value. Note however that in a general case of K-Means, the input image would be a vector image and therefore the means will be vectors of the same dimension as the image pixels.

for (unsigned k = 0; k < numberOfInitialClasses; ++k)
{
const double userProvidedInitialMean = atof(argv[k + argoffset]);
}

The itk::ScalarImageKmeansImageFilter is predefined for producing an 8 bits scalar image as output. This output image contains labels associated to each one of the classes in the K-Means algorithm. In the following lines we use the OutputImageType in order to instantiate the type of a otb::ImageFileWriter. Then create one, and connect it to the output of the classification filter.

typedef KMeansFilterType::OutputImageType OutputImageType;

typedef otb::ImageFileWriter<OutputImageType> WriterType;

WriterType::Pointer writer = WriterType::New();

writer->SetInput(kmeansFilter->GetOutput());

writer->SetFileName(outputImageFileName);

We are now ready for triggering the execution of the pipeline. This is done by simply invoking the Update() method in the writer. This call will propagate the update request to the reader and then to the classifier.

try
{
writer->Update();
}
catch (itk::ExceptionObject& excp)
{
std::cerr << "Problem encountered while writing ";
std::cerr << " image file : " << argv[2] << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}

At this point the classification is done, the labeled image is saved in a file, and we can take a look at the means that were found as a result of the model estimation performed inside the classifier filter.

KMeansFilterType::ParametersType estimatedMeans =
kmeansFilter->GetFinalMeans();

const unsigned int numberOfClasses = estimatedMeans.Size();

for (unsigned int i = 0; i < numberOfClasses; ++i)
{
std::cout << "cluster[" << i << "] ";
std::cout << "    estimated mean : " << estimatedMeans[i] << std::endl;
}

Figure 19.4 illustrates the effect of this filter with three classes. The means can be estimated by ScalarImageKmeansModelEstimator.cxx.

The source code for this example can be found in the file
Examples/Classification/ScalarImageKmeansModelEstimator.cxx.

This example shows how to compute the KMeans model of an Scalar Image.

The itk::Statistics::KdTreeBasedKmeansEstimator is used for taking a scalar image and applying the K-Means algorithm in order to define classes that represents statistical distributions of intensity values in the pixels. One of the drawbacks of this technique is that the spatial distribution of the pixels is not considered at all. It is common therefore to combine the classification resulting from K-Means with other segmentation techniques that will use the classification as a prior and add spatial information to it in order to produce a better segmentation.

// Create a List from the scalar image

// Define the Measurement vector type from the AdaptorType

// Create the K-d tree structure
typedef itk::Statistics::WeightedCentroidKdTreeGenerator<
TreeGeneratorType;

TreeGeneratorType::Pointer treeGenerator = TreeGeneratorType::New();

treeGenerator->SetBucketSize(16);
treeGenerator->Update();

typedef TreeGeneratorType::KdTreeType                         TreeType;
typedef itk::Statistics::KdTreeBasedKmeansEstimator<TreeType> EstimatorType;

EstimatorType::Pointer estimator = EstimatorType::New();

const unsigned int numberOfClasses = 4;

EstimatorType::ParametersType initialMeans(numberOfClasses);
initialMeans[0] = 25.0;
initialMeans[1] = 125.0;
initialMeans[2] = 250.0;

estimator->SetParameters(initialMeans);

estimator->SetKdTree(treeGenerator->GetOutput());
estimator->SetMaximumIteration(200);
estimator->SetCentroidPositionChangesThreshold(0.0);
estimator->StartOptimization();

EstimatorType::ParametersType estimatedMeans = estimator->GetParameters();

for (unsigned int i = 0; i < numberOfClasses; ++i)
{
std::cout << "cluster[" << i << "] " << std::endl;
std::cout << "    estimated mean : " << estimatedMeans[i] << std::endl;
}
##### General approach

The source code for this example can be found in the file
Examples/Classification/KMeansImageClassificationExample.cxx.

The K-Means classification proposed by ITK for images is limited to scalar images and is not streamed. In this example, we show how the use of the otb::KMeansImageClassificationFilter allows for a simple implementation of a K-Means classification application. We will start by including the appropirate header file.

#include "otbKMeansImageClassificationFilter.h"

We will assume double precision input images and will also define the type for the labeled pixels.

const unsigned int Dimension = 2;
typedef double         PixelType;
typedef unsigned short LabeledPixelType;

Our classifier will be genric enough to be able to process images with any number of bands. We read the images as otb::VectorImages. The labeled image will be a scalar image.

typedef otb::VectorImage<PixelType, Dimension>  ImageType;
typedef otb::Image<LabeledPixelType, Dimension> LabeledImageType;

We can now define the type for the classifier filter, which is templated over its input and output image types.

typedef otb::KMeansImageClassificationFilter<ImageType, LabeledImageType>
ClassificationFilterType;
typedef ClassificationFilterType::KMeansParametersType KMeansParametersType;

And finally, we define the reader and the writer. Since the images to classify can be very big, we will use a streamed writer which will trigger the streaming ability of the classifier.

typedef otb::ImageFileWriter<LabeledImageType> WriterType;

We instantiate the classifier and the reader objects and we set their parameters. Please note the call of the GenerateOutputInformation() method on the reader in order to have available the information about the input image (size, number of bands, etc.) without needing to actually read the image.

ClassificationFilterType::Pointer filter = ClassificationFilterType::New();

The classifier needs as input the centroids of the classes. We declare the parameter vector, and we read the centroids from the arguments of the program.

const unsigned int sampleSize =
ClassificationFilterType::MaxSampleDimension;
const unsigned int   parameterSize = nbClasses  sampleSize;
KMeansParametersType parameters;

parameters.SetSize(parameterSize);
parameters.Fill(0);

for (unsigned int i = 0; i < nbClasses; ++i)
{
for (unsigned int j = 0; j <
{
parameters[i  sampleSize + j] =
atof(argv[4 + i
+ j]);
}
}

std::cout << "Parameters: " << parameters << std::endl;

We set the parameters for the classifier, we plug the pipeline and trigger its execution by updating the output of the writer.

filter->SetCentroids(parameters);

WriterType::Pointer writer = WriterType::New();
writer->SetInput(filter->GetOutput());
writer->SetFileName(outfname);
writer->Update();

#### 19.1.3 Bayesian Plug-In Classifier

The source code for this example can be found in the file
Examples/Classification/BayesianPluginClassifier.cxx.

In this example, we present a system that places measurement vectors into two Gaussian classes. The Figure 19.5 shows all the components of the classifier system and the data flow. This system differs with the previous k-means clustering algorithms in several ways. The biggest difference is that this classifier uses the itk::Statistics::GaussianDensityFunctions as membership functions instead of the itk::Statistics::EuclideanDistance. Since the membership function is different, the membership function requires a different set of parameters, mean vectors and covariance matrices. We choose the itk::Statistics::MeanCalculator (sample mean) and the itk::Statistics::CovarianceCalculator (sample covariance) for the estimation algorithms of the two parameters. If we want more robust estimation algorithm, we can replace these estimation algorithms with more alternatives without changing other components in the classifier system.

It is a bad idea to use the same sample for test and training (parameter estimation) of the parameters. However, for simplicity, in this example, we use a sample for test and training.

We use the itk::Statistics::ListSample as the sample (test and training). The itk::Vector is our measurement vector class. To store measurement vectors into two separate sample containers, we use the itk::Statistics::Subsample objects.

#include "itkVector.h"
#include "itkListSample.h"
#include "itkSubsample.h"

The following two files provides us the parameter estimation algorithms.

#include "itkMeanCalculator.h"
#include "itkCovarianceCalculator.h"

The following files define the components required by ITK statistical classification framework: the decision rule, the membership function, and the classifier.

#include "itkMaximumRatioDecisionRule.h"
#include "itkGaussianDensityFunction.h"
#include "itkSampleClassifier.h"

We will fill the sample with random variables from two normal distribution using the itk::Statistics::NormalVariateGenerator.

#include "itkNormalVariateGenerator.h"

Since the NormalVariateGenerator class only supports 1-D, we define our measurement vector type as a one component vector. We then, create a ListSample object for data inputs.

We also create two Subsample objects that will store the measurement vectors in sample into two separate sample containers. Each Subsample object stores only the measurement vectors belonging to a single class. This class sample will be used by the parameter estimation algorithms.

typedef itk::Vector<double,
1>
MeasurementVectorType;
typedef itk::Statistics::ListSample<MeasurementVectorType> SampleType;
SampleType::Pointer sample = SampleType::New();
sample->SetMeasurementVectorSize(1);   // length of measurement vectors
// in the sample.

typedef itk::Statistics::Subsample<SampleType> ClassSampleType;
std::vector<ClassSampleType::Pointer> classSamples;
for (unsigned int i = 0; i < 2; ++i)
{
classSamples.push_back(ClassSampleType::New());
classSamples[i]->SetSample(sample);
}

The following code snippet creates a NormalVariateGenerator object. Since the random variable generator returns values according to the standard normal distribution (the mean is zero, and the standard deviation is one) before pushing random values into the sample, we change the mean and standard deviation. We want two normal (Gaussian) distribution data. We have two for loops. Each for loop uses different mean and standard deviation. Before we fill the sample with the second distribution data, we call Initialize(random seed) method, to recreate the pool of random variables in the normalGenerator. In the second for loop, we fill the two class samples with measurement vectors using the AddInstance() method.

To see the probability density plots from the two distributions, refer to Figure 19.3.

typedef itk::Statistics::NormalVariateGenerator NormalGeneratorType;
NormalGeneratorType::Pointer normalGenerator = NormalGeneratorType::New();

normalGenerator->Initialize(101);

MeasurementVectorType          mv;
double                         mean = 100;
double                         standardDeviation = 30;
SampleType::InstanceIdentifier id = 0UL;
for (unsigned int i = 0; i < 100; ++i)
{
mv.Fill((normalGenerator->GetVariate()  standardDeviation) + mean);
sample->PushBack(mv);
++id;
}

normalGenerator->Initialize(3024);
mean = 200;
standardDeviation = 30;
for (unsigned int i = 0; i < 100; ++i)
{
mv.Fill((normalGenerator->GetVariate()  standardDeviation) + mean);
sample->PushBack(mv);
++id;
}

In the following code snippet, notice that the template argument for the MeanCalculator and CovarianceCalculator is ClassSampleType (i.e., type of Subsample) instead of SampleType (i.e., type of ListSample). This is because the parameter estimation algorithms are applied to the class sample.

typedef itk::Statistics::MeanCalculator<ClassSampleType> MeanEstimatorType;
typedef itk::Statistics::CovarianceCalculator<ClassSampleType>
CovarianceEstimatorType;

std::vector<MeanEstimatorType::Pointer>       meanEstimators;
std::vector<CovarianceEstimatorType::Pointer> covarianceEstimators;

for (unsigned int i = 0; i < 2; ++i)
{
meanEstimators.push_back(MeanEstimatorType::New());
meanEstimators[i]->SetInputSample(classSamples[i]);
meanEstimators[i]->Update();

covarianceEstimators.push_back(CovarianceEstimatorType::New());
covarianceEstimators[i]->SetInputSample(classSamples[i]);
covarianceEstimators[i]->SetMean(meanEstimators[i]->GetOutput());
covarianceEstimators[i]->Update();
}

We print out the estimated parameters.

for (unsigned int i = 0; i < 2; ++i)
{
std::cout << "class[" << i << "] " << std::endl;
std::cout << "    estimated mean : "
<< ⋆(meanEstimators[i]->GetOutput())
<< "    covariance matrix : "
<< ⋆(covarianceEstimators[i]->GetOutput()) << std::endl;
}

After creating a SampleClassifier object and a MaximumRatioDecisionRule object, we plug in the decisionRule and the sample to the classifier. Then, we specify the number of classes that will be considered using the SetNumberOfClasses() method.

The MaximumRatioDecisionRule requires a vector of a priori probability values. Such a priori probability will be the P(ωi) of the following variation of the Bayes decision rule:

 (19.1)

The remainder of the code snippet shows how to use user-specified class labels. The classification result will be stored in a MembershipSample object, and for each measurement vector, its class label will be one of the two class labels, 100 and 200 (unsigned int).

typedef itk::Statistics::GaussianDensityFunction
<MeasurementVectorType> MembershipFunctionType;
typedef itk::MaximumRatioDecisionRule DecisionRuleType;
DecisionRuleType::Pointer decisionRule = DecisionRuleType::New();

DecisionRuleType::APrioriVectorType aPrioris;
aPrioris.push_back(classSamples[0]->GetTotalFrequency()
/ sample->GetTotalFrequency());
aPrioris.push_back(classSamples[1]->GetTotalFrequency()
/ sample->GetTotalFrequency());
decisionRule->SetAPriori(aPrioris);

typedef itk::Statistics::SampleClassifier<SampleType> ClassifierType;
ClassifierType::Pointer classifier = ClassifierType::New();

classifier->SetDecisionRule((itk::DecisionRuleBase::Pointer) decisionRule);
classifier->SetSample(sample);
classifier->SetNumberOfClasses(2);

std::vector<unsigned int> classLabels;
classLabels.resize(2);
classLabels[0] = 100;
classLabels[1] = 200;
classifier->SetMembershipFunctionClassLabels(classLabels);

The classifier is almost ready to perform the classification except that it needs two membership functions that represent the two clusters.

In this example, we can imagine that the two clusters are modeled by two Euclidean distance functions. The distance function (model) has only one parameter, the mean (centroid) set by the SetOrigin() method. To plug-in two distance functions, we call the AddMembershipFunction() method. Then invocation of the Update() method will perform the classification.

std::vector<MembershipFunctionType::Pointer> membershipFunctions;
for (unsigned int i = 0; i < 2; ++i)
{
membershipFunctions.push_back(MembershipFunctionType::New());
membershipFunctions[i]->SetMean(meanEstimators[i]->GetOutput());
membershipFunctions[i]->
SetCovariance(covarianceEstimators[i]->GetOutput());
}

classifier->Update();

The following code snippet prints out pairs of a measurement vector and its class label in the sample.

ClassifierType::OutputType membershipSample =
classifier->GetOutput();
ClassifierType::OutputType::ConstIterator iter = membershipSample->Begin();

while (iter != membershipSample->End())
{
std::cout << "measurement vector = " << iter.GetMeasurementVector()
<< "class label = " << iter.GetClassLabel() << std::endl;
++iter;
}

#### 19.1.4 Expectation Maximization Mixture Model Estimation

The source code for this example can be found in the file
Examples/Classification/ExpectationMaximizationMixtureModelEstimator.cxx.

In this example, we present ITK’s implementation of the expectation maximization (EM) process to generate parameter estimates for a two Gaussian component mixture model.

The Bayesian plug-in classifier example (see Section 19.1.3) used two Gaussian probability density functions (PDF) to model two Gaussian distribution classes (two models for two class). However, in some cases, we want to model a distribution as a mixture of several different distributions. Therefore, the probability density function (p(x)) of a mixture model can be stated as follows :

 (19.2)

where i is the index of the component, c is the number of components, αi is the proportion of the component, and fi is the probability density function of the component.

Now the task is to find the parameters(the component PDF’s parameters and the proportion values) to maximize the likelihood of the parameters. If we know which component a measurement vector belongs to, the solutions to this problem is easy to solve. However, we don’t know the membership of each measurement vector. Therefore, we use the expectation of membership instead of the exact membership. The EM process splits into two steps:

1. E step: calculate the expected membership values for each measurement vector to each classes.
2. M step: find the next parameter sets that maximize the likelihood with the expected membership values and the current set of parameters.

The E step is basically a step that calculates the a posteriori probability for each measurement vector.

The M step is dependent on the type of each PDF. Most of distributions belonging to exponential family such as Poisson, Binomial, Exponential, and Normal distributions have analytical solutions for updating the parameter set. The itk::Statistics::ExpectationMaximizationMixtureModelEstimator class assumes that such type of components.

In the following example we use the itk::Statistics::ListSample as the sample (test and training). The itk::Vector::is our measurement vector class. To store measurement vectors into two separate sample container, we use the itk::Statistics::Subsample objects.

#include "itkVector.h"
#include "itkListSample.h"

The following two files provide us the parameter estimation algorithms.

#include "itkGaussianMixtureModelComponent.h"
#include "itkExpectationMaximizationMixtureModelEstimator.h"

We will fill the sample with random variables from two normal distribution using the itk::Statistics::NormalVariateGenerator.

#include "itkNormalVariateGenerator.h"

Since the NormalVariateGenerator class only supports 1-D, we define our measurement vector type as a one component vector. We then, create a ListSample object for data inputs.

We also create two Subsample objects that will store the measurement vectors in the sample into two separate sample containers. Each Subsample object stores only the measurement vectors belonging to a single class. This class sample will be used by the parameter estimation algorithms.

unsigned int numberOfClasses = 2;
typedef itk::Vector<double,
1>
MeasurementVectorType;
typedef itk::Statistics::ListSample<MeasurementVectorType> SampleType;
SampleType::Pointer sample = SampleType::New();
sample->SetMeasurementVectorSize(1);   // length of measurement vectors
// in the sample.

The following code snippet creates a NormalVariateGenerator object. Since the random variable generator returns values according to the standard normal distribution (the mean is zero, and the standard deviation is one) before pushing random values into the sample, we change the mean and standard deviation. We want two normal (Gaussian) distribution data. We have two for loops. Each for loop uses different mean and standard deviation. Before we fill the sample with the second distribution data, we call Initialize() method to recreate the pool of random variables in the normalGenerator. In the second for loop, we fill the two class samples with measurement vectors using the AddInstance() method.

To see the probability density plots from the two distribution, refer to Figure 19.3.

typedef itk::Statistics::NormalVariateGenerator NormalGeneratorType;
NormalGeneratorType::Pointer normalGenerator = NormalGeneratorType::New();

normalGenerator->Initialize(101);

MeasurementVectorType mv;
double                mean = 100;
double                standardDeviation = 30;
for (unsigned int i = 0; i < 100; ++i)
{
mv[0] = (normalGenerator->GetVariate()  standardDeviation) + mean;
sample->PushBack(mv);
}

normalGenerator->Initialize(3024);
mean = 200;
standardDeviation = 30;
for (unsigned int i = 0; i < 100; ++i)
{
mv[0] = (normalGenerator->GetVariate()  standardDeviation) + mean;
sample->PushBack(mv);
}

In the following code snippet notice that the template argument for the MeanCalculator and CovarianceCalculator is ClassSampleType (i.e., type of Subsample) instead of SampleType (i.e., type of ListSample). This is because the parameter estimation algorithms are applied to the class sample.

typedef itk::Array<double> ParametersType;
ParametersType params(2);

std::vector<ParametersType> initialParameters(numberOfClasses);
params[0] = 110.0;
params[1] = 800.0;
initialParameters[0] = params;

params[0] = 210.0;
params[1] = 850.0;
initialParameters[1] = params;

typedef itk::Statistics::GaussianMixtureModelComponent<SampleType>
ComponentType;

std::vector<ComponentType::Pointer> components;
for (unsigned int i = 0; i < numberOfClasses; ++i)
{
components.push_back(ComponentType::New());
(components[i])->SetSample(sample);
(components[i])->SetParameters(initialParameters[i]);
}

We run the estimator.

typedef itk::Statistics::ExpectationMaximizationMixtureModelEstimator<
SampleType> EstimatorType;
EstimatorType::Pointer estimator = EstimatorType::New();

estimator->SetSample(sample);
estimator->SetMaximumIteration(200);

itk::Array<double> initialProportions(numberOfClasses);
initialProportions[0] = 0.5;
initialProportions[1] = 0.5;

estimator->SetInitialProportions(initialProportions);

for (unsigned int i = 0; i < numberOfClasses; ++i)
{
(components[i]).GetPointer());
}

estimator->Update();

We then print out the estimated parameters.

for (unsigned int i = 0; i < numberOfClasses; ++i)
{
std::cout << "Cluster[" << i << "]" << std::endl;
std::cout << "    Parameters:" << std::endl;
std::cout << "         " << (components[i])->GetFullParameters()
<< std::endl;
std::cout << "    Proportion: ";
std::cout << "         " << (⋆estimator->GetProportions())[i] << std::endl;
}

#### 19.1.5 Classification using Markov Random Fields

Markov Random Fields are probabilistic models that use the statistical dependency between pixels in a neighborhood to infeer the value of a give pixel.

##### ITK framework

The itk::Statistics::MRFImageFilter uses the maximum a posteriori (MAP) estimates for modeling the MRF. The object traverses the data set and uses the model generated by the Mahalanobis distance classifier to get the the distance between each pixel in the data set to a set of known classes, updates the distances by evaluating the influence of its neighboring pixels (based on a MRF model) and finally, classifies each pixel to the class which has the minimum distance to that pixel (taking the neighborhood influence under consideration). The energy function minimization is done using the iterated conditional modes (ICM) algorithm [13].

The source code for this example can be found in the file
Examples/Classification/ScalarImageMarkovRandomField1.cxx.

This example shows how to use the Markov Random Field approach for classifying the pixel of a scalar image.

The itk::Statistics::MRFImageFilter is used for refining an initial classification by introducing the spatial coherence of the labels. The user should provide two images as input. The first image is the one to be classified while the second image is an image of labels representing an initial classification.

The following headers are related to reading input images, writing the output image, and making the necessary conversions between scalar and vector images.

#include "otbImage.h"
#include "itkFixedArray.h"
#include "otbImageFileWriter.h"
#include "itkScalarToArrayCastImageFilter.h"

The following headers are related to the statistical classification classes.

#include "itkMRFImageFilter.h"
#include "itkDistanceToCentroidMembershipFunction.h"
#include "itkMinimumDecisionRule.h"
#include "itkImageClassifierBase.h"

First we define the pixel type and dimension of the image that we intend to classify. With this image type we can also declare the otb::ImageFileReader needed for reading the input image, create one and set its input filename.

typedef unsigned char PixelType;
const unsigned int Dimension = 2;

typedef otb::Image<PixelType, Dimension> ImageType;

As a second step we define the pixel type and dimension of the image of labels that provides the initial classification of the pixels from the first image. This initial labeled image can be the output of a K-Means method like the one illustrated in section 19.1.2.

typedef unsigned char LabelPixelType;

typedef otb::Image<LabelPixelType, Dimension> LabelImageType;

Since the Markov Random Field algorithm is defined in general for images whose pixels have multiple components, that is, images of vector type, we must adapt our scalar image in order to satisfy the interface expected by the MRFImageFilter. We do this by using the itk::ScalarToArrayCastImageFilter. With this filter we will present our scalar image as a vector image whose vector pixels contain a single component.

typedef itk::FixedArray<LabelPixelType, 1> ArrayPixelType;

typedef otb::Image<ArrayPixelType, Dimension> ArrayImageType;

typedef itk::ScalarToArrayCastImageFilter<
ImageType, ArrayImageType> ScalarToArrayFilterType;

ScalarToArrayFilterType::Pointer
scalarToArrayFilter = ScalarToArrayFilterType::New();

With the input image type ImageType and labeled image type LabelImageType we instantiate the type of the itk::MRFImageFilter that will apply the Markov Random Field algorithm in order to refine the pixel classification.

typedef itk::MRFImageFilter<ArrayImageType, LabelImageType> MRFFilterType;

MRFFilterType::Pointer mrfFilter = MRFFilterType::New();

mrfFilter->SetInput(scalarToArrayFilter->GetOutput());

We set now some of the parameters for the MRF filter. In particular, the number of classes to be used during the classification, the maximum number of iterations to be run in this filter and the error tolerance that will be used as a criterion for convergence.

mrfFilter->SetNumberOfClasses(numberOfClasses);
mrfFilter->SetMaximumNumberOfIterations(numberOfIterations);
mrfFilter->SetErrorTolerance(1e-7);

The smoothing factor represents the tradeoff between fidelity to the observed image and the smoothness of the segmented image. Typical smoothing factors have values between 1 5. This factor will multiply the weights that define the influence of neighbors on the classification of a given pixel. The higher the value, the more uniform will be the regions resulting from the classification refinement.

mrfFilter->SetSmoothingFactor(smoothingFactor);

Given that the MRF filter needs to continually relabel the pixels, it needs access to a set of membership functions that will measure to what degree every pixel belongs to a particular class. The classification is performed by the itk::ImageClassifierBase class, that is instantiated using the type of the input vector image and the type of the labeled image.

typedef itk::ImageClassifierBase<
ArrayImageType,
LabelImageType>   SupervisedClassifierType;

SupervisedClassifierType::Pointer classifier =
SupervisedClassifierType::New();

The classifier needs a decision rule to be set by the user. Note that we must use GetPointer() in the call of the SetDecisionRule() method because we are passing a SmartPointer, and smart pointer cannot perform polymorphism, we must then extract the raw pointer that is associated to the smart pointer. This extraction is done with the GetPointer() method.

typedef itk::MinimumDecisionRule DecisionRuleType;

DecisionRuleType::Pointer classifierDecisionRule = DecisionRuleType::New();

classifier->SetDecisionRule(classifierDecisionRule.GetPointer());

We now instantiate the membership functions. In this case we use the itk::Statistics::DistanceToCentroidMembershipFunction class templated over the pixel type of the vector image, which in our example happens to be a vector of dimension 1.

typedef itk::Statistics::DistanceToCentroidMembershipFunction<
ArrayPixelType>
MembershipFunctionType;

typedef MembershipFunctionType::Pointer MembershipFunctionPointer;

double             meanDistance = 0;
vnl_vector<double> centroid(1);
for (unsigned int i = 0; i < numberOfClasses; ++i)
{
MembershipFunctionPointer membershipFunction =
MembershipFunctionType::New();

centroid[0] = atof(argv[i + numberOfArgumentsBeforeMeans]);

membershipFunction->SetCentroid(centroid);

meanDistance += static_cast<double> (centroid[0]);
}
meanDistance /= numberOfClasses;

and we set the neighborhood radius that will define the size of the clique to be used in the computation of the neighbors’ influence in the classification of any given pixel. Note that despite the fact that we call this a radius, it is actually the half size of an hypercube. That is, the actual region of influence will not be circular but rather an N-Dimensional box. For example, a neighborhood radius of 2 in a 3D image will result in a clique of size 5x5x5 pixels, and a radius of 1 will result in a clique of size 3x3x3 pixels.

We should now set the weights used for the neighbors. This is done by passing an array of values that contains the linear sequence of weights for the neighbors. For example, in a neighborhood of size 3x3x3, we should provide a linear array of 9 weight values. The values are packaged in a std::vector and are supposed to be double. The following lines illustrate a typical set of values for a 3x3x3 neighborhood. The array is arranged and then passed to the filter by using the method SetMRFNeighborhoodWeight().

std::vector<double> weights;
weights.push_back(1.5);
weights.push_back(2.0);
weights.push_back(1.5);
weights.push_back(2.0);
weights.push_back(0.0); // This is the central pixel
weights.push_back(2.0);
weights.push_back(1.5);
weights.push_back(2.0);
weights.push_back(1.5);

We now scale weights so that the smoothing function and the image fidelity functions have comparable value. This is necessary since the label image and the input image can have different dynamic ranges. The fidelity function is usually computed using a distance function, such as the itk::DistanceToCentroidMembershipFunction or one of the other membership functions. They tend to have values in the order of the means specified.

double totalWeight = 0;
for (std::vector<double>::const_iterator wcIt = weights.begin();
wcIt != weights.end(); ++wcIt)
{
totalWeight += wcIt;
}
for (std::vector<double>::iterator wIt = weights.begin();
wIt != weights.end(); wIt++)
{
wIt = static_cast<double> ((⋆wIt)  meanDistance / (2  totalWeight));
}

mrfFilter->SetMRFNeighborhoodWeight(weights);

Finally, the classifier class is connected to the Markov Random Fields filter.

mrfFilter->SetClassifier(classifier);

The output image produced by the itk::MRFImageFilter has the same pixel type as the labeled input image. In the following lines we use the OutputImageType in order to instantiate the type of a otb::ImageFileWriter. Then create one, and connect it to the output of the classification filter after passing it through an intensity rescaler to rescale it to an 8 bit dynamic range

typedef MRFFilterType::OutputImageType OutputImageType;
typedef otb::ImageFileWriter<OutputImageType> WriterType;

WriterType::Pointer writer = WriterType::New();

writer->SetInput(intensityRescaler->GetOutput());

writer->SetFileName(outputImageFileName);

We are now ready for triggering the execution of the pipeline. This is done by simply invoking the Update() method in the writer. This call will propagate the update request to the reader and then to the MRF filter.

try
{
writer->Update();
}
catch (itk::ExceptionObject& excp)
{
std::cerr << "Problem encountered while writing ";
std::cerr << " image file : " << argv[2] << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}

Figure 19.6 illustrates the effect of this filter with four classes. In this example the filter was run with a smoothing factor of 3. The labeled image was produced by ScalarImageKmeansClassifier.cxx and the means were estimated by ScalarImageKmeansModelEstimator.cxx described in section 19.1.2. The obtained result can be compared with the one of figure 19.4 to see the interest of using the MRF approach in order to ensure the regularization of the classified image.

##### OTB framework

The ITK approach was considered not to be flexible enough for some remote sensing applications. Therefore, we decided to implement our own framework.

The source code for this example can be found in the file
Examples/Markov/MarkovClassification1Example.cxx.

This example illustrates the details of the otb::MarkovRandomFieldFilter. This filter is an application of the Markov Random Fields for classification, segmentation or restauration.

This example applies the otb::MarkovRandomFieldFilter to classify an image into four classes defined by their mean and variance. The optimization is done using an Metropolis algorithm with a random sampler. The regularization energy is defined by a Potts model and the fidelity by a Gaussian model.

The first step toward the use of this filter is the inclusion of the proper header files.

#include "otbMRFEnergyPotts.h"
#include "otbMRFEnergyGaussianClassification.h"
#include "otbMRFOptimizerMetropolis.h"
#include "otbMRFSamplerRandom.h"

Then we must decide what pixel type to use for the image. We choose to make all computations with double precision. The labelled image is of type unsigned char which allows up to 256 different classes.

const unsigned int Dimension = 2;

typedef double                                   InternalPixelType;
typedef unsigned char                            LabelledPixelType;
typedef otb::Image<InternalPixelType, Dimension> InputImageType;
typedef otb::Image<LabelledPixelType, Dimension> LabelledImageType;

We define a reader for the image to be classified, an initialization for the classification (which could be random) and a writer for the final classification.

typedef otb::ImageFileWriter<LabelledImageType> WriterType;

WriterType::Pointer writer = WriterType::New();

const char  inputFilename  = argv[1];
const char  outputFilename = argv[2];

writer->SetFileName(outputFilename);

Finally, we define the different classes necessary for the Markov classification. A otb::MarkovRandomFieldFilter is instanciated, this is the main class which connect the other to do the Markov classification.

typedef otb::MarkovRandomFieldFilter
<InputImageType, LabelledImageType> MarkovRandomFieldFilterType;

An otb::MRFSamplerRandomMAP, which derives from the otb::MRFSampler, is instanciated. The sampler is in charge of proposing a modification for a given site. The otb::MRFSamplerRandomMAP, randomly pick one possible value according to the MAP probability.

typedef otb::MRFSamplerRandom<InputImageType, LabelledImageType> SamplerType;

An otb::MRFOptimizerMetropoli, which derives from the otb::MRFOptimizer, is instanciated. The optimizer is in charge of accepting or rejecting the value proposed by the sampler. The otb::MRFSamplerRandomMAP, accept the proposal according to the variation of energy it causes and a temperature parameter.

typedef otb::MRFOptimizerMetropolis OptimizerType;

Two energy, deriving from the otb::MRFEnergy class need to be instanciated. One energy is required for the regularization, taking into account the relashionship between neighborhing pixels in the classified image. Here it is done with the otb::MRFEnergyPotts which implement a Potts model.

The second energy is for the fidelity to the original data. Here it is done with an otb::MRFEnergyGaussianClassification class, which defines a gaussian model for the data.

typedef otb::MRFEnergyPotts
<LabelledImageType, LabelledImageType>  EnergyRegularizationType;
typedef otb::MRFEnergyGaussianClassification
<InputImageType, LabelledImageType>  EnergyFidelityType;

The different filters composing our pipeline are created by invoking their New() methods, assigning the results to smart pointers.

MarkovRandomFieldFilterType::Pointer markovFilter =
MarkovRandomFieldFilterType::New();
EnergyRegularizationType::Pointer energyRegularization =
EnergyRegularizationType::New();
EnergyFidelityType::Pointer energyFidelity = EnergyFidelityType::New();
OptimizerType::Pointer      optimizer = OptimizerType::New();
SamplerType::Pointer        sampler = SamplerType::New();

Parameter for the otb::MRFEnergyGaussianClassification class, meand and standard deviation are created.

unsigned int nClass = 4;
energyFidelity->SetNumberOfParameters(2  nClass);
EnergyFidelityType::ParametersType parameters;
parameters.SetSize(energyFidelity->GetNumberOfParameters());
parameters[0] = 10.0; //Class 0 mean
parameters[1] = 10.0; //Class 0 stdev
parameters[2] = 80.0; //Class 1 mean
parameters[3] = 10.0; //Class 1 stdev
parameters[4] = 150.0; //Class 2 mean
parameters[5] = 10.0; //Class 2 stdev
parameters[6] = 220.0; //Class 3 mean
parameters[7] = 10.0; //Class 3 stde
energyFidelity->SetParameters(parameters);

Parameters are given to the different class an the sampler, optimizer and energies are connected with the Markov filter.

OptimizerType::ParametersType param(1);
param.Fill(atof(argv[5]));
optimizer->SetParameters(param);
markovFilter->SetNumberOfClasses(nClass);
markovFilter->SetMaximumNumberOfIterations(atoi(argv[4]));
markovFilter->SetErrorTolerance(0.0);
markovFilter->SetLambda(atof(argv[3]));

markovFilter->SetEnergyRegularization(energyRegularization);
markovFilter->SetEnergyFidelity(energyFidelity);
markovFilter->SetOptimizer(optimizer);
markovFilter->SetSampler(sampler);

The pipeline is connected. An itk::RescaleIntensityImageFilter rescale the classified image before saving it.

typedef itk::RescaleIntensityImageFilter
<LabelledImageType, LabelledImageType> RescaleType;
RescaleType::Pointer rescaleFilter = RescaleType::New();
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(255);

rescaleFilter->SetInput(markovFilter->GetOutput());

writer->SetInput(rescaleFilter->GetOutput());

Finally, the pipeline execution is trigerred.

writer->Update();

Figure 19.8 shows the output of the Markov Random Field classification after 20 iterations with a random sampler and a Metropolis optimizer.

The source code for this example can be found in the file
Examples/Markov/MarkovClassification2Example.cxx.

Using a similar structure as the previous program and the same energy function, we are now going to slightly alter the program to use a different sampler and optimizer. The proposed sample is proposed randomly according to the MAP probability and the optimizer is the ICM which accept the proposed sample if it enable a reduction of the energy.

First, we need to include header specific to these class:

#include "otbMRFSamplerRandomMAP.h"
#include "otbMRFOptimizerICM.h"

And to declare these new type:

typedef otb::MRFSamplerRandomMAP<InputImageType,
LabelledImageType> SamplerType;
//   typedef otb::MRFSamplerRandom< InputImageType, LabelledImageType> SamplerType;
typedef otb::MRFOptimizerICM OptimizerType;

As the otb::MRFOptimizerICM does not have any parameters, the call to optimizer->SetParameters() must be removed

Apart from these, no further modification is required.

Figure 19.9 shows the output of the Markov Random Field classification after 5 iterations with a MAP random sampler and an ICM optimizer.

The source code for this example can be found in the file
Examples/Markov/MarkovClassification3Example.cxx.

This example illustrates the details of the MarkovRandomFieldFilter by using the Fisher distribution to model the likelihood energy. This filter is an application of the Markov Random Fields for classification.

This example applies the MarkovRandomFieldFilter to classify an image into four classes defined by their Fisher distribution parameters L, M and mu. The optimization is done using a Metropolis algorithm with a random sampler. The regularization energy is defined by a Potts model and the fidelity or likelihood energy is modelled by a Fisher distribution. The parameter of the Fisher distribution was determined for each class in a supervised step. ( See the File OtbParameterEstimatioOfFisherDistribution ) This example is a contribution from Jan Wegner.

Then we must decide what pixel type to use for the image. We choose to make all computations with double precision. The labeled image is of type unsigned char which allows up to 256 different classes.

const unsigned int Dimension = 2;
typedef double        InternalPixelType;
typedef unsigned char LabelledPixelType;

typedef otb::Image<InternalPixelType, Dimension>    InputImageType;
typedef otb::Image<LabelledPixelType, Dimension>    LabelledImageType;

We define a reader for the image to be classified, an initialization for the classification (which could be random) and a writer for the final classification.

typedef otb::ImageFileWriter< LabelledImageType >  WriterType;

WriterType::Pointer writer = WriterType::New();

Finally, we define the different classes necessary for the Markov classification. A MarkovRandomFieldFilter is instantiated, this is the main class which connect the other to do the Markov classification.

typedef otb::MarkovRandomFieldFilter
<InputImageType, LabelledImageType> MarkovRandomFieldFilterType;

An MRFSamplerRandomMAP, which derives from the MRFSampler, is instanciated. The sampler is in charge of proposing a modification for a given site. The MRFSamplerRandomMAP, randomly pick one possible value according to the MAP probability.

typedef otb::MRFSamplerRandom< InputImageType, LabelledImageType> SamplerType;

An MRFOptimizerMetropolis, which derives from the MRFOptimizer, is instanciated. The optimizer is in charge of accepting or rejecting the value proposed by the sampler. The MRFSamplerRandomMAP, accept the proposal according to the variation of energy it causes and a temperature parameter.

typedef otb::MRFOptimizerMetropolis OptimizerType;

Two energy, deriving from the MRFEnergy class need to be instantiated. One energy is required for the regularization, taking into account the relationship between neighboring pixels in the classified image. Here it is done with the MRFEnergyPotts, which implements a Potts model.

The second energy is used for the fidelity to the original data. Here it is done with a MRFEnergyFisherClassification class, which defines a Fisher distribution to model the data.

typedef otb::MRFEnergyPotts
<LabelledImageType, LabelledImageType>  EnergyRegularizationType;
typedef otb::MRFEnergyFisherClassification
<InputImageType, LabelledImageType>     EnergyFidelityType;

The different filters composing our pipeline are created by invoking their New() methods, assigning the results to smart pointers.

MarkovRandomFieldFilterType::Pointer markovFilter = MarkovRandomFieldFilterType::New();
EnergyRegularizationType::Pointer energyRegularization = EnergyRegularizationType::New();
EnergyFidelityType::Pointer energyFidelity = EnergyFidelityType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
SamplerType::Pointer sampler = SamplerType::New();

Parameter for the MRFEnergyFisherClassification class are created. The shape parameters M, L and the weighting parameter mu are computed in a supervised step

unsigned int nClass =4;
energyFidelity->SetNumberOfParameters(3⋆nClass);
EnergyFidelityType::ParametersType parameters;
parameters.SetSize(energyFidelity->GetNumberOfParameters());
//Class 0
parameters[0] =         12.353042;    //Class 0 mu
parameters[1] =         2.156422;       //Class 0 L
parameters[2] =         4.920403;       //Class 0 M
//Class 1
parameters[3] =         72.068291;    //Class 1 mu
parameters[4] =         11.000000;   //Class 1 L
parameters[5] =         50.950001;    //Class 1 M
//Class 2
parameters[6] =         146.665985;   //Class 2 mu
parameters[7] =         11.000000;   //Class 2 L
parameters[8] =         50.900002;    //Class 2 M
//Class 3
parameters[9]  =      200.010132;     //Class 3 mu
parameters[10] =      11.000000;      //Class 3 L
parameters[11] =      50.950001;      //Class 3 M

energyFidelity->SetParameters(parameters);

Parameters are given to the different classes and the sampler, optimizer and energies are connected with the Markov filter.

OptimizerType::ParametersType param(1);
param.Fill(atof(argv[6]));
optimizer->SetParameters(param);
markovFilter->SetNumberOfClasses(nClass);
markovFilter->SetMaximumNumberOfIterations(atoi(argv[5]));
markovFilter->SetErrorTolerance(0.0);
markovFilter->SetLambda(atof(argv[4]));

markovFilter->SetEnergyRegularization(energyRegularization);
markovFilter->SetEnergyFidelity(energyFidelity);
markovFilter->SetOptimizer(optimizer);
markovFilter->SetSampler(sampler);

The pipeline is connected. An itkRescaleIntensityImageFilter rescales the classified image before saving it.

typedef itk::RescaleIntensityImageFilter
< LabelledImageType, LabelledImageType > RescaleType;
RescaleType::Pointer rescaleFilter = RescaleType::New();
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(255);

rescaleFilter->SetInput( markovFilter->GetOutput() );

writer->SetInput( rescaleFilter->GetOutput() );
writer->Update();

We can now create an image file writer and save the image.

typedef otb::ImageFileWriter<RGBImageType> WriterRescaledType;

WriterRescaledType::Pointer writerRescaled = WriterRescaledType::New();

writerRescaled->SetFileName( outputRescaledImageFileName );
writerRescaled->SetInput( colormapper->GetOutput() );

writerRescaled->Update();

Figure 19.10 shows the output of the Markov Random Field classification into four classes using the Fisher-distribution as likelihood term.

The source code for this example can be found in the file
Examples/Markov/MarkovRegularizationExample.cxx.

This example illustrates the use of the otb::MarkovRandomFieldFilter. to regularize a classification obtained previously by another classifier. Here we will apply the regularization to the output of an SVM classifier presented in 19.3.5.

The reference image and the starting image are both going to be the original classification. Both regularization and fidelity energy are defined by Potts model.

The convergence of the Markov Random Field is done with a random sampler and a Metropolis model as in example 1. As you should get use to the general program structure to use the MRF framework, we are not going to repeat the entire example. However, remember you can find the full source code for this example in your OTB source directory.

To find the number of classes available in the original image we use the itk::LabelStatisticsImageFilter and more particularly the method GetNumberOfLabels().

typedef itk::LabelStatisticsImageFilter
<LabelledImageType, LabelledImageType> LabelledStatType;
LabelledStatType::Pointer labelledStat = LabelledStatType::New();
labelledStat->Update();

unsigned int nClass = labelledStat->GetNumberOfLabels();

Figure 19.11 shows the output of the Markov Random Field regularization on the classification output of another method.

### 19.2 Statistical Segmentations

#### 19.2.1 Stochastic Expectation Maximization

The Stochastic Expectation Maximization (SEM) approach is a stochastic version of the EM mixture estimation seen on section 19.1.4. It has been introduced by [22] to prevent convergence of the EM approach from local minima. It avoids the analytical maximization issued by integrating a stochastic sampling procedure in the estimation process. It induces an almost sure (a.s.) convergence to the algorithm.

From the initial two step formulation of the EM mixture estimation, the SEM may be decomposed into 3 steps:

1. E-step, calculates the expected membership values for each measurement vector to each classes.
2. S-step, performs a stochastic sampling of the membership vector to each classes, according to the membership values computed in the E-step.
3. M-step, updates the parameters of the membership probabilities (parameters to be defined through the class itk::Statistics::ModelComponentBase and its inherited classes).

The implementation of the SEM has been turned to a contextual SEM in the sense where the evaluation of the membership parameters is conditioned to membership values of the spatial neighborhood of each pixels.

The source code for this example can be found in the file
Examples/Learning/SEMModelEstimatorExample.cxx.

In this example, we present OTB’s implementation of SEM, through the class otb::SEMClassifier. This class performs a stochastic version of the EM algorithm, but instead of inheriting from itk::ExpectationMaximizationMixtureModelEstimator, we choosed to inherit from .html >itk::Statistics::ListSample< TSample >, in the same way as otb::SVMClassifier.

The program begins with otb::VectorImage and outputs itb::Image. Then appropriate header files have to be included:

#include "otbImage.h"
#include "otbVectorImage.h"
#include "otbImageFileWriter.h"

otb::SEMClassifier performs estimation of mixture to fit the initial histogram. Actually, mixture of Gaussian pdf can be performed. Those generic pdf are treated in otb::Statistics::ModelComponentBase. The Gaussian model is taken in charge with the class otb::Statistics::GaussianModelComponent.

#include "otbGaussianModelComponent.h"
#include "otbSEMClassifier.h"

Input/Output images type are define in a classical way. In fact, a itk::VariableLengthVector is to be considered for the templated MeasurementVectorType, which will be used in the ListSample interface.

typedef double PixelType;

typedef otb::VectorImage<PixelType, 2>  ImageType;

typedef otb::Image<unsigned char, 2>          OutputImageType;
typedef otb::ImageFileWriter<OutputImageType> WriterType;

Once the input image is opened, the classifier may be initialised by SmartPointer.

typedef otb::SEMClassifier<ImageType, OutputImageType> ClassifType;
ClassifType::Pointer classifier = ClassifType::New();

Then, it follows, classical initializations of the pipeline.

classifier->SetNumberOfClasses(numberOfClasses);
classifier->SetMaximumIteration(numberOfIteration);
classifier->SetNeighborhood(neighborhood);
classifier->SetTerminationThreshold(terminationThreshold);

When an initial segmentation is available, the classifier may use it as image (of type OutputImageType) or as a itk::SampleClassifier result (of type .html >itk::Statistics::MembershipSample< SampleType >).

if (fileNameImgInit != NULL)
{
}

By default, otb::SEMClassifier performs initialization of ModelComponentBase by as many instanciation of otb::Statistics::GaussianModelComponent as the number of classes to estimate in the mixture. Nevertheless, the user may add specific distribution into the mixture estimation. It is permited by the use of AddComponent for the given class number and the specific distribution.

typedef ClassifType::ClassSampleType ClassSampleType;
typedef otb::Statistics::GaussianModelComponent<ClassSampleType>
GaussianType;

for (int i = 0; i < numberOfClasses; ++i)

Once the pipeline is instanciated. The segmentation by itself may be launched by using the Update function.

try
{
classifier->Update();
}

The segmentation may outputs a result of type .html >itk::Statistics::MembershipSample< SampleType > as it is the case for the otb::SVMClassifier. But when using GetOutputImage the output is directly an Image.

Only for visualization purposes, we choose to rescale the image of classes before saving it to a file. We will use the itk::RescaleIntensityImageFilter for this purpose.

typedef itk::RescaleIntensityImageFilter<OutputImageType,

rescaler->SetOutputMinimum(itk::NumericTraits<unsigned char>::min());
rescaler->SetOutputMaximum(itk::NumericTraits<unsigned char>::max());

rescaler->SetInput(classifier->GetOutputImage());

WriterType::Pointer writer = WriterType::New();
writer->SetFileName(fileNameOut);
writer->SetInput(rescaler->GetOutput());
writer->Update();

Figure 19.12 shows the result of the SEM segmentation with 4 different classes and a contextual neighborhood of 3 pixels.

As soon as the segmentation is performed by an iterative stochastic process, it is worth verifying the output status: does the segmentation ends when it has converged or just at the limit of the iteration numbers.

std::cerr << "Program terminated with a ";
if (classifier->GetTerminationCode() ==
ClassifType::CONVERGED) std::cerr << "converged ";
else std::cerr << "not-converged ";
std::cerr << "code...\n";

The text output gives for each class the parameters of the pdf (e.g. mean of each component of the class and there covariance matrix, in the case of a Gaussian mixture model).

classifier->Print(std::cerr);

### 19.3 Support Vector Machines

Kernel based learning methods in general and the Support Vector Machines (SVM) in particular, have been introduced in the last years in learning theory for classification and regression tasks, [136]. SVM have been successfully applied to text categorization, [74], and face recognition, [105]. Recently, they have been successfully used for the classification of hyperspectral remote-sensing images, [15].

Simply stated, the approach consists in searching for the separating surface between 2 classes by the determination of the subset of training samples which best describes the boundary between the 2 classes. These samples are called support vectors and completely define the classification system. In the case where the two classes are nonlinearly separable, the method uses a kernel expansion in order to make projections of the feature space onto higher dimensionality spaces where the separation of the classes becomes linear.

#### 19.3.1 Mathematical formulation

This section reminds the basic principles of SVM learning and classification. A good tutorial on SVM can be found in, [17].

We have N samples represented by the couple (yi,xi),i = 1N where yi ∈{−1,+1} is the class label and xi n is the feature vector of dimension n. A classifier is a function

where α are the classifier parameters. The SVM finds the optimal separating hyperplane which fulfills the following constraints :
• The samples with labels +1 and 1 are on different sides of the hyperplane.
• The distance of the closest vectors to the hyperplane is maximised. These are the support vectors (SV) and this distance is called the margin.

The separating hyperplane has the equation

with w being its normal vector and x being any point of the hyperplane. The orthogonal distance to the origin is given by . Vectors located outside the hyperplane have either w x +b > 0 or w x +b < 0.

Therefore, the classifier function can be written as

The SVs are placed on two hyperplanes which are parallel to the optimal separating one. In order to find the optimal hyperplane, one sets w and b :

Since there must not be any vector inside the margin, the following constraint can be used:

which can be rewritten as

The orthogonal distances of the 2 parallel hyperplanes to the origin are and . Therefore the modulus of the margin is equal to and it has to be maximised.

Thus, the problem to be solved is:

• Find w and b which minimise
• under the constraint : yi(w xi+b) 1  i = 1N.

This problem can be solved by using the Lagrange multipliers with one multiplier per sample. It can be shown that only the support vectors will have a positive Lagrange multiplier.

In the case where the two classes are not exactly linearly separable, one can modify the constraints above by using

If ξi > 1, one considers that the sample is wrong. The function which has then to be minimised is w2 +C;, where C is a tolerance parameter. The optimisation problem is the same than in the linear case, but one multiplier has to be added for each new constraint ξi 0.

If the decision surface needs to be non-linear, this solution cannot be applied and the kernel approach has to be adopted.

One drawback of the SVM is that, in their basic version, they can only solve two-class problems. Some works exist in the field of multi-class SVM (see [5141], and the comparison made by [60]), but they are not used in our system.

You have to be aware that to achieve better convergence of the algorithm it is strongly advised to normalize feature vector components in the [1;1] interval.

For problems with N > 2 classes, one can choose either to train N SVM (one class against all the others), or to train N ×(N 1) SVM (one class against each of the others). In the second approach, which is the one that we use, the final decision is taken by choosing the class which is most often selected by the whole set of SVM.

#### 19.3.2 Learning With PointSets

The source code for this example can be found in the file
Examples/Learning/SVMPointSetModelEstimatorExample.cxx.

This example illustrates the use of the otb::SVMPointSetModelEstimator in order to perform the SVM learning from an itk::PointSet data structure.

The first step required to use this filter is to include its header file.

#include "otbSVMPointSetModelEstimator.h"

In the framework of supervised learning and classification, we will always use feature vectors for the characterization of the classes. On the other hand, the class labels are scalar values. Here, we start by defining the type of the features as the PixelType, which will be used to define the feature VectorType. We also declare the type for the labels.

typedef float                  PixelType;
typedef std::vector<PixelType> VectorType;
typedef int                    LabelPixelType;

We can now proceed to define the point sets used for storing the features and the labels.

typedef itk::PointSet<VectorType,  Dimension> FeaturePointSetType;

typedef itk::PointSet<LabelPixelType,  Dimension> LabelPointSetType;

FeaturePointSetType::Pointer fPSet = FeaturePointSetType::New();
LabelPointSetType::Pointer   lPSet = LabelPointSetType::New();

We will need to get access to the data stored in the point sets, so we define the appropriate for the points and the points containers used by the point sets (see the section 5.2 for more information oin haw to use point sets).

typedef FeaturePointSetType::PointType FeaturePointType;
typedef LabelPointSetType::PointType   LabelPointType;

typedef FeaturePointSetType::PointsContainer FeaturePointsContainer;
typedef LabelPointSetType::PointsContainer   LabelPointsContainer;

FeaturePointsContainer::Pointer fCont = FeaturePointsContainer::New();
LabelPointsContainer::Pointer   lCont = LabelPointsContainer::New();

We need now to build the training set for the SVM learning. In this simple example, we will build a SVM who classes points depending on which side of the line x = y they are located. We start by generating 500 random points.

int lowest = 0;
int range = 1000;

for (unsigned int pointId = 0; pointId < 500; pointId++)
{

FeaturePointType fP;
LabelPointType   lP;

int x_coord = lowest + static_cast<int>(range  (rand() / (RAND_MAX + 1.0)));
int y_coord = lowest + static_cast<int>(range  (rand() / (RAND_MAX + 1.0)));

We set the coordinates of the points. They are the same for the feature vector and for the label.

fP[0] = x_coord;
fP[1] = y_coord;

lP[0] = x_coord;
lP[1] = y_coord;

We push the features in the vector after a normalization which is useful for SVM convergence.

VectorType feature;
feature.push_back(static_cast<PixelType>((x_coord  1.0 - lowest) / range));
feature.push_back(static_cast<PixelType>((y_coord  1.0 - lowest) / range));

We decide on the label for each point.

LabelPixelType label;

if (x_coord < y_coord) label = -1;
else label = 1;

And we insert the points in the points containers.

fCont->InsertElement(pointId, fP);
fPSet->SetPointData(pointId, feature);

lCont->InsertElement(pointId, lP);
lPSet->SetPointData(pointId, label);

}

After the loop, we set the points containers to the point sets.

fPSet->SetPoints(fCont);
lPSet->SetPoints(lCont);

Up to now, we have only prepared the data for the SVM learning. We can now create the SVM model estimator. This class is templated over the feature and the label point set types.

typedef otb::SVMPointSetModelEstimator<FeaturePointSetType,
LabelPointSetType>   EstimatorType;

EstimatorType::Pointer estimator = EstimatorType::New();

The next step consists in setting the point sets for the estimator and the number of classes for the model. The feture point set is set using the SetInputPointSet and the label point set is set with the SetTrainingPointSet method.

estimator->SetInputPointSet(fPSet);
estimator->SetTrainingPointSet(lPSet);

The model estimation is triggered by calling the Update method.

estimator->Update();

Finally, we can save the result of the learning to a file.

estimator->SaveModel("svm_model.svm");

The otb::otbSVMModel class provides several accessors in order to get some information about the result of the learning step. For instance, one can get the number of support vectors kept to define the separation surface by using the GetNumberOfSupportVectors(). This can be very useful to detect some kind of overlearning (the number of support vectors is close to the number of examples). One can also get the SVs themselves by calling the GetSupportVectors(). The α values for the support vectors can be accessed by using the GetAlpha() method. Finally the Evaluate() method will return the result of the classification of a sample and the EvaluateHyperplaneDistance() will return the distance of the sample to the separating surface (or surfaces in the case of multi-class problems).

#### 19.3.3 PointSet Classification

The source code for this example can be found in the file
Examples/Learning/SVMPointSetClassificationExample.cxx.

This example illustrates the use of the otb::SVMClassifier class for performing SVM classification on pointsets. The first thing to do is include the header file for the class. Since the otb::SVMClassifier takes itk::ListSamples as input, the class itk::PointSetToListAdaptor is needed.

We start by including the needed header files.

#include "itkListSample.h"
#include "otbSVMClassifier.h"

In the framework of supervised learning and classification, we will always use feature vectors for the characterization of the classes. On the other hand, the class labels are scalar values. Here, we start by defining the type of the features as the PixelType, which will be used to define the feature VectorType. We also declare the type for the labels.

typedef float InputPixelType;

typedef std::vector<InputPixelType> InputVectorType;
typedef int                         LabelPixelType;

We can now proceed to define the point sets used for storing the features and the labels.

typedef itk::PointSet<InputVectorType,  Dimension> MeasurePointSetType;

typedef itk::PointSet<LabelPixelType,  Dimension> LabelPointSetType;

We will need to get access to the data stored in the point sets, so we define the appropriate for the points and the points containers used by the point sets (see the section 5.2 for more information on how to use point sets).

typedef MeasurePointSetType::PointType MeasurePointType;
typedef LabelPointSetType::PointType   LabelPointType;

typedef MeasurePointSetType::PointsContainer MeasurePointsContainer;
typedef LabelPointSetType::PointsContainer   LabelPointsContainer;

MeasurePointSetType::Pointer    tPSet = MeasurePointSetType::New();
MeasurePointsContainer::Pointer tCont = MeasurePointsContainer::New();

We need now to build the test set for the SVM. In this simple example, we will build a SVM who classes points depending on which side of the line x = y they are located. We start by generating 500 random points.

int lowest = 0;
int range = 1000;

for (pointId = 0; pointId < 100; pointId++)
{

MeasurePointType tP;

int x_coord = lowest + static_cast<int>(range  (rand() / (RAND_MAX + 1.0)));
int y_coord = lowest + static_cast<int>(range  (rand() / (RAND_MAX + 1.0)));

std::cout << "coords : " << x_coord << " " << y_coord << std::endl;
tP[0] = x_coord;
tP[1] = y_coord;

We push the features in the vector after a normalization which is useful for SVM convergence.

InputVectorType measure;
measure.push_back(static_cast<InputPixelType>((x_coord  1.0 -
lowest) / range));
measure.push_back(static_cast<InputPixelType>((y_coord  1.0 -
lowest) / range));

And we insert the points in the points container.

tCont->InsertElement(pointId, tP);
tPSet->SetPointData(pointId, measure);

}

After the loop, we set the points container to the point set.

tPSet->SetPoints(tCont);

Once the pointset is ready, we must transform it to a sample which is compatible with the classification framework. We will use a itk::Statistics::PointSetToListAdaptor for this task. This class is templated over the point set type used for storing the measures.

SampleType;
SampleType::Pointer sample = SampleType::New();

After instantiation, we can set the point set as an imput of our sample adaptor.

sample->SetPointSet(tPSet);

Now, we need to declare the SVM model which is to be used by the classifier. The SVM model is templated over the type of value used for the measures and the type of pixel used for the labels.

typedef otb::SVMModel<SampleType::MeasurementVectorType::ValueType,
LabelPixelType> ModelType;

ModelType::Pointer model = ModelType::New();

After instantiation, we can load a model saved to a file (see section 19.3.2 for an example of model estimation and storage to a file).

We have now all the elements to create a classifier. The classifier is templated over the sample type (the type of the data to be classified) and the label type (the type of the output of the classifier).

typedef otb::SVMClassifier<SampleType, LabelPixelType> ClassifierType;

ClassifierType::Pointer classifier = ClassifierType::New();

We set the classifier parameters : number of classes, SVM model, the sample data. And we trigger the classification process by calling the Update method.

int numberOfClasses = model->GetNumberOfClasses();
classifier->SetNumberOfClasses(numberOfClasses);
classifier->SetModel(model);
classifier->SetSample(sample.GetPointer());
classifier->Update();

After the classification step, we usually want to get the results. The classifier gives an output under the form of a sample list. This list supports the classical STL iterators.

ClassifierType::OutputType membershipSample =
classifier->GetOutput();

ClassifierType::OutputType::ConstIterator m_iter =
membershipSample->Begin();
ClassifierType::OutputType::ConstIterator m_last =
membershipSample->End();

We will iterate through the list, get the labels and compute the classification error.

double error = 0.0;
pointId = 0;
while (m_iter != m_last)
{

We get the label for each point.

ClassifierType::ClassLabelType label = m_iter.GetClassLabel();

And we compare it to the corresponding one of the test set.

InputVectorType measure;

tPSet->GetPointData(pointId, &measure);

ClassifierType::ClassLabelType expectedLabel;
if (measure[0] < measure[1]) expectedLabel = -1;
else expectedLabel = 1;

double dist = fabs(measure[0] - measure[1]);

if (label != expectedLabel) error++;

std::cout << int(label) << "/" << int(expectedLabel) << " --- " << dist <<
std::endl;

++pointId;
++m_iter;
}

std::cout << "Error = " << error / pointId << " % " << std::endl;

#### 19.3.4 Learning With Images

The source code for this example can be found in the file
Examples/Learning/SVMImageModelEstimatorExample.cxx.

This example illustrates the use of the otb::SVMImageModelEstimator class. This class allows the estimation of a SVM model (supervised learning) from a feature image and an image of labels. In this example, we will train an SVM to separate between water and non-water pixels by using the RGB values only. The images used for this example are shown in figure 19.13.

The first thing to do is include the header file for the class.

#include "otbSVMImageModelEstimator.h"

We define the types for the input and training images. Even if the input image will be an RGB image, we can read it as a 3 component vector image. This simplifies the interfacing with OTB’s SVM framework.

typedef unsigned char InputPixelType;
const unsigned int Dimension = 2;

typedef otb::VectorImage<InputPixelType,  Dimension> InputImageType;

typedef otb::Image<InputPixelType,  Dimension> TrainingImageType;

The otb::SVMImageModelEstimator class is templated over the input (features) and the training (labels) images.

typedef otb::SVMImageModelEstimator<InputImageType,
TrainingImageType>   EstimatorType;

As usual, we define the readers for the images.

We read the images. It is worth to note that, in order to ensure the pipeline coherence, the output of the objects which preceed the model estimator in the pipeline, must be up to date, so we call the corresponding Update methods.

We can now instantiate the model estimator and set its parameters.

EstimatorType::Pointer svmEstimator = EstimatorType::New();

The model estimation procedure is triggered by calling the estimator’s Update method.

svmEstimator->Update();

Finally, the estimated model can be saved to a file for later use.

svmEstimator->SaveModel(outputModelFileName);

#### 19.3.5 Image Classification

The source code for this example can be found in the file
Examples/Learning/SVMImageClassificationExample.cxx.

This example illustrates the use of the otb::SVMClassifier class for performing SVM classification on images. In this example, we will use an SVM model estimated in the example of section 19.3.4 to separate between water and non-water pixels by using the RGB values only. The images used for this example are shown in figure 19.13. The first thing to do is include the header file for the class. Since the otb::SVMClassifier takes itk::ListSamples as input, the class itk::PointSetToListAdaptor is needed.

In the framework of supervised learning and classification, we will always use feature vectors for the characterization of the classes. On the other hand, the class labels are scalar values. Here, we start by defining the type of the features as the PixelType, which will be used to define the feature VectorType. We also declare the type for the labels.

typedef double                 PixelType;
typedef std::vector<PixelType> VectorType;
typedef int                    LabelPixelType;

We can now proceed to define the image type used for storing the features. We also define the reader.

typedef otb::Image<itk::FixedArray<PixelType, 3>,
Dimension>          InputImageType;

We can now read the image by calling the Update method of the reader.

The image has now to be transformed to a sample which is compatible with the classification framework. We will use a itk::Statistics::ImageToListAdaptor for this task. This class is templated over the image type used for storing the measures.

SampleType::Pointer sample = SampleType::New();

After instantiation, we can set the image as an imput of our sample adaptor.

Now, we need to declare the SVM model which is to be used by the classifier. The SVM model is templated over the type of value used for the measures and the type of pixel used for the labels.

typedef otb::SVMModel<PixelType, LabelPixelType> ModelType;

ModelType::Pointer model = ModelType::New();

After instantiation, we can load a model saved to a file (see section 19.3.4 for an example of model estimation and storage to a file.

We have now all the elements to create a classifier. The classifier is templated over the sample type (the type of the data to be classified) and the label type (the type of the output of the classifier).

typedef otb::SVMClassifier<SampleType, LabelPixelType> ClassifierType;

ClassifierType::Pointer classifier = ClassifierType::New();

We set the classifier parameters : number of classes, SVM model, the sample data. And we trigger the classification process by calling the Update method.

int numberOfClasses = model->GetNumberOfClasses();
classifier->SetNumberOfClasses(numberOfClasses);
classifier->SetModel(model);
classifier->SetSample(sample.GetPointer());
classifier->Update();

After the classification step, we usually want to get the results. The classifier gives an output under the form of a sample list. This list supports the classical STL iterators. Therefore, we will create an output image and fill it up with the results of the classification. The pixel type of the output image is the same as the one used for the labels.

typedef ClassifierType::ClassLabelType         OutputPixelType;
typedef otb::Image<OutputPixelType, Dimension> OutputImageType;

OutputImageType::Pointer outputImage = OutputImageType::New();

We allocate the memory for the output image using the information from the input image.

typedef itk::Index<Dimension>       myIndexType;
typedef itk::Size<Dimension>        mySizeType;
typedef itk::ImageRegion<Dimension> myRegionType;

mySizeType size;

myIndexType start;
start[0] = 0;
start[1] = 0;

myRegionType region;
region.SetIndex(start);
region.SetSize(size);

outputImage->SetRegions(region);
outputImage->Allocate();

We can now declare the interators on the list that we get at the output of the classifier as well as the iterator to fill the output image.

ClassifierType::OutputType membershipSample =
classifier->GetOutput();
ClassifierType::OutputType::ConstIterator m_iter =
membershipSample->Begin();
ClassifierType::OutputType::ConstIterator m_last =
membershipSample->End();

typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
OutputIteratorType outIt(outputImage,
outputImage->GetBufferedRegion());

outIt.GoToBegin();

We will iterate through the list, get the labels and assign pixel values to the output image.

while (m_iter != m_last && !outIt.IsAtEnd())
{
outIt.Set(m_iter.GetClassLabel());
++m_iter;
++outIt;
}

Only for visualization purposes, we choose to rescale the image of classes before saving it to a file. We will use the itk::RescaleIntensityImageFilter for this purpose.

typedef otb::Image<unsigned char, Dimension> FileImageType;

typedef itk::RescaleIntensityImageFilter<OutputImageType,

rescaler->SetOutputMinimum(itk::NumericTraits<unsigned char>::min());
rescaler->SetOutputMaximum(itk::NumericTraits<unsigned char>::max());

rescaler->SetInput(outputImage);

We can now create an image file writer and save the image.

typedef otb::ImageFileWriter<FileImageType> WriterType;

WriterType::Pointer writer = WriterType::New();

writer->SetFileName(outputFilename);
writer->SetInput(rescaler->GetOutput());

writer->Update();

Figure 19.14 shows the result of the SVM classification.

The source code for this example can be found in the file
Examples/Learning/SVMImageEstimatorClassificationMultiExample.cxx.

This example illustrates the OTB’s multi-class SVM capabilities. The theory behind this kind of classification is out of the scope of this guide. In OTB, the multi-class SVM classification is used in the same way as the two-class one. Figure 19.15 shows the image to be classified and the associated ground truth, which is composed of 4 classes.

The following header files are needed for the program:

#include "otbSVMImageModelEstimator.h"
#include "itkListSample.h"
#include "otbSVMClassifier.h"

We define the types for the input and training images. Even if the input image will be an RGB image, we can read it as a 3 component vector image. This simplifies the interfacing with OTB’s SVM framework.

typedef unsigned short InputPixelType;
const unsigned int Dimension = 2;

typedef otb::VectorImage<InputPixelType, Dimension> InputImageType;

typedef otb::Image<InputPixelType,  Dimension> TrainingImageType;

The otb::SVMImageModelEstimator class is templated over the input (features) and the training (labels) images.

typedef otb::SVMImageModelEstimator<InputImageType,
TrainingImageType>   EstimatorType;

As usual, we define the readers for the images.

We read the images. It is worth to note that, in order to ensure the pipeline coherence, the output of the objects which preceed the model estimator in the pipeline, must be up to date, so we call the corresponding Update methods.

We can now instantiate the model estimator and set its parameters.

EstimatorType::Pointer svmEstimator = EstimatorType::New();

The model estimation procedure is triggered by calling the estimator’s Update method.

svmEstimator->Update();

We can now proceed to the image classification. We start by declaring the type of the image to be classified. ITK’s classification framework needs the type of the pixel to be of fixed type, so we declare the following types.

typedef otb::Image<itk::FixedArray<InputPixelType, 3>,
Dimension>          ClassifyImageType;

We can now read the image by calling the Update method of the reader.

The image has now to be transformed to a sample which is compatible with the classification framework. We will use a itk::Statistics::ImageToListAdaptor for this task. This class is templated over the image type used for storing the measures.

SampleType::Pointer sample = SampleType::New();

After instantiation, we can set the image as an imput of our sample adaptor.

Now, we need to declare the SVM model which is to be used by the classifier. The SVM model is templated over the type of value used for the measures and the type of pixel used for the labels. The model is obtained from the model estimator by calling the GetModel method.

typedef InputPixelType LabelPixelType;

typedef otb::SVMModel<InputPixelType, LabelPixelType> ModelType;

ModelType::Pointer model = svmEstimator->GetModel();

We have now all the elements to create a classifier. The classifier is templated over the sample type (the type of the data to be classified) and the label type (the type of the output of the classifier).

typedef otb::SVMClassifier<SampleType, LabelPixelType> ClassifierType;

ClassifierType::Pointer classifier = ClassifierType::New();

We set the classifier parameters : number of classes, SVM model, the sample data. And we trigger the classification process by calling the Update method.

int numberOfClasses = model->GetNumberOfClasses();
classifier->SetNumberOfClasses(numberOfClasses);
classifier->SetModel(model);
classifier->SetSample(sample.GetPointer());
classifier->Update();

After the classification step, we usually want to get the results. The classifier gives an output under the form of a sample list. This list supports the classical STL iterators. Therefore, we will create an output image and fill it up with the results of the classification. The pixel type of the output image is the same as the one used for the labels.

typedef ClassifierType::ClassLabelType         OutputPixelType;
typedef otb::Image<OutputPixelType, Dimension> OutputImageType;

OutputImageType::Pointer outputImage = OutputImageType::New();

We allocate the memory for the output image using the information from the input image.

typedef itk::Index<Dimension>       myIndexType;
typedef itk::Size<Dimension>        mySizeType;
typedef itk::ImageRegion<Dimension> myRegionType;

mySizeType size;

myIndexType start;
start[0] = 0;
start[1] = 0;

myRegionType region;
region.SetIndex(start);
region.SetSize(size);

outputImage->SetRegions(region);
outputImage->Allocate();
std::cout << "---" << std::endl;

We can now declare the interators on the list that we get at the output of the classifier as well as the iterator to fill the output image.

ClassifierType::OutputType membershipSample =
classifier->GetOutput();
ClassifierType::OutputType::ConstIterator m_iter =
membershipSample->Begin();
ClassifierType::OutputType::ConstIterator m_last =
membershipSample->End();

typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
OutputIteratorType outIt(outputImage,
outputImage->GetBufferedRegion());

outIt.GoToBegin();

We will iterate through the list, get the labels and assign pixel values to the output image.

while (m_iter != m_last && !outIt.IsAtEnd())
{
outIt.Set(m_iter.GetClassLabel());
++m_iter;
++outIt;
}
std::cout << "---" << std::endl;

Only for visualization purposes, we choose a color mapping to the image of classes before saving it to a file. The itk::Functor::ScalarToRGBPixelFunctor class is a special function object designed to hash a scalar value into an itk::RGBPixel. Plugging this functor into the itk::UnaryFunctorImageFilter creates an image filter for that converts scalar images to RGB images.

typedef itk::RGBPixel<unsigned char> RGBPixelType;
typedef otb::Image<RGBPixelType, 2>  RGBImageType;
typedef itk::Functor::ScalarToRGBPixelFunctor<unsigned long>
ColorMapFunctorType;
typedef itk::UnaryFunctorImageFilter<OutputImageType,
RGBImageType,
ColorMapFunctorType> ColorMapFilterType;
ColorMapFilterType::Pointer colormapper = ColorMapFilterType::New();

colormapper->SetInput(outputImage);

We can now create an image file writer and save the image.

typedef otb::ImageFileWriter<RGBImageType> WriterRescaledType;

WriterRescaledType::Pointer writerRescaled = WriterRescaledType::New();

writerRescaled->SetFileName(outputRescaledImageFileName);
writerRescaled->SetInput(colormapper->GetOutput());

writerRescaled->Update();

Figure 19.16 shows the result of the SVM classification.

#### 19.3.6 Generic Kernel SVM

OTB has developed a specific interface for user-defined kernels. A function k(,) is considered to be a kernel when:

 ∀g(⋅) ∈2(Rn) so that ∫ g(x)2dx be finite, (19.3) then ∫ k(x,y)g(x)g(y)dxdy ≥ 0,
which is known as the Mercer condition.

When defined through the OTB, a kernel is a class that inherits from GenericKernelFunctorBase. Several virtual functions have to be overloaded:

• The Evaluate function, which implements the behavior of the kernel itself. For instance, the classical linear kernel could be re-implemented with:
double
MyOwnNewKernel
::Evaluate ( const svm_node ⋆ x, const svm_node ⋆ y,
const svm_parameter & param ) const
{
return this->dot(x,y);
}

This simple example shows that the classical dot product is already implemented into otb::GenericKernelFunctorBase::dot() as a protected function.

• The Update() function which synchronizes local variables and their integration into the initial SVM procedure. The following examples will show the way to use it.

Some pre-defined generic kernels have already been implemented in OTB:

##### Learning with User Defined Kernels

The source code for this example can be found in the file
Examples/Learning/SVMGenericKernelImageModelEstimatorExample.cxx.

This example illustrates the modifications to be added to the use of otb::SVMImageModelEstimator in order to add a user defined kernel. This initial program has been explained in section 19.3.4.

The first thing to do is to include the header file for the new kernel.

#include "otbSVMImageModelEstimator.h"
#include "otbMixturePolyRBFKernelFunctor.h"

Once the otb::SVMImageModelEstimator is instanciated, it is possible to add the new kernel and its parameters.

Then in addition to the initial code:

EstimatorType::Pointer svmEstimator = EstimatorType::New();

svmEstimator->SetSVMType(C_SVC);

The instanciation of the kernel is to be implemented. The kernel which is used here is a linear combination of a polynomial kernel and an RBF one. It is written as

with k1(x,y) = d and k2(x,y) = exp. Then, the specific parameters of this kernel are:
• Mixture (μ),
• GammaPoly (γ1),
• CoefPoly (c0),
• DegreePoly (d),
• GammaRBF (γ2).

Their instanciations are achieved through the use of the SetValue function.

otb::MixturePolyRBFKernelFunctor myKernel;
myKernel.SetValue("Mixture", 0.5);
myKernel.SetValue("GammaPoly", 1.0);
myKernel.SetValue("CoefPoly", 0.0);
myKernel.SetValue("DegreePoly", 1);
myKernel.SetValue("GammaRBF", 1.5);
myKernel.Update();

Once the kernel’s parameters are affected and the kernel updated, the connection to otb::SVMImageModelEstimator takes place here.

svmEstimator->SetKernelFunctor(&myKernel);
svmEstimator->SetKernelType(GENERIC);

The model estimation procedure is triggered by calling the estimator’s Update method.

svmEstimator->Update();

The rest of the code remains unchanged...

svmEstimator->SaveModel(outputModelFileName);

In the file outputModelFileName a specific line will appear when using a generic kernel. It gives the name of the kernel and its parameters name and value.

##### Classification with user defined kernel

The source code for this example can be found in the file
Examples/Learning/SVMGenericKernelImageClassificationExample.cxx.

This example illustrates the modifications to be added to use the otb::SVMClassifier class for performing SVM classification on images with a user-defined kernel. In this example, we will use an SVM model estimated in the previous section to separate between water and non-water pixels by using the RGB values only. The images used for this example are shown in figure 19.13. The first thing to do is include the header file for the class as well as the header of the appropriated kernel to be used.

#include "otbSVMClassifier.h"
#include "otbMixturePolyRBFKernelFunctor.h"

We need to declare the SVM model which is to be used by the classifier. The SVM model is templated over the type of value used for the measures and the type of pixel used for the labels.

typedef otb::SVMModel<PixelType, LabelPixelType> ModelType;
ModelType::Pointer model = ModelType::New();

After instantiation, we can load a model saved to a file (see section 19.3.4 for an example of model estimation and storage to a file).

When using a user defined kernel, an explicit instanciation has to be performed.

otb::MixturePolyRBFKernelFunctor myKernel;
model->SetKernelFunctor(&myKernel);

Then, the rest of the classification program remains unchanged.

#### 19.3.7 Multi-band, streamed classification

The source code for this example can be found in the file
Examples/Classification/SVMImageClassifierExample.cxx.

In previous examples, we have used the otb::SVMClassifier, which uses the ITK classification framework. This good for compatibility with the ITK framework, but introduces the limitations of not being able to use streaming and being able to know at compilation time the number of bands of the image to be classified. In OTB we have avoided this limitation by developing the otb::SVMImageClassificationFilter. In this example we will illustrate its use. We start by including the appropriate header file.

#include "otbSVMImageClassificationFilter.h"

We will assume double precision input images and will also define the type for the labeled pixels.

const unsigned int Dimension = 2;
typedef double         PixelType;
typedef unsigned short LabeledPixelType;

Our classifier will be generic enough to be able to process images with any number of bands. We read the images as otb::VectorImages. The labeled image will be a scalar image.

typedef otb::VectorImage<PixelType, Dimension>  ImageType;
typedef otb::Image<LabeledPixelType, Dimension> LabeledImageType;

We can now define the type for the classifier filter, which is templated over its input and output image types.

typedef otb::SVMImageClassificationFilter
<ImageType, LabeledImageType>  ClassificationFilterType;
typedef ClassificationFilterType::ModelType ModelType;

And finally, we define the reader and the writer. Since the images to classify can be very big, we will use a streamed writer which will trigger the streaming ability of the classifier.

typedef otb::ImageFileWriter<LabeledImageType> WriterType;

We instantiate the classifier and the reader objects and we set the existing SVM model obtained in a previous training step.

ClassificationFilterType::Pointer filter = ClassificationFilterType::New();

ModelType::Pointer model = ModelType::New();

filter->SetModel(model);

We plug the pipeline and trigger its execution by updating the output of the writer.

WriterType::Pointer writer = WriterType::New();
writer->SetInput(filter->GetOutput());
writer->SetFileName(outfname);
writer->Update();

#### 19.3.8 Classification map regularization

The source code for this example can be found in the file
Examples/Classification/ClassificationMapRegularizationExample.cxx.

After having generated a classification map, it is possible to regularize such a labeled image in order to obtain more homogeneous areas, which makes the interpretation of its classes easier. For this purpose, the otb::NeighborhoodMajorityVotingImageFilter was implemented. Like a morphological filter, this filter uses majority voting in a ball shaped neighborhood in order to set each pixel of the classification map to the more representative label value in its neighborhood.

In this example we will illustrate its use. We start by including the appropriate header file.

#include "otbNeighborhoodMajorityVotingImageFilter.h"

Since the input image is a classification map, we will assume a single band input image for which each pixel value is a label coded on 8 bits as an integer between 0 and 255.

typedef unsigned char IOLabelPixelType; // 8 bits
const unsigned int Dimension = 2;

Thus, both input and output images are single band labeled images, which are composed of the same type of pixels in this example (unsigned char).

typedef otb::Image<IOLabelPixelType, Dimension> IOLabelImageType;

We can now define the type for the neighborhood majority voting filter, which is templated over its input and output images types and over its structuring element type. Choosing only the input image type in the template of this filter induces that, both input and output images types are the same and that the structuring element is a ball ( itk::BinaryBallStructuringElement).

// Neighborhood majority voting filter type
typedef otb::NeighborhoodMajorityVotingImageFilter<IOLabelImageType>
NeighborhoodMajorityVotingFilterType;

Since the otb::NeighborhoodMajorityVotingImageFilter is a neighborhood based image filter, it is necessary to set the structuring element which will be used for the majority voting process. By default, the structuring element is a ball ( itk::BinaryBallStructuringElement) with a radius defined by two sizes (respectively along X and Y). Thus, it is possible to handle anisotropic structuring elements such as ovals.

// Binary ball Structuring Element type
typedef NeighborhoodMajorityVotingFilterType::KernelType StructuringType;

Finally, we define the reader and the writer.

typedef otb::ImageFileWriter<IOLabelImageType> WriterType;

We instantiate the otb::NeighborhoodMajorityVotingImageFilter and the reader objects.

// Neighborhood majority voting filter
NeighborhoodMajorityVotingFilterType::Pointer NeighMajVotingFilter;
NeighMajVotingFilter = NeighborhoodMajorityVotingFilterType::New();

The ball shaped structuring element seBall is instantiated and its two radii along X and Y are initialized.

StructuringType seBall;

seBall.CreateStructuringElement();

Then, this ball shaped neighborhood is used as the kernel structuring element for the otb::NeighborhoodMajorityVotingImageFilter.

NeighMajVotingFilter->SetKernel(seBall);

Not classified input pixels are assumed to have the noDataValue label and will keep this label in the output image.

NeighMajVotingFilter->SetLabelForNoDataPixels(noDataValue);

Moreover, since the majority voting regularization may lead to not unique majority labels in the neighborhood, it is important to define which behaviour the filter must have in this case. For this purpose, a Boolean parameter is used in the filter to choose whether pixels with more than one majority class are set to undecidedValue (true), or to their Original labels (false = default value) in the output image.

NeighMajVotingFilter->SetLabelForUndecidedPixels(undecidedValue);

if (KeepOriginalLabelBoolStr.compare("true") == 0)
{
NeighMajVotingFilter->SetKeepOriginalLabelBool(true);
}
else
{
NeighMajVotingFilter->SetKeepOriginalLabelBool(false);
}

We plug the pipeline and trigger its execution by updating the output of the writer.

WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputFileName);
writer->SetInput(NeighMajVotingFilter->GetOutput());
writer->Update();

### 19.4 Kohonen’s Self Organizing Map

The Self Organizing Map, SOM, introduced by Kohonen is a non-supervised neural learning algorithm. The map is composed of neighboring cells which are in competition by means of mutual interactions and they adapt in order to match characteristic patterns of the examples given during the learning. The SOM is usually on a plane (2D).

The algorithm implements a nonlinear projection from a high dimensional feature space to a lower dimension space, usually 2D. It is able to find the correspondence between a set of structured data and a network of much lower dimension while keeping the topological relationships existing in the feature space. Thanks to this topological organization, the final map presents clusters and their relationships.

#### 19.4.1 The algorithm

Kohonen’s SOM is usually represented as an array of cells where each cell is, i, associated to a feature (or weight) vector mi = T n (figure 19.17).

A cell (or neuron) in the map is a good detector for a given input vector x = T n if the latter is close to the former. This distance between vectors can be represented by the scalar product xT mi, but for most of the cases other distances can be used, as for instance the Euclidean one. The cell having the weight vector closest to the input vector is called the winner.

##### Learning

The goal of the learning step is to get a map which is representative of an input example set. It is an iterative procedure which consists in passing each input example to the map, testing the response of each neuron and modifying the map to get it closer to the examples.

Algorithm 1 SOM learning:

1. t = 0.
2. Initialize the weight vectors of the map (randomly, for instance).
3. While t < number of iterations, do:
1. k = 0.
2. While k < number of examples, do:
1. Find the vector mi(t) which minimizes the distance d(xk,mi(t))
2. For a neighborhood Nc(t) around the winner cell, apply the transformation:
 (19.4)

3. k = k+1
3. t = t +1.

In 19.4, β(t) is a decreasing function with the geometrical distance to the winner cell. For instance:

 (19.5)

with β0(t) and σ(t) decreasing functions with time and r the cell coordinates in the output – map – space.

Therefore the algorithm consists in getting the map closer to the learning set. The use of a neighborhood around the winner cell allows the organization of the map into areas which specialize in the recognition of different patterns. This neighborhood also ensures that cells which are topologically close are also close in terms of the distance defined in the feature space.

#### 19.4.2 Building a color table

The source code for this example can be found in the file
Examples/Learning/SOMExample.cxx.

This example illustrates the use of the otb::SOM class for building Kohonen’s Self Organizing Maps.

We will use the SOM in order to build a color table from an input image. Our input image is coded with 3×8 bits and we would like to code it with only 16 levels. We will use the SOM in order to learn which are the 16 most representative RGB values of the input image and we will assume that this is the optimal color table for the image.

The first thing to do is include the header file for the class. We will also need the header files for the map itself and the activation map builder whose utility will be explained at the end of the example.

#include "otbSOMMap.h"
#include "otbSOM.h"
#include "otbSOMActivationBuilder.h"

Since the otb::SOM class uses a distance, we will need to include the header file for the one we want to use

#include "itkEuclideanDistance.h"

The Self Organizing Map itself is actually an N-dimensional image where each pixel contains a neuron. In our case, we decide to build a 2-dimensional SOM, where the neurons store RGB values with floating point precision.

const unsigned int Dimension = 2;
typedef double                                 PixelType;
typedef otb::VectorImage<PixelType, Dimension> ImageType;
typedef ImageType::PixelType                   VectorType;

The distance that we want to apply between the RGB values is the Euclidean one. Of course we could choose to use other type of distance, as for instance, a distance defined in any other color space.

typedef itk::Statistics::EuclideanDistance<VectorType> DistanceType;

We can now define the type for the map. The otb::SOMMap::class is templated over the neuron type – PixelType here –, the distance type and the number of dimensions. Note that the number of dimensions of the map could be different from the one of the images to be processed.

typedef otb::SOMMap<VectorType, DistanceType, Dimension> MapType;

We are going to perform the learning directly on the pixels of the input image. Therefore, the image type is defined using the same pixel type as we used for the map. We also define the type for the imge file reader.

Since the otb::SOM class works on lists of samples, it will need to access the input image through an adaptor. Its type is defined as follows:

typedef itk::Statistics::ListSample<VectorType> SampleListType;

We can now define the type for the SOM, which is templated over the input sample list and the type of the map to be produced and the two functors that hold the training behavior.

typedef otb::Functor::CzihoSOMLearningBehaviorFunctor
LearningBehaviorFunctorType;
typedef otb::Functor::CzihoSOMNeighborhoodBehaviorFunctor
NeighborhoodBehaviorFunctorType;
typedef otb::SOM<SampleListType, MapType,
LearningBehaviorFunctorType, NeighborhoodBehaviorFunctorType>
SOMType;

As an alternative to standart SOMType, one can decide to use an otb::PeriodicSOM, which behaves like otb::SOM but is to be considered to as a torus instead of a simple map. Hence, the neighborhood behavior of the winning neuron does not depend on its location on the map... otb::PeriodicSOM is defined in otbPeriodicSOM.h.

We can now start building the pipeline. The first step is to instantiate the reader and pass its output to the adaptor.

SampleListType::Pointer sampleList = SampleListType::New();

GetBufferedRegion());
imgIter.GoToBegin();

GetBufferedRegion());
imgIterEnd.GoToEnd();

do
{
sampleList->PushBack(imgIter.Get());
++imgIter;
}
while (imgIter != imgIterEnd);

We can now instantiate the SOM algorithm and set the sample list as input.

SOMType::Pointer som = SOMType::New();
som->SetListSample(sampleList);

We use a SOMType::SizeType array in order to set the sizes of the map.

SOMType::SizeType size;
size[0] = sizeX;
size[1] = sizeY;
som->SetMapSize(size);

The initial size of the neighborhood of each neuron is set in the same way.

The other parameters are the number of iterations, the initial and the final values for the learning rate – β – and the maximum initial value for the neurons (the map will be randomly initialized).

som->SetNumberOfIterations(nbIterations);
som->SetBetaInit(betaInit);
som->SetBetaEnd(betaEnd);
som->SetMaxWeight(static_cast<PixelType>(initValue));

Now comes the intialization of the functors.

LearningBehaviorFunctorType learningFunctor;
som->SetBetaFunctor(learningFunctor);

NeighborhoodBehaviorFunctorType neighborFunctor;
som->SetNeighborhoodSizeFunctor(neighborFunctor);
som->Update();

Finally, we set up the las part of the pipeline where the plug the output of the SOM into the writer. The learning procedure is triggered by calling the Update() method on the writer. Since the map is itself an image, we can write it to disk with an otb::ImageFileWriter.

Figure 19.18 shows the result of the SOM learning. Since we have performed a learning on RGB pixel values, the produced SOM can be interpreted as an optimal color table for the input image. It can be observed that the obtained colors are topologically organised, so similar colors are also close in the map. This topological organisation can be exploited to further reduce the number of coding levels of the pixels without performing a new learning: we can subsample the map to get a new color table. Also, a bilinear interpolation between the neurons can be used to increase the number of coding levels.

We can now compute the activation map for the input image. The activation map tells us how many times a given neuron is activated for the set of examples given to the map. The activation map is stored as a scalar image and an integer pixel type is usually enough.

typedef unsigned char OutputPixelType;

typedef otb::Image<OutputPixelType, Dimension> OutputImageType;
typedef otb::ImageFileWriter<OutputImageType>  ActivationWriterType;

In a similar way to the otb::SOM class the otb::SOMActivationBuilder is templated over the sample list given as input, the SOM map type and the activation map to be built as output.

typedef otb::SOMActivationBuilder<SampleListType, MapType,
OutputImageType> SOMActivationBuilderType;

We instantiate the activation map builder and set as input the SOM map build before and the image (using the adaptor).

SOMActivationBuilderType::Pointer somAct
= SOMActivationBuilderType::New();
somAct->SetInput(som->GetOutput());
somAct->SetListSample(sampleList);
somAct->Update();

The final step is to write the activation map to a file.

if (actMapFileName != NULL)
{
ActivationWriterType::Pointer actWriter = ActivationWriterType::New();
actWriter->SetFileName(actMapFileName);

The righthand side of figure 19.18 shows the activation map obtained.

#### 19.4.3 SOM Classification

The source code for this example can be found in the file
Examples/Learning/SOMClassifierExample.cxx.

This example illustrates the use of the otb::SOMClassifier class for performing a classification using an existing Kohonen’s Self Organizing. Actually, the SOM classification consists only in the attribution of the winner neuron index to a given feature vector.

We will use the SOM created in section 19.4.2 and we will assume that each neuron represents a class in the image.

The first thing to do is include the header file for the class.

#include "otbSOMClassifier.h"

As for the SOM learning step, we must define the types for the otb::SOMMap, and therefore, also for the distance to be used. We will also define the type for the SOM reader, which is actually an otb::ImageFileReader::which the appropiate image type.

typedef itk::Statistics::EuclideanDistance<PixelType>   DistanceType;
typedef otb::SOMMap<PixelType, DistanceType, Dimension> SOMMapType;

The classification will be performed by the otb::SOMClassifier::, which, as most of the classifiers, works on itk::Statistics::ListSamples. In order to be able to perform an image classification, we will need to use the itk::Statistics::ImageToListAdaptor which is templated over the type of image to be adapted. The SOMClassifier is templated over the sample type, the SOMMap type and the pixel type for the labels.

typedef itk::Statistics::ListSample<PixelType> SampleType;
typedef otb::SOMClassifier<SampleType, SOMMapType, LabelPixelType>
ClassifierType;

The result of the classification will be stored on an image and saved to a file. Therefore, we define the types needed for this step.

typedef otb::Image<LabelPixelType, Dimension> OutputImageType;
typedef otb::ImageFileWriter<OutputImageType> WriterType;

We can now start reading the input image and the SOM given as inputs to the program. We instantiate the readers as usual.

The conversion of the input data from image to list sample is easily done using the adaptor.

SampleType::Pointer sample = SampleType::New();

GetLargestPossibleRegion());

it.GoToBegin();

while (!it.IsAtEnd())
{
sample->PushBack(it.Get());
++it;
}

The classifier can now be instantiated. The input data is set by using the SetSample() method and the SOM si set using the SetMap() method. The classification is triggered by using the Update() method.

ClassifierType::Pointer classifier = ClassifierType::New();
classifier->SetSample(sample.GetPointer());
classifier->Update();

Once the classification has been performed, the sample list obtained at the output of the classifier must be converted into an image. We create the image as follows :

OutputImageType::Pointer outputImage = OutputImageType::New();
outputImage->Allocate();

We can now get a pointer to the classification result.

ClassifierType::OutputType membershipSample = classifier->GetOutput();

And we can declare the iterators pointing to the front and the back of the sample list.

ClassifierType::OutputType::ConstIterator m_iter =  membershipSample->Begin();
ClassifierType::OutputType::ConstIterator m_last =  membershipSample->End();

We also declare an itk::ImageRegionIterator::in order to fill the output image with the class labels.

typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;

OutputIteratorType outIt(outputImage, outputImage->GetLargestPossibleRegion());

We iterate through the sample list and the output image and assign the label values to the image pixels.

outIt.GoToBegin();

while (m_iter != m_last && !outIt.IsAtEnd())
{
outIt.Set(m_iter.GetClassLabel());
++m_iter;
++outIt;
}

Finally, we write the classified image to a file.

WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputFilename);
writer->SetInput(outputImage);
writer->Update();

Figure 19.19 shows the result of the SOM classification.

#### 19.4.4 Multi-band, streamed classification

The source code for this example can be found in the file
Examples/Classification/SOMImageClassificationExample.cxx.

In previous examples, we have used the otb::SOMClassifier, which uses the ITK classification framework. This good for compatibility with the ITK framework, but introduces the limitations of not being able to use streaming and being able to know at compilation time the number of bands of the image to be classified. In OTB we have avoided this limitation by developing the otb::SOMImageClassificationFilter. In this example we will illustrate its use. We start by including the appropriate header file.

#include "otbSOMImageClassificationFilter.h"

We will assume double precision input images and will also define the type for the labeled pixels.

const unsigned int Dimension = 2;
typedef double         PixelType;
typedef unsigned short LabeledPixelType;

Our classifier will be genric enough to be able to process images with any number of bands. We read the images as otb::VectorImages. The labeled image will be a scalar image.

typedef otb::VectorImage<PixelType, Dimension>  ImageType;
typedef otb::Image<LabeledPixelType, Dimension> LabeledImageType;

We can now define the type for the classifier filter, which is templated over its input and output image types and the SOM type.

typedef otb::SOMMap<ImageType::PixelType> SOMMapType;
typedef otb::SOMImageClassificationFilter<ImageType,
LabeledImageType,
SOMMapType>
ClassificationFilterType;

And finally, we define the readers (for the input image and theSOM) and the writer. Since the images, to classify can be very big, we will use a streamed writer which will trigger the streaming ability of the classifier.

typedef otb::ImageFileWriter<LabeledImageType> WriterType;

We instantiate the classifier and the reader objects and we set the existing SOM obtained in a previous training step.

ClassificationFilterType::Pointer filter = ClassificationFilterType::New();