# Classification of Cell Images Using MPEG-7-influenced Descriptors and
Support Vector Machines in Cell Morphology

Tobias Abenius
<EMAIL_ADDRESS>
###### Abstract

Counting and classifying blood cells is an important diagnostic tool in
medicine. Support Vector Machines are increasingly popular and efficient and
could replace artificial neural network systems. Here a method to classify
blood cells is proposed using SVM. A set of statistics on images are
implemented in C++. The MPEG-7 descriptors Scalable Color Descriptor, Color
Structure Descriptor, Color Layout Descriptor and Homogeneous Texture
Descriptor are extended in size and combined with textural features
corresponding to textural properties perceived visually by humans. From a set
of images of human blood cells these statistics are collected. A SVM is
implemented and trained to classify the cell images. The cell images come from
a CellaVision™ DM-96 machine which classify cells from images from microscopy.
The output images and classification of the CellaVision™ machine is taken as
ground truth, a truth that is 90-95% correct. The problem is divided in two —
the primary and the simplified. The primary problem is to classify the same
classes as the CellaVision™ machine. The simplified problem is to differ
between the five most common types of white blood cells. An encouraging result
is achieved in both cases — error rates of 10.8% and 3.1% — considering that
the SVM is misled by the errors in ground truth. Conclusion is that further
investigation of performance is worthwhile.

###### Sammanfattning

Att räkna och klassificera blodceller är ett viktigt diagnostiskt redskap inom
läkarvetenskapen. Support Vector Machines är effektiva, ökar i popularitet och
kan ersätta artificiella neurala nätverkssystem. Här föreslås en metod för att
klassificera blodceller m.h.a. SVM. En mängd statistika på bilder
implementeras i C++. De s.k. MPEG-7 descriptors Scalable Color Descriptor,
Color Structure Descriptor, Color Layout Descriptor och Homogeneous Texture
Descriptor utvidgas i storlek och kombineras med textur-mått motsvarande
textur-egenskaper som uppfattas visuellt av människor. Från en mängd bilder av
mänskliga blodceller samlas dessa mått. En SVM implementeras och tränas att
klassificera cellbilderna. Cellbilderna kommer från en CellaVision™ DM-96 som
klassificerar celler från mikroskoperade bilder. Bilderna och dess klasser
från en CellaVision™ DM-96-maskin tas som facit, ett facit som är 90-95%
korrekt. Problemet delas i två — det primära och det förenklade. Det primära
problemet är att skilja mellan de klasser som CellaVision™s maskin gör. Det
förenklade problemet är att skilja mellan de fem vanligaste typerna av vita
blodkroppar. Ett glädjande resultat uppnås i båda fallen — felfrekvenser om
10,8% och 3,1% — med tanke på att SVM missleddes av felen i det tagna facitet.
Slutsatsen är att vidare studier angående prestanda är lönsamma.

veelo companion subsection

firstpageempty firstpage  Examensarbete för 30 hp
Institutionen för datavetenskap, Naturvetenskapliga fakulteten, Lunds
universitet

Thesis for a diploma in Computer Science, 30 ECTS credits
Department of Computer Science, Faculty of Science, Lund University
firstpage0.85 adjustwidth*-

adjustwidth*-

—

Klassificering av cellbilder med hjälp av MPEG-7-inspirerade mått och support
vector machines i cellmorfologi

to Britta,
to my family

subsection

###### Contents

  1. Acknowledgments
  2. 1 Introduction
  3. 2 Background
    1. 1 Support Vector Machines
      1. 1.1 Supervised Learning
      2. 1.2 Linear Learning Machines
      3. 1.3 Maximum Margin Classifier
      4. 1.4 Optimization Theory
      5. 1.5 The Kernel Trick
      6. 1.6 Gradient Ascent
      7. 1.7 Multiclass SVM
    2. 2 Features
      1. 2.1 Scalable Color Descriptor
      2. 2.2 Color Structure Descriptor
      3. 2.3 Color Layout Descriptor
      4. 2.4 Homogeneous Texture Descriptor
      5. 2.5 Visual Texture Features
        1. 2.5.1 Neighborhood Gray-Tone Difference Matrix
        2. 2.5.2 Coarseness
        3. 2.5.3 Contrast
        4. 2.5.4 Busyness
        5. 2.5.5 Complexity
        6. 2.5.6 Texture strength
    3. 3 Fast 2D Convolution
    4. 4 Scaling data
  4. 3 Material and Methods
    1. 5 Material
    2. 6 Implementation details
      1. 6.1 Support Vector Machine
        1. 6.1.1 A Stochastic Gradient Ascent Variant
        2. 6.1.2 Multiclass SVM
      2. 6.2 Features
        1. 6.2.1 Scalable Color Descriptor
        2. 6.2.2 Color Structure Descriptor
        3. 6.2.3 Color Layout Descriptor
        4. 6.2.4 Homogeneous Texture Descriptor
        5. 6.2.5 Neighborhood Gray-Tone Difference Matrix
      3. 6.3 Convolution
      4. 6.4 Data View
  5. 4 Experimental Setup and Results
    1. 7 Experimental Setup
      1. 7.1 Performance test method
      2. 7.2 Description of the simplified problem
    2. 8 Results
      1. 8.1 Primary Problem
      2. 8.2 Simplified Problem
  6. 5 Discussion
  7. 6 Software Usage
    1. 6.A train – Train a model
    2. 6.B cellfeatures – Generate examples from the cell database
    3. 6.C jpeg_genfeature – Feature generation from images
    4. 6.D predict – Predicting a set of features
    5. 6.E extractcelltype – Extract a class of images from the cell database
    6. 6.F extractcellid – Extract given instances from the cell database
    7. 6.G extractcellinfo – Extract statistics of instances from the cell database
    8. 6.H tolibsvm – Save cell features in libSVM format

###### List of Tables

  1. 1.1 Abundance of different types of white blood cells (leukocytes) in healthy humans
  2. 3.1 Cell types classified in the data set
  3. 4.1 Cell types left in the simplified problem
  4. 4.2 SVM cell classifier results for the primary problem
  5. 4.3 Confusion Matrix for the primary problem
  6. 4.4 SVM cell classifier results for the simplified problem
  7. 4.5 Confusion matrix for the simplified problem

## Acknowledgments

First of all thanks to Doc. Christian Balkenius and Doc. Jacek Malec for
supervising my thesis. To Dr. Ferenc Belik for managing all practical details.
to Doc. B.S. Manjunath for inspiration and lending me figure 2.3. To Sebastian
Ganslandt for initial chats about support vector machines and thesis ideas. To
all others that made this work possible.

##  1 Introduction

After the introduction of MPEG-7 descriptors by the Movie Producers Expert
Group (MPEG) committeemanjunath2001ColorTextureDescriptors  it is interesting
to see how these features perform in the field of machine learning. In this
thesis a subset of them will be tested on the problem of classifying different
cell types, i.e. cell morphology, by using Support Vector Machines.

In medicine, more specifically the fields of hematology and infectious
diseases, classifying different kinds of blood cells can be used as a tool in
diagnosis — by counting certain cells’ relative frequencies and comparing to
what is normal, conclusions can be made about possible diagnosis.

Type | Approx. Abundance
---|---
neutrophil granulocytes | 70%
eosinophil granulocytes | 1-6%
basophil granulocytes | 0.01-0.3%
lymphocyte | 20-40%
monocytes | 3-8%
Table 1.1: Abundance of different types of white blood cells (leukocytes) in
healthy humans

[Neutrophil Granulocyte, segmented (class 1)] [Neutrophil Granulocyte, band
(class 6)] [Eosinophil Granulocyte (class 2)] [Basophil Granulocyte (class 3)]
[Lymphocyte (class 4)] [Monocyte (class 5)]

Figure 1.1: Some typical images of common white blood cells

Classifying cells using microscopy is used to classify infectious diseases by
determining the relative amount of cells called neutrophils compared to the
amount of cells called lymphocytes. Typical relative frequencies of the cells
are found in table 1.1. Typical images of some common cells are found in
figure 1.1.

Another method used is flow cytometry where receptors on the cells are colored
and the different types of cells are counted. Flow cytometry uses a
complicated and expensive apparatus while microscopy is very cheap.

However, microscopy is personnel intensive, many cells are hard to classify
even for human experts, often several experts are needed to be certain. To be
able to classify cells, great efforts of training are required, even more, to
sustain competence, regular frequent work is required. This competence is
impossible to sustain at small clinics or in the countryside especially in
developing countries. Instead, samples have to be sent to hematology labs.

As processing power becomes cheaper and machine learning and computer vision
algorithms grow better, machines can help less experienced personnel or give
preliminary results while waiting for definite results.

The problem this thesis try to investigate is how well these different types
of white blood cells can be classified using a Support Vector Machine and a
set of measures on the images, called features.

There has been a lot of hype about Support Vector Machines since its
introduction in the 1990’s. SVM is applied within a broad range of fields,
from bioinformaticsLengauerBioinformaticsPredictionOfHIVCoreceptorUsage  to
food engineeringDuMultiClassificationOfPizzaUsingComputerVisionAndSVM , iris
recognitionIrisRecognitionBasedOnScoreLevelFusionByUsingSVM , texture
classification and object
recognitionZhangLocalFeaturesAndKernelsForClassificationOfTextureAndObjectCategories
. It is now one of standard tools available for machine learning—A recent
search for “Support Vector Machine” (SVM) gave 6 394 articles compared to 17
893 for “Artificial Neural Network” (ANN) which has existed for much longer.
That is why my supervisor and I chose to work with SVM.

The SVM is trained with measures of the cell images, called features or
descriptors. These are values that describe the essence of an image. In this
thesis I will describe and implement a subset of the color and texture
descriptors found in the MPEG-7 standard with minor variance. I chose to work
with MPEG-7 as a guide because of the MPEG committee’s well known expertise.

The MPEG committee developed e.g. the audio compression techniques used in
MPEG-1 Layer 3 (MP3), the video compression used in e.g. DVDs (MPEG-4) and
MPEG-7. The committee consists of experts from a broad range of areas that
deal with digital information.urlMPEG7

MPEG-7 identify several descriptors which has proved useful in the Color and
Texture Core Experimentsmanjunath2001ColorTextureDescriptors  while developing
of the standard. They have proved useful in image browsing, search and
retrievalmanjunath2000ATextureDescriptorForBrowsingAndSimilarityRetrieval  as
well as in image
classificationSpyrouFuzzySVMForImageClassificationFusingMPEG7VisualDescriptors
. Color histogram based features has been successful both in image
retrievalSergyanColorHistogramFeaturesBasedImageClassificationInContentBasedImageRetrievalSystems
and image
classificationSergyanColorHistogramFeaturesBasedImageClassificationInContentBasedImageRetrievalSystems
, ChapelleSVMForHistogramBasedImageClassification ,
BarlaOldFashionedStateOfTheArtImageClassification  systems. Texture features
like Gabor Wavelet Filter Bank used in MPEG-7 has been successfully applied to
irisIrisRecognitionBasedOnScoreLevelFusionByUsingSVM  and facial
expressionBuciuICAAndGaborRepresentationForFacialExpressionRecognition
recognition.

##  2 Background

###  1 Support Vector Machines

In this section I will briefly introduce Support Vector Machines from a
theoretical perspective. Further introduction may be found in Bishop’s
book[Bishop, , chapters 6,7 and E]. If more substance is wanted I recommend
reading the whole book by Christianini and Shawe-TaylorNello . The very
thorough coverage of the topic by its original implementor Vapnik in his
bookVapnik , sometimes called the bible, was often an additional useful source
for me.

####  1.1 Supervised Learning

Supervised learning is a kind of machine learning where the machine is fed
with examples, i.e. instances of data tied to their class. The machine is told
what class an instance belongs to.

The task that a learning machine performs is to recognize an element
$\boldsymbol{\mathbf{x}}\in\mathcal{X}$ as a member of a class — to classify
it. These classes are called destination values and I use the notation
$y\in\mathcal{Y}$. In the binary case for example $\mathcal{Y}=\\{-1,+1\\}$.
The task would then be to construct a function such that
$d(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{\alpha}})=y$, given
$\boldsymbol{\mathbf{\alpha}}$ is the information the machine has previously
gathered during the training process. During training, the machine observes a
tuple of pairs

$\displaystyle
S=\big{(}(\boldsymbol{\mathbf{x}}_{1},y_{1}),\ldots,(\boldsymbol{\mathbf{x}}_{\ell},y_{\ell})\big{)}\subseteq(\mathcal{X}\times\mathcal{Y})^{\ell},$

which is called the training set, and produces parameters
$\boldsymbol{\mathbf{\alpha}}\in\mathbb{R}^{n}$ deduced from this
information.Nello

####  1.2 Linear Learning Machines

Imagine the space $\mathcal{X}$ which has $n$ dimensions. To be able to
classify instances into the two classes labeled positive, $y=+1$, or negative,
$y=-1$, a hyperplane, i.e. an affine subspace of dimension $n-1$, must be
found that separates the instances of the respective classes from each other.
If such a hyperplane exists, the data is said to be linearly separable.

Imagine a two-dimensional coordinate system in which the instances are placed.
If a straight line can be placed between the two classes of instances, the
data is linearly separable. That straight line is a hyperplane of dimension 1.
The generalized hyperplane of dimension $n-1$ is defined by the equation

$\displaystyle\langle\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{x}}\rangle+b=0.$

The normal vector $\boldsymbol{\mathbf{w}}$ is orthogonal to the hyperplane
and the bias $b$ is the hyperplane’s offset from the origin.

Now consider the function

$\displaystyle
f(\boldsymbol{\mathbf{x}})=\langle\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{x}}\rangle+b=\sum_{i=1}^{n}w_{i}x_{i}+b$
(2.1)

Where: $\boldsymbol{\mathbf{x}}$ – instance

$\boldsymbol{\mathbf{w}}$ – coefficients learned

$b$ – system bias
It will tell whether an instance is above or below the hyperplane. This is
similar to linear regression in statistics.

A decision function for the binary classification case then becomes

$\displaystyle d(\boldsymbol{\mathbf{x}})$
$\displaystyle=\mathrm{sgn}(f(\boldsymbol{\mathbf{x}}))$
$\displaystyle\mathrm{sgn}(a)$ $\displaystyle=\begin{cases}-1,&\;a<0\\\
+1,&\;a\geq 0\end{cases}$

An example of an iterative algorithm that find the vector
$\boldsymbol{\mathbf{w}}$ from a set of
$\boldsymbol{\mathbf{x}}\in\mathcal{X}$ is Rosenblatt’s perceptron which was
the first and simplest type of an Artificial Neural Networks (ANN). It is
guaranteed to converge if the data is linearly separable. This criterion could
also be written

$\displaystyle\exists\boldsymbol{\mathbf{w}}\forall i:\gamma_{i}$
$\displaystyle=y_{i}(\langle\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{x}}_{i}\rangle+b)>0,$
$\displaystyle i\in[0,\ell),$

i.e. all instances are classified correctly. The quantity $\gamma_{i}$ is
called the margin as it specifies how far from the hyperplane an instance is.
If $\boldsymbol{\mathbf{w}}$ and $b$ are normalized, to
$\frac{\boldsymbol{\mathbf{w}}}{\lVert\boldsymbol{\mathbf{w}}\rVert}$ and
$\frac{b}{\lVert\boldsymbol{\mathbf{w}}\rVert}$, then the margin is called the
geometric margin which measures the euclidean distances of the points
$\boldsymbol{\mathbf{x}}$ to the hyperplane. The closest point, the
$\boldsymbol{\mathbf{x}}_{i}$ with minimal $\gamma_{i}$, define the margin of
a hyperplane which is a stripe of empty space where no instances are. If the
data is not linearly separable $\exists i:\gamma_{i}\leq 0$.Nello , Bishop

####  1.3 Maximum Margin Classifier

The task of a maximum margin classifier is to maximize the margin which can be
motivated, using statistical learning theory, gives the least generalization
error.

The maximum margin solution, the optimal $\boldsymbol{\mathbf{w}}$ and $b$, is
found by solving

$\displaystyle\operatorname*{arg\
max}_{\boldsymbol{\mathbf{w}},b}\left\\{\min_{i}\frac{\gamma_{i}}{\lVert\boldsymbol{\mathbf{w}}\rVert}\right\\}=\operatorname*{arg\
max}_{\boldsymbol{\mathbf{w}},b}\left\\{\frac{1}{\lVert\boldsymbol{\mathbf{w}}\rVert}\min_{i}y_{i}(\langle\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{x}}_{i}\rangle+b)\right\\}$

To solve this first rescale
$\boldsymbol{\mathbf{w}}\rightarrow\kappa\boldsymbol{\mathbf{w}}$ and
$b\rightarrow\kappa b$. The distance to the hyperplane is still the same
$\min_{i}\gamma_{i}$. Then set

$\displaystyle\gamma_{j}=y_{j}(\langle\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{x}}_{j}\rangle+b=1$

for the point $\boldsymbol{\mathbf{x}}_{j}$ that is closest to the hyperplane.
All points will then have $\gamma_{i}\geq 1$ and since the minimum
$\gamma_{j}=1$ all that have to be done is to maximize
$\lVert\boldsymbol{\mathbf{w}}\rVert^{-1}$ or minimize
$\lVert\boldsymbol{\mathbf{w}}\rVert^{2}$. The problem that is left is to

$\begin{split}\text{find}&\quad\operatorname*{arg\
min}_{\boldsymbol{\mathbf{w}},b}\frac{\lVert\boldsymbol{\mathbf{w}}\rVert^{2}}{2},\\\
\text{subject to}&\quad\gamma_{i}\geq 1,\end{split}$ (2.2)

which is much easier. This problem is what is called a quadratic programming
problem and can be solved using the theory of optimization theory and Lagrange
Multipliers.Nello , Bishop

####  1.4 Optimization Theory

The theory on Lagrangian multipliers states that to

$\begin{split}\text{optimize}\quad&f(\boldsymbol{\mathbf{x}})\\\ \text{subject
to}\quad&g(\boldsymbol{\mathbf{x}})\geq 0\\\ \end{split}$

one should optimize the Lagrangian function

$\begin{split}&L(\boldsymbol{\mathbf{x}},\alpha)=f(\boldsymbol{\mathbf{x}})+\alpha
g(\boldsymbol{\mathbf{x}})\\\ \text{subject
to}\quad&g(\boldsymbol{\mathbf{x}})\geq 0\\\ &\alpha\geq 0\\\ &\alpha
g(\boldsymbol{\mathbf{x}})=0.\end{split}$

These conditions are known as the Karush-Kuhn-Tucker(KKT) conditions. More
generally, to add more constraints $g_{j}(\boldsymbol{\mathbf{x}})$, replace
the $\alpha g(\boldsymbol{\mathbf{x}})$ with a linear combination of all
Lagrange multipliers $\alpha_{j}$ and their corresponding functions
$g_{j}(\boldsymbol{\mathbf{x}})$Bishop :

$\begin{split}\text{optimize}\quad\qquad\\!&L(\boldsymbol{\mathbf{x}},\\{\alpha_{j}\\})=f(\boldsymbol{\mathbf{x}})+\sum_{j=1}^{J}\alpha_{j}g_{j}(\boldsymbol{\mathbf{x}})\\\
\text{subject to}\quad\forall j:\ &g_{j}(\boldsymbol{\mathbf{x}})\geq 0\\\
&\alpha_{j}\geq 0\\\ &\alpha_{j}g_{j}(\boldsymbol{\mathbf{x}})=0.\end{split}$

In order to quickly find a solution to (2.2) it can now be rewritten as the
Lagrangian function

$\displaystyle
L(\boldsymbol{\mathbf{w}},b,\alpha)=\underbrace{\frac{1}{2}\lVert\boldsymbol{\mathbf{w}}\rVert^{2}}_{f(\boldsymbol{\mathbf{x}})}-\sum_{i=1}^{\ell}\alpha_{i}\underbrace{(y_{i}(\langle\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{x}}_{i}\rangle+b)-1)}_{g_{i}(\boldsymbol{\mathbf{x}})}.$

The constraint function is negative because we are minimizing wrt
$\lVert\boldsymbol{\mathbf{w}}\rVert$ and $b$ while maximizing wrt
$\boldsymbol{\mathbf{\alpha}}$. To finally arrive at what is called the dual
representation of the maximum margin problem the derivatives of $L$ wrt to
$\boldsymbol{\mathbf{w}}$ and $b$, are set to $0$. Maximizing this dual
representation,

$\begin{split}&W(\boldsymbol{\mathbf{\alpha}})={\tilde{L}}(\boldsymbol{\mathbf{\alpha}})=\sum_{i=1}^{\ell}\alpha_{i}-\frac{1}{2}\sum_{i=1}^{\ell}\sum_{j=1}^{\ell}\alpha_{i}\alpha_{j}y_{i}y_{j}\langle\boldsymbol{\mathbf{x}}_{i},\boldsymbol{\mathbf{x}}_{j}\rangle,\\\
\text{by finding}\quad&\boldsymbol{\mathbf{\alpha}},\\\ \text{subject
to}\quad&\forall i:\alpha_{i}\geq 0,\\\
&\sum_{i=1}^{\ell}\alpha_{i}y_{i}=0,\end{split}$ (2.3)

will construct the maximal margin classifier.Nello , Bishop , Vapnik

The instances that have a corresponding $\alpha_{i}>0$ are called support
vectors. That is because they lie on the margin. They are thus used in the
decision function.

Note how the input variables $\boldsymbol{\mathbf{x}}_{i}$ are only used in an
inner product which let the SVM avoid the curse of dimensionality caused by a
data set with instances of too high dimension.Nello

####  1.5 The Kernel Trick

The Kernel Trick is used implicitly in Support Vector Machines but it has also
been tried out in e.g. RBF Networks, which is a kind of ANN.Bishop

The inner product used in the dual optimization problem can be a linear one.
Though it will not separate the instances fully when the dataset is not
linearly separable, data must be mapped to another space where it is.

A non-linear feature function $\phi(\boldsymbol{\mathbf{x}})$ can do such a
mapping. However, there is no need to know the feature function explicitly, it
is easier to define it implicitly via a Mercer Kernel.Nello

A complete, normed space with an inner product is called a Hilbert Space One
of the beauties of Hilbert spaces lies in that any given function in the
$L_{2}$ space could be approximated infinitely well in the
$\lVert\cdot\rVert_{2}$ and represented by an infinite linear combination of
some coefficients and some basis functions. An example of this is the Fourier
Series using Fourier coefficients and the Dirichlet Kernel Functions
$\\{e^{-ikx}\\}_{k}$.

A special kind of Hilbert spaces are the ones which are called Reproducing
Kernel Hilbert spaces. A function
$\langle\boldsymbol{\mathbf{x}}_{i},\boldsymbol{\mathbf{x}}_{j}\rangle=K(\boldsymbol{\mathbf{x}}_{i},\boldsymbol{\mathbf{x}}_{j})=\phi(\boldsymbol{\mathbf{x}}_{i})\phi(\boldsymbol{\mathbf{x}}_{j})$
is called a kernel when it satisfies the criteria in Mercer’s Theorem.

A Mercer kernel $K$ is defined as an inner product on elements of some space
$\mathcal{X}$.Nello  An inner product is a function that is a positive-
definite sesqui-linear111anti-linear in the second argument and linear in the
first form. In the $\mathbb{R}$ case this becomes a function

$\displaystyle\langle\cdot,\cdot\rangle$
$\displaystyle:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ such that
$\displaystyle K(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{z}})$
$\displaystyle=\langle\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{z}}\rangle={\langle\boldsymbol{\mathbf{z}},\boldsymbol{\mathbf{x}}\rangle}={K(\boldsymbol{\mathbf{z}},\boldsymbol{\mathbf{x}})}$
(Symmetry) $\displaystyle
K(a\boldsymbol{\mathbf{x}}+b\boldsymbol{\mathbf{y}},c\boldsymbol{\mathbf{z}})$
$\displaystyle=ab{c}\big{(}K(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{z}})+K(\boldsymbol{\mathbf{y}},\boldsymbol{\mathbf{z}})\big{)}$
(Bilinearity)
$\displaystyle\forall\boldsymbol{\mathbf{x}}:K(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{x}})$
$\displaystyle\geq 0$ (Positivity) $\displaystyle
K(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{x}})$
$\displaystyle=0\iff\boldsymbol{\mathbf{x}}=\boldsymbol{\mathbf{0}}$
(Definiteness) A Mercer kernel also have non-negative eigenvalues
$\lambda_{i}$ of the Gram matrix $\boldsymbol{\mathbf{G}}$ since it’s defined
as a Hermitian matrix $\displaystyle\forall i:\lambda_{i}$ $\displaystyle\geq
0|\boldsymbol{\mathbf{G}}$ (Positive semi-definite Gram matrix)
$\displaystyle\boldsymbol{\mathbf{G}}$
$\displaystyle=\Big{(}K\big{(}{\boldsymbol{\mathbf{x}}}_{i},{\boldsymbol{\mathbf{x}}}_{j}\big{)}\Big{)}_{i,j\in[1,\ell]^{2}}$
(2.4)

Note that the elements of the space $\mathcal{X}$ do not need to be real
vectors as they will be in this context, they could also be e.g. strings of
symbols as well. As soon as a symmetric sesqui-linear positive-definite
function could be defined on the elements of the space $\mathcal{X}$, the
space becomes an inner product space and the Support Vector Machine will do
its job.Nello

Here are some commonly used Mercer kernels defined on
$\mathbb{R}^{n}\times\mathbb{R}^{n}$Nello , Bishop , Vapnik :

$\displaystyle\langle\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}\rangle_{Linear}$
$\displaystyle={\boldsymbol{\mathbf{x}}}^{\mathsf{T}}\boldsymbol{\mathbf{y}}$
(Linear, dot product, kernel)
$\displaystyle\langle\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}\rangle_{Poly}$
$\displaystyle=\Big{(}{\boldsymbol{\mathbf{x}}}^{\mathsf{T}}\boldsymbol{\mathbf{y}}+1\Big{)}^{d}$
(Complete Polynomial of degree $d$)
$\displaystyle\langle\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}\rangle_{RBF}$
$\displaystyle=\exp\left(-\frac{1}{2\sigma^{2}}\lVert\boldsymbol{\mathbf{x}}-\boldsymbol{\mathbf{y}}\rVert\right)$
(Gaussian, Radial Basis Function)
$\displaystyle\langle\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}\rangle_{MLP}$
$\displaystyle=\tanh({\boldsymbol{\mathbf{x}}}^{\mathsf{T}}\boldsymbol{\mathbf{y}}+b)$
(Multilayer perceptron, for some $b$) the norm used in RBF is usually the
euclidean distance, $p=2$ below $\displaystyle\lVert x-y\rVert_{L^{p}}$
$\displaystyle=\Big{(}\sum_{i}\lvert x_{i}-y_{i}\rvert^{p}\Big{)}^{1/p}$
($L^{p}$ distance)

####  1.6 Gradient Ascent

An easy approach to find coefficients $\boldsymbol{\mathbf{\alpha}}$ is to
update them in the direction of the gradient of the objective function
$W(\boldsymbol{\mathbf{\alpha}})$,

$\displaystyle\frac{\partial
W(\boldsymbol{\mathbf{\alpha}})}{\partial\alpha_{i}}$
$\displaystyle=1-y_{i}\sum_{j=1}^{\ell}\alpha_{j}y_{j}\langle\boldsymbol{\mathbf{x}}_{i},\boldsymbol{\mathbf{x}}_{j}\rangle.$
To maximize the objective function $W(\boldsymbol{\mathbf{\alpha}})$ one could
just iterate $\displaystyle\alpha_{i}^{\prime}$
$\displaystyle\leftarrow\alpha_{i}+\eta\frac{\partial
W(\boldsymbol{\mathbf{\alpha}})}{\partial\alpha_{i}}.$

Where: $\eta$ – the learning rate
It is shown e.g. in Nello’s book that setting
$\eta=\frac{1}{K(\boldsymbol{\mathbf{x}}_{i},\boldsymbol{\mathbf{x}}_{j})}$
maximizes the gain if the $\alpha_{i}\in[0,C],C\in\mathbb{R}$ and that
convergence is guaranteed if the hyperplane exists.Nello

####  1.7 Multiclass SVM

There are three major methods for training a set of classifiers to be able to
classify several classesHsuLinMultiClass , i.e. $|\mathcal{Y}|=k>2$.

In the one-against-the-rest method $k$ binary classifiers are created where
classifier $i\in[0,k)$ is told that all examples with class $i$ are positive
and the rest are negative. When predicting which class
$\boldsymbol{\mathbf{x}}$ belongs to all classifiers are tested and the one
which gave the highest certainty wins.

In the one-against-one method $k(k-1)/2$ binary classifiers are created such
that all 2-combinations of classes $i,j$ have a corresponding classifier.

$\displaystyle C^{n}_{2}={n\choose
2}=\frac{n!}{2!(n-2)!}=\frac{n(n-1)(n-2)!}{2(n-2)!}=\frac{n(n-1)}{2}$

The prediction is then done by voting, all binary classifiers vote on their
respective class $i$ or $j$. The class with the highest vote wins, this
approach is called the ”Max Wins” strategy.

Direct Acyclic Graph SVM (DAGSVM) is the third method. It uses the same
training method as one-against-one but a different decision mechanism. The
classifiers are placed in a rooted DAG with the classifiers as internal nodes
and the classes as leaves. Starting at the root a binary decision means move
either left or right. When a leaf is reached the decision is
done.HsuLinMultiClass

###  2 Features

Features, or descriptors, try to take useful information out of an image —
color distribution, measures on edges and texture properties. They capture
information in a more condensed and efficient way than by just using the color
values in each pixel.

These descriptors are also scale invariant — it does not matter which size the
images have. This is necessary as the images have different sizes.

Scalable Color Descriptor, Color Structure Descriptor and Color Layout
Descriptor are the three color descriptors that I describe below and that are
implemented in the project. After the description of those come descriptions
of two texture descriptors. One of them is similar to the Homogeneous Texture
Descriptor from MPEG-7. Another set of descriptors, named Visual Texture
Features, is from an article by Amadasum and King which describe computational
measures which approximate how humans perceive
texture.Amadasun1989TexturalFeaturesCorrespondingToTexturalProperties

####  2.1 Scalable Color Descriptor

The HSV space is uniformly quantized into a 3D histogram of 256 bins. Hue is
divided into 16 levels, Saturation into 4 and Value into 4. In the MPEG-7
specification the $16\times 4\times 4=256$ bins are truncated to a 11-bit
integer mapped to a non-linear 4-bit representation and then encoded using a
Haar transform to drastically reduce space footprint. The scalability in this
descriptor comes from the ability to choose how many Haar coefficients to
store, see an article by Manjunath et al. for more
details.manjunath2001ColorTextureDescriptors

####  2.2 Color Structure Descriptor

To express local color structure in an image this descriptor slides an
$8\times 8$-structuring element across the image counting in how many of these
elements each color exists. By this technique one can differ between the
images in figure 2.1.

This descriptor is scale invariant as the structuring elements spatial extent
scale with the image size. The structure element uses replacement sub-sampling
if the image is larger than $256\times 256$ pixels. If e.g. a $512\times 512$
image is processed every other row and column will represent the image and the
rest of the $2\times 2$ areas are thrown away. More generally

$\displaystyle p$
$\displaystyle=\max\\{0,\mathrm{round}(0.5\log_{2}(WH)-8)\\}$ (2.5)
$\displaystyle K$ $\displaystyle=2^{p},\,\,E=8K$ (2.6)

Where: $E\times E$ – spatial extent of the structuring element

$K$ – sub-sampling factor
Each bin in the generated histogram represents the number of occasions a
structuring element is found to contain the color associated with the bin.

[Blackness $\pi 2^{2}$]  [Blackness $4\pi 1^{2}$]

Figure 2.1: These images contain the same amount of black and would yield an
identical color histogram but a different color structure descriptor.

####  2.3 Color Layout Descriptor

This is kind of a low-pass filter capturing spatial information. Again it is
inspired by the MPEG-7 specification. The image is first divided in $8\times
8$ blocks. Then interpolation sub-sampling222the average of all pixels
involved in the block represent the whole block as opposed to replacement sub-
sampling where a single pixel represent the whole block is applied, i.e.
calculating the average color in each block, giving one representative color
for each block. A 2D discrete cosine transform (DCT-II) is performed on the
resulting $8\times 8$ matrix. Low-frequency coefficients are selected using
zigzag scanning order, see figure 2.2. In MPEG-7 the 6 first Y, the 3 first of
U and V coefficients are extracted.

$\displaystyle\begin{pmatrix}1&3&4&10&11\\\ 2&5&9&12&19\\\ 6&8&13&18&20\\\
7&14&17&21&24\\\ 15&16&22&23&25\\\ \end{pmatrix}$ Figure 2.2: Zigzag scan
order of a $5\times 5$ matrix

####  2.4 Homogeneous Texture Descriptor

Gabor wavelets have proved to be the best set of features compared to pyramid-
structured wavelet transform (PWT), tree-structured wavelet transform (TWT)
and multi-resolution simultaneous autoregressive model (MR-SAR) based
descriptors.manjunath1996TextureFeaturesForBrowsingAndRetrievalOfImageData
They are used in the MPEG-7 Homogeneous Texture Descriptor (HTD).

Gabor wavelets are a family of modulated Gaussians, they form a complete basis
set implying that, any given function $f(\cdot,\cdot)$ can be expanded in
terms of these basis functions. However, as they are not orthonormal, there is
redundant information present in a set of coefficients. To decrease that
redundancy I follow the strategy used by Manjunath et al., that is aligning
the Gaussians such that their half-peaks meet like in figure
2.3.manjunath2000ATextureDescriptorForBrowsingAndSimilarityRetrieval

Figure 2.3: $U_{hi}=0.4,U_{lo}=0.05,S=5,K=6$

. To achieve this we first make a change of variables. The Gaussian is a
Gaussian in both frequency and space domains. The width of the Gaussian in the
frequency domain ($\sigma_{u},\sigma_{v}$) is inversely related to the
Gaussian in the space domain ($\sigma_{x},\sigma_{y}$). In other words, the
wider the Gaussian, the narrower its bandwidth.wikiGaussianFunction ,
Wallis_Linear_Models_Of_Simple_Cells_Mammal_Vision_Model

$\displaystyle\sigma_{x}$
$\displaystyle=\frac{1}{2\pi\sigma_{u}},\quad\\!\sigma_{y}=\frac{1}{2\pi\sigma_{v}}$

These parameters are needed for scaling

$\displaystyle a$
$\displaystyle=(U_{hi}/U_{lo})^{1/(S-1)},\quad\\!\sigma_{u}=\frac{(a-1)U_{hi}}{(a+1)\sqrt{2\ln
2}},$ $\displaystyle\sigma_{v}$
$\displaystyle=\tan\left(\frac{\pi}{2K}\right)\\!\left[U_{hi}-2\ln
2\left(\frac{\sigma^{2}_{u}}{U_{hi}}\right)\middle]\\!\middle[2\ln
2-\left(\frac{(2\ln
2)\sigma_{u}}{U_{hi}}\right)^{\\!\\!2\,}\right]^{\frac{1}{2}}$

Where: $U_{lo}\in\mathbb{R}$ – lower center frequency of interest

$U_{hi}\in\mathbb{R}$ – upper center frequency of interest

$m\in[0,S)\subset\mathbb{Z}^{+}$ – scale index

$S\subset\mathbb{N}$ – number of scales

$a>1\in\mathbb{R}$ – scale factor
For different orientations the image needs to be rotated before filtering and
scaling wrt $a$.

$\displaystyle x^{\prime}$ $\displaystyle=a^{-m}(x\cos\theta+y\sin\theta)$
$\displaystyle y^{\prime}$ $\displaystyle=a^{-m}(-x\sin\theta+y\cos\theta)$
$\displaystyle\theta$ $\displaystyle=n\pi/K$

Where: $n\in[0,K)\subset\mathbb{Z}^{+}$ – orientation index

$K\in\mathbb{N}$ – number of orientations

$\theta\in[0,\pi)$ – orientation angle

The generated filter bank are matrices that should be convoluted with the
image

$\displaystyle I^{\prime}=I*G$

Where: $*$ – the convolution operator
See section 3 for details about 2D convolution. In figure 3.1 images of the
Gabor wavelet filter bank kernels of different orientations are presented.

In MPEG-7, rotation invariance is achieved in this descriptor, by rotating the
features in the direction of the dominant direction.

####  2.5 Visual Texture Features

The features described in the article by Amadasun and King are implemented.
These are features corresponding to properties of texture that humans can
perceive. In the article measures of coarseness, contrast, busyness,
complexity and strength are introduced and compared by rank with how humans
sensed ten natural textures from the widely used Brodatz’s album. I give here
a very brief overview of the proposed measures. They all use a column vector
called neighborhood gray-tone difference matrix
(NGTDM).Amadasun1989TexturalFeaturesCorrespondingToTexturalProperties

#####  2.5.1 Neighborhood Gray-Tone Difference Matrix

In a pixel $p$ with coordinates $\langle k,l\rangle$ neighborhood of size $d$,
i.e. of the square surrounding a pixel, but without the center pixel the mean
is calculated.

$\displaystyle\begin{split}\bar{A}_{p}=\bar{A}(k,l)=\frac{1}{W-1}\left[\sum_{m=-d}^{d}\sum_{n=-d}^{d}f(k+m,l+n)\right],\\\
\quad(m,n)\neq(0,0)\end{split}$ (2.7)

Where: $W=(2d+1)²$
The $i$th entry in the NGTDM is a sum of deviations from the mean of the
center pixel, only concerning those pixels in the image which do not lie in
the peripheral regions of width $d$.

$\displaystyle s(i)=\begin{cases}\displaystyle\sum_{p\in
N_{i}}\left|i-\bar{A}_{p}\right|,&\;\text{there is a pixel with gray-tone
}i\\\ 0,&\;\text{otherwise}\end{cases}$ (2.8)

Where: $N_{i}$ – the pixels with gray-tone $i$

$G_{h}$ – the largest gray-tone
The relative frequency, i.e. the probability of occurrence, of different gray-
tones is calculated as:
Amadasun1989TexturalFeaturesCorrespondingToTexturalProperties

$\displaystyle p_{i}$ $\displaystyle=\lvert N_{i}\rvert/n,$ $\displaystyle n$
$\displaystyle=(width-2d)(height-2d).$ (2.9)

Note that (2.9) allows a rectangular region of interest as opposed to the
square regions used in the article by Amadasun and King, and that $n$ replaces
$n^{2}$ in the formulas.

#####  2.5.2 Coarseness

Coarseness is a measure of how rough a surface is, e.g. how large particles it
is composed of.

$\displaystyle f_{cos}=\left[\epsilon+\sum_{i=0}^{G_{h}}p_{i}s(i)\right]^{-1}$

This is (inversely) a weighted sum of the deviations from the center pixels
wrt the surrounding pixels. The small value $\epsilon$ is to cope with
division by $0$.

#####  2.5.3 Contrast

High contrast means the intensity difference between neighboring regions is
large.

$\displaystyle f_{con}$
$\displaystyle=\left[\frac{1}{N_{g}(N_{g}-1)}\sum_{i=0}^{G_{h}}\sum_{j=0}^{G_{h}}p_{i}p_{j}(i-j)²\middle]\middle[\frac{1}{n}\sum_{i=0}^{G_{h}}s(i)\right]$
$\displaystyle N_{g}$ $\displaystyle=\sum_{i=0}^{G_{h}}Q_{i}$ $\displaystyle
Q_{i}$ $\displaystyle=\begin{cases}1,&\;\text{if }p_{i}\neq 0\\\
0,&\;\text{otherwise}\end{cases}$

Where: $N_{g}$ – the number of different gray-tones present in the image
The first factor is used to reflect the dynamic range of gray scale weighted
with the product of relative frequencies of the two gray-tone values under
consideration. The second factor increases with the amount of local variation
in intensity.

#####  2.5.4 Busyness

A busy texture is one where the spatial frequency of intensity changes are
high.

$\begin{split}f_{bus}=\left.\sum_{i=0}^{G_{h}}p_{i}s(i)\middle/\sum_{i=0}^{G_{h}}\sum_{j=i}^{G_{h}}ip_{i}-jp_{j}\right.\\!\\!,\\\
p_{i}\neq 0,\>p_{j}\neq 0\end{split}$

The numerator is a measure of the spatial rate of change in intensity,
inversely related to coarseness. The denominator is a summation of the
magnitude of differences between the different gray-tone values. This formula
differs slightly from the one described in the article by Amadasun and
KingAmadasun1989TexturalFeaturesCorrespondingToTexturalProperties  — I’m
certain there’s a typo in that formula making it always zero.

#####  2.5.5 Complexity

Complexity means high information content. This could mean many primitives or
patches, especially if they have different average intensity.

$\displaystyle f_{com}=\sum_{i=0}^{G_{h}}\sum_{j=0}^{G_{h}}\left.\frac{\lvert
i-j\rvert}{n(p_{i}+p_{j})}\Big{(}p_{i}s(i)+p_{j}s(j)\Big{)}\right.$

An elaborate description of this formula (and the others in this section) are
found in the article by Amadasun and
KingAmadasun1989TexturalFeaturesCorrespondingToTexturalProperties .

#####  2.5.6 Texture strength

A strong texture is generally referred to as strong if its building blocks are
easily definable and clearly visible. Such texture tend to look attractive.
However a strong texture is difficult to define
conciselyAmadasun1989TexturalFeaturesCorrespondingToTexturalProperties . It is
defined as

$\displaystyle
f_{str}=\frac{\displaystyle\sum_{\begin{subarray}{c}i=0\end{subarray}}^{G_{h}}\sum_{\begin{subarray}{c}j=0\end{subarray}}^{G_{h}}\left(p_{i}+p_{j}\middle)\middle(i-j\right)^{2}}{\displaystyle\epsilon+\sum_{i=0}^{G_{h}}s(i)}\\!,\qquad{p_{i}\neq
0,\>p_{j}\neq 0}.$

Where the numerator is a factor stressing the differences between intensity
levels, and therefore may reflect intensity differences between adjacent
primitives. The probabilities $p_{.}$ tend to be high for large primitives.
The denominator would be small for coarse texture and high for busy or fine
textures considering the definition in (2.8).

###  3 Fast 2D Convolution

Two-dimensional discrete convolution in the spatial domain is defined as

$\displaystyle(f*g)[n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{m=-\infty}^{\infty}f[m]\cdot
g[n-m].$

By the Circular Convolution TheoremwikiCircularConvolution  this can instead
be done in the frequency domain considering

$\displaystyle\mathcal{F}\\{f*g\\}=\mathcal{F}\\{f\\}\cdot\mathcal{F}\\{g\\}$
(2.10)

Where: $*$ – the convolution operator

$\mathcal{F}\\{\cdot\\}$ – the Fourier Transform (FT)
First apply FT to image and to convolution kernel, then multiply the two
matrices element-wise. To get the filtered image just apply inverse FT.

For this to work the kernel has to be placed in a matrix the same size as the
image, wrapped around the origin333origin aka DC component, zero frequency,
which in FFTW is at position $\langle 0,0\rangle$, like in figure 2.4. Also,
there are border cases in the image, it has to be padded with wraparound
pixels.Convolution2DNVIDIA

$\displaystyle\overbrace{\begin{pmatrix}11&12&13&14&15\\\ 21&22&23&24&25\\\
31&32&33&34&35\\\ 41&42&43&44&45\\\ 51&52&53&54&55\\\
\end{pmatrix}}^{\text{Example $5\times 5$
kernel}}\xrightarrow{\text{Layout}}\overbrace{\begin{pmatrix}33&34&35&0&\dots&0&31&32\\\
43&44&45&0&\dots&0&41&42\\\ 53&54&55&0&\dots&0&51&52\\\ 0&0&0&0&\dots&0&0&0\\\
[5]{8}\\\ 0&0&0&0&\dots&0&0&0\\\ 13&14&15&0&\dots&0&11&12\\\
23&24&25&0&\dots&0&21&22\end{pmatrix}}^{\text{Image height $\times$ Image
width matrix}}$ Figure 2.4: How to make sure the kernel wraps around the
origin in frequency space

###  4 Scaling data

Scaling is very important. If scaling is not applied to all features a feature
with a larger numeric range may dominate others with smaller numeric range.

$\displaystyle\displaystyle range$ $\displaystyle=\max_{i}x_{i}-\min_{i}x_{i}$
$\displaystyle midrange$
$\displaystyle=\left(\displaystyle\max_{i}x_{i}+\min_{i}x_{i}\middle)\right/\\!2$
$\displaystyle x^{\prime}_{i}$
$\displaystyle=\begin{cases}\displaystyle\frac{x_{i}-midrange}{range/2},&range\neq
0\\\ 0,&range=0\end{cases}$

Where: $x_{i}$ – feature value of example $i$

$i\in[0,\ell)\subset\mathbb{Z}^{+}$

##  3 Material and Methods

###  5 Material

Blood samples were taken from four individuals. The cells were photographed on
a CellaVision™ DM-96. The width of the images lies in the range $[119,267]$.
The height of the images lies in the range $[119,258]$. On average an image is
about $139\times 139$ pixels. This correspond to about $13.7$
$\umu\mathrm{m}$.

The cells are normal, e.g. there are no cancer cells or malaria infected
cells. There are very few (2) blast cells indicating the only possible cancer
type would be lymphoma, i.e. a cancer in the lymph nodes.

The cells were classified on the CellaVision™ DM-96 and its result was taken
as ground truth. The machine is 90% to 95% correct depending on the
individual. The cell types of the data set are given in table 3.1. Typical
relative frequencies of the cells are found in table 1.1. Typical images of
some common cells are found in figure 1.1.

Class No. | Class Name
---|---
1 | neutrophil granulocytes, segmented
6 | neutrophil granulocytes, band
2 | eosinophil granulocytes
3 | basophil granulocytes
4 | lymphocytes
7 | lymphocytes, variants
5 | monocytes
9 | myelocytes
10 | meta-myelocytes
11 | blast, immature cell
21 | artifacts
24 | broken cell
25 | thrombocytes (platelets)
29 | clots of thrombocytes
Table 3.1: Cell types classified in the data set

From the set of images of the cells a range of descriptors, or features, were
extracted. A set of features extracted from a single image, called instance or
example, is denoted $\boldsymbol{\mathbf{x}}$ and the space of all possible
features is denoted $\mathcal{X}$.

A Support Vector Machine (SVM) was trained using the set of features
described.

###  6 Implementation details

####  6.1 Support Vector Machine

The SVM was written in C++ within the Boost C++ Libraries framework. The Gram
matrix $\boldsymbol{\mathbf{G}}$, defined in (2.4), the output of the kernel
function, is cached in memory to dramatically reduce running time.

#####  6.1.1 A Stochastic Gradient Ascent Variant

Stochastic gradient ascent differs from ordinary gradient ascent in that the
coefficients $\alpha_{i}$ updated are used right away, instead of in the next
iteration. In this project a variant of the stochastic gradient ascent method
of training a SVM were implemented.

The coefficients $\boldsymbol{\mathbf{\alpha}}_{KKT}$ that invalidate the
Karush-Kuhn-Tucker (KKT) conditions are selected first for update. They are
likely the ones that will affect the solution most rapid. When these satisfies
the KKT conditions, or when no progress has been made in some iterations, the
greater problem of updating all coefficients $\boldsymbol{\mathbf{\alpha}}$ is
considered.

#####  6.1.2 Multiclass SVM

I use the one-against-the-rest methodHsuLinMultiClass  because it is the
simplest and it has similar precision to the latter twoVapnik ,
HsuLinMultiClass . The latter two are however faster to train because they can
train all the classifiers at once.Nello

####  6.2 Features

#####  6.2.1 Scalable Color Descriptor

In MPEG-7 the 3D color histogram bins are reduced in size by truncation and
encoding (see 2.1). To release the SVM from this hassle it receives the values
as ordinary real values representing the relative frequency of color channel
values. The bounded time complexity to calculate this descriptor is
$O(3W\\!H)$.

#####  6.2.2 Color Structure Descriptor

This is implemented by calculating a histogram for each structuring element
and then summing over all structuring elements

$\displaystyle
h(m)=\sum_{i=1}^{\frac{W-8K}{K}}\sum_{j=1}^{\frac{H-8K}{K}}\min\\{1,h_{s_{i,j}}(c_{m})\\}$
(3.1)

Where: $m$ – bin index in the final histogram

$c_{m}$ – quantized color level

$h_{s_{i,j}}$ – histogram for structuring element $\langle i,j\rangle$
Calculating this descriptor is much more expensive than Scalable Color
Descriptor described in section 2.1, $O(\frac{(w-8k)(h-8k)}{k}8^{2})$ for each
channel, this is more than a 30-fold increase on a $640\times 480$ image
compared to the above.

#####  6.2.3 Color Layout Descriptor

The Discrete Cosine Transform of type DCT-II is calculated using the software
library FFTW3 (Fastest Fourier Transform in the West). The zigzag scanning
order described in figure 2.2 is implemented as an C++ STL iterator using the
simple algorithm presented in listing 1. A wider low pass band is used than in
MPEG-7. The 10 first Y (6 in MPEG-7), the 5 first of U and V (3) coefficients
are extracted.

⬇

x = 0; y = 0; forward = true;

value_type get_current() { return source(x,y); }

void next() {

if (forward)

if (y < length-1) {

y ++; x –;

if (x < 0) {

x = 0;

forward = false;

}

} else

if (y == length-1) {

x ++;

forward = false;

}

else

if (x < length-1) {

x ++; y –;

if (y < 0) {

y = 0;

forward = true;

}

} else

if (x == length-1) {

y ++;

forward = true;

}

}

Listing 1: Simplified source for the implemented zigzag order on a
length$\times$length square matrix

#####  6.2.4 Homogeneous Texture Descriptor

By symmetry the filter might as well be rotated instead of the image and since
that is more efficient that is what is done. The bandwidth $b$ is set to 1
octave by relation (3.2) and setting $\sigma=\sigma_{x}$

$\displaystyle\frac{\sigma}{\lambda}$
$\displaystyle=\frac{1}{\pi}\sqrt{\frac{\ln
2}{2}}\frac{2^{b}+1}{2^{b}-1}\approx 0.5622$ (3.2)

In MPEG-7 rotation invariance in this descriptor is achieved by rotating the
features in the direction of the dominant direction. This is not implemented
in this project.

In figure 3.1 images of the Gabor wavelet filter bank kernels of different
orientations are presented.

[$\theta=0°$]  [$\theta=36°$]  [$\theta=72°$]

[$\theta=108°$]  [$\theta=144°$]  [$\theta=180°$]

Figure 3.1: Gabor Filter bank at scale = $S-1$ at different orientations. Gray
areas are the ones with zero magnitude, darker is negative, lighter is
positive

#####  6.2.5 Neighborhood Gray-Tone Difference Matrix

The $\bar{A}$ used in the Neighborhood Gray-Tone Difference Matrix (2.7) can
be divided into subproblems which do not need to be calculated every time. By
keeping the center value $(m,n)=(0,0)$ in the sum (not writing out
normalization)

$\displaystyle A^{\prime}(k,l)$
$\displaystyle=\sum_{m=-d}^{d}\sum_{n=-d}^{d}f(k+m,l+n),$ it can also be
written as $\displaystyle A^{\prime}(k,l)$
$\displaystyle=\begin{cases}\begin{split}\underbrace{A^{\prime}(k,l-1)}_{\text{above}}+\phantom{\qquad\text{or
as }}\\\ \sum_{m=-d}^{d}f(k+m,l+d)-f(k+m,l-d-1)\qquad\text{or as
}\end{split}\\\ \begin{split}\underbrace{A^{\prime}(k-1,l)}_{\text{to the
left}}+\\\ \sum_{n=-d}^{d}f(k+d,l+n)-f(k-d-1,l+n).\end{split}\end{cases}$

Given the value above or the value to the left the others can be calculated
faster.

To find all $\bar{A}$ first fill in a table with all $A^{\prime}$, from left
to right, top-down. Then for all positions remove the center value and make
sure the accumulated value is correctly normalized. The time complexity is
thereby reduced from $O(d^{2})$ per pixel to $O(d)$ per pixel.

####  6.3 Convolution

Using the method for convolution described in section 3 is much more efficient
than the naive approach of doing the calculations in the spatial domain. It
reduces the complexity from $O(K²)$ per pixel, where $K$ is the size of the
convolution kernel, to $O(\log N)$, where the image is $N\times N$ in size and
$N=2^{k},k\in\mathbb{Z}^{+}$. The last requirement make sure that the much
more efficient Fast Fourier Transform (FFT) can be used instead of a normal
Discrete Fourier Transform (DFT).

With the largest kernel used, $K²=91²=8281$, and a $1000\times 1000$ image,
$\log 1000\approx 6.9$, a thousandfold speed-up can be achieved.

These figures are however for FFT on matrices of size $N=2^{k}$. Padding to
the next larger 2-power is not implemented since the software library used for
FFT, called FFTW444Heavily used library with an impressing architecture, used
in e.g. Matlab (Fastest Fourier Transform in the West) supports other sizes
too and still provides great speed.

####  6.4 Data View

The classifiers view data. Rather than giving them the data structure holding
data directly an abstraction was built named DataView. The abstraction was
realized in 11 classes which are found together with their base abstract class
in figure 3.2. The derived classes can all be used transparently releasing the
classifier and data set loader from the tasks of the views.

Figure 3.2: Abstract class (interface) to data views and their realizations

These three views below contain pointers to the real data.

  * DataSetView

view of data represented by a DataSet instance

  * ExampleView

view of data represented by a vector of Example instances

  * ArrayView

view of data from an boost::Array, convenient for the unit tests concerning
views

The views below contain other views and just map their values. They are often
chained together to get the wanted view.

  * DataViewScaled

view the features as if they were in the range $[-1,1]$, avoids feature-wise
bias, see section 4

  * DataViewRange

selected only a subset of the examples, used in e.g. cross-validation

  * DataViewConcat

view two views as if they were one, also used in cross-validation

  * DataViewShuffle

shuffle the order of examples. It is of course not wanted to split an ordered
set and train on the first part and test on the other, a class may then be
present only in the latter

  * DataViewClassMapLinear

if e.g. only classes $\\{0,3,42,\ldots\\}$ exists it is convenient if they can
be represented by $\\{0,1,\ldots\\}$

  * DataViewClassMapBinary

one class given is said to be positive, all other is said to be negative. Used
in multiclass classifier

  * DataViewClassJoin

join groups of classes into new classes

  * DataViewClassRemove

view with a class removed

##  4 Experimental Setup and Results

###  7 Experimental Setup

The CellaVision™ DM-96 machine achieves an error rate of approx. 5-10%
depending on individual. Thus there are errors in the ground truth.

I have divided the problem in two parts.

  * •

the primary problem — the SVM should classify all classes present in the data
set.

  * •

the simplified problem — some classes are merged and others are removed.

####  7.1 Performance test method

In both cases 2-fold cross-validation is used to test performance. This means
that two models will be trained. In the first, half the data set is the
training set and the other half is the test set. In the other, the roles of
the subsets are swapped. This way both halves will act as both training and
test sets.

####  7.2 Description of the simplified problem

Class 1 and 6, Neutrophil granulocytes, segmented and band variants are merged
to form class 30. Even human experts have approx. 25% error rate on these. It
is often more a matter of opinion than of objective decision.

Class 4 and 7, Lymphocytes and their variants, are joined. The variants are
rather uncommon, there are only 8 instances in the dataset, compared to 160 of
Lymphocytes. Due to the skew distribution these are merged to form class 31.

The following classes are removed. Class 0 are unidentified objects, it is a
very heterogeneous group but there are only 6 of them. Class 21 are artifacts,
random garbage, they are removed. Class 24 are broken cells, there are only 7
of them. Class 25 and 29 are thrombocytes and clots of them, i.e. platelets.
Since they aren’t even white blood cells they are removed. Class 11, called
blast is a kind of immature cell which would be interesting to classify but
there are only two of them so they are removed as well. Class 9 and 10 are
myelocytes and meta-myelocytes, which are a development stage of different
granulocytes. There can be e.g. eosinophilic myelocytes and basophilic
myelocytes. In the dataset they are also too rare to train a general
classifier. There are only a total of 4 myelocytes in this heterogeneous
group. All classes that are left are presented again in table 4.1.

Class No. | Class Name
---|---
30 (1+6) | neutrophil granulocytes
2 | eosinophil granulocytes
3 | basophil granulocytes
31 (4+7) | lymphocytes and variants
5 | monocytes
Table 4.1: Cell types left in the simplified problem

###  8 Results

####  8.1 Primary Problem

The error rate in the primary problem is 10.8%. The type of kernel function
that was the most successful was the Polynomial kernel. This is compared to
the slightly better result using libSVM, 9.6%. See table 4.2.

Most confusion occurs between classes 1 (segmented neutrophil granulocytes)
and 6 (band neutrophil granulocytes). Much confusion is also present when
recognizing class 3 (basophil granulocytes) — they are often (2 of their total
of 7) misclassified as class 1 (segmented neutrophil granulocytes), which is a
very large group.

####  8.2 Simplified Problem

In the simplified problem the error rate is 3.1%. Also in this problem the
most successful kernel was the Polynomial kernel. This is compared to the
better result using libSVM, 2.3%. See table 4.4.

In the simplified problem most confusion (by number) occurs between class 5
(monocytes) and the new class 31 (lymphocytes). By percentage the largest
confusion occurs between class 30 (segmented and band neutrophil granulocytes)
and class 3 (basophil granulocytes). Class 3 have only 8 instances of which 3
were misclassified as 30.

Implemented SVM Results
---
Kernel Type |  Error Rate (%) | Parameters
| Total | Max | Min |
RBF with $L^{2}$ norm | 11,5385 | 12,0192 | 11,0577 | $\sigma²=20$
RBF with $L^{2}$ norm | 11,5385 | 12,0192 | 11,0577 | $\sigma²=22$
Polynomial | 11,7788 | 12,5 | 11,0577 | $d=2$
Polynomial | 11,0577 | 11,5385 | 10,5769 | $d=3$
Polynomial | 11,5385 | 12,0192 | 11,0577 | $d=4$
Polynomial | 10,8173 | 11,0577 | 10,5769 | $d=5$
Polynomial | 11,2981 | 12,0192 | 10,5769 | $d=6$
libSVM Results
RBF | 9,5923 |  | $C=512,\gamma^{-1}=8192$
Table 4.2: SVM cell classifier results for the primary problem Number of
Confusions
---
| Guessed Class |
Class | $(n)$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 11 | 21 | 24
0 | (4) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
1 | (205) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 1 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
2 | (14) | 1 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
3 | (7) | $\cdot$ | 2 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
4 | (104) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 2 | $\cdot$ | 1 | $\cdot$
5 | (32) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 2 | $\cdot$ | $\cdot$ | 1 | $\cdot$ | $\cdot$ | $\cdot$
6 | (12) | 1 | 6 | $\cdot$ | $\cdot$ | $\cdot$ | 1 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
7 | (6) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 1 | 1 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
11 | (1) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 1 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
21 | (31) | $\cdot$ | $\cdot$ | 1 | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
24 | (1) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 1 | $\cdot$
Table 4.3: Confusion Matrix for the primary problem Implemented SVM Results
---
Kernel Type |  Error Rate (%) | Parameters
| Total | Max | Min
RBF | 3,64109 | 4,42708 | 2,86458 | $C=128,\sigma^{2}=16$
RBF | 3,25098 | 4,16667 | 2,34375 | $C=512,\sigma^{2}=128$
Polynomial | 3,12094 | 3,38542 | 2,864 | $d=3$
Polynomial | 3,51105 | 3,90625 | 3,125 | $d=5$
libSVM results
RBF | 2,470 |  | $C=8,\gamma^{-1}=128$
Polynomial | 2,3407 |  | $C=8,\gamma^{-1}=128,d=3$
Polynomial | 3,5111 |  | $C=8,\gamma^{-1}=128,d=5$
Table 4.4: SVM cell classifier results for the simplified problem Number of
Confusions
---
| Guessed class |
Class | $(n)$ | 2 | 3 | 5 | 30 | 31
2 | (20) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$
3 | (8) | $\cdot$ | $\cdot$ | $\cdot$ | 3 | $\cdot$
5 | (56) | $\cdot$ | $\cdot$ | $\cdot$ | $\cdot$ | 2
30 | (517) | $\cdot$ | $\cdot$ | 1 | $\cdot$ | $\cdot$
31 | (168) | $\cdot$ | $\cdot$ | 7 | $\cdot$ | $\cdot$
Table 4.5: Confusion matrix for the simplified problem

##  5 Discussion

The accuracy achieved in the primary problem was 89.2% and in the simplified
problem 96.9%. I regard these results as good when compared to CellaVision™
DM-96’s result of the primary problem, 90-95%. One has to consider that there
are errors in the ground truth misleading the SVM. Thus, it is uncertain
whether the results are better than the DM-96 or worse. Because the DM-96 has
an error rate of about 5-10% a 0% error rate in the primary problem would mean
something like 5-10% error, while a 5% error could possibly mean 0-15% error.

I conclude that using the combination of MPEG-7 descriptors and visual texture
features in combination with SVM to classify cells is good but need further
investigation to find out how good. A more comprehensive study could
investigate whether a set of SVM or ANN variants perform better on the set of
features implemented or on the set of features developed at CellaVision™.

I would like to stress that using a SVM instead of an Artificial Neural
Network as in the CellaVision™ DM-96 machine is more statistically rigor —
Confidence intervals of the classifier can be found, which to my knowledge is
impossible in ANNs. In medicine it is important to know the strength of the
method used.

It would be very interesting to test the features on the real training set
they have developed at CellaVision™. The company has a training set of
thousands of cells classified by field experts. Some cell images required five
experts to be certain of the cell type. Without the errors in the ground truth
the results could possibly compete with the CellaVision™machine.

The result of the primary problem states that the most confused instances are
those that are guessed to be a segmented neutrophil (class 1) but that are a
band neutrophil (class 6) in the ground truth. These two often look very
similar, humans often have different opinions about which class cells are.
Also the CellaVision™ DM-96 have problems with these classes indicating that
there are several errors in the ground truth. The errors in ground truth
probably mislead the SVM. There are only 12 cell images of class 6 of which
some have the wrong class and there are 205 images of class 1, of which not
all are truly class 1. This situation pushes bias to the larger class.

In the result of the simplified problem most confusion occur between the
monocytes (class 5) and the lymphocytes (class 31). This was expected as they
are very hard to classify for both humans and the CellaVision™ DM-96. Late in
writing this thesis I discovered that there is great discrepancy in size of
these two types of cells. The discrepancy indicates that the size of the cells
could be used as a feature too.

#### Practical Use

Even the simplified problem would give useful information when applied to
medicine. Standard measures used in diagnosis involve counting the total
number of white blood cells, leukocytes, determining the distribution of
lymphocytes and granulocytes and determining the number of monocytes.

Malaria infected and cancer cells look different compared to healthy blood
cells. It would be interesting to test the features on these kind of cells to
be able to classify them as well.

#### Runtime Performance

To increase cache performance in the Color Structure Descriptor (section 2.2)
it would be wise to first extract all sub-samples i.e. the representative
color for each $K\times K$ area as the other pixels aren’t used. They will
otherwise quickly fill up the cache during memory pre-fetch. Now, the sub-
samples are viewed using a sub sampling view present in Generic Image Library
(boost::gil), contributed to Boost by Adobe. The views in GIL are virtual,
meaning they only keep information about offset calculations — no data is
duplicated.

The 2D convolution was first done in the spatial domain but I soon realized it
was way to slow with my bigger Gabor filter kernels of which the largest are
$91²$ pixels big. Instead the calculations are done in the frequency domain
which is much faster, see sections 3 and 6.3.

To improve performance of the Gabor Wavelet Filter further the kernels should
of course be kept in memory when generating features of many images, however
they are not.

To improve SVM training performance the Gradient Ascent training algorithm
must be replaced or at least improved. The algorithm implemented divide the
problem into a subproblem where the coefficients violating the KKT conditions
are first optimized. This is a heuristic called chunking in the
literatureNello . By using this, fewer elements of the Gram matrix, and their
corresponding support vectors, need to be kept in memory. This is something I
don’t take advantage of because I had enough memory for my purposes. By
refining chunking into decomposition where a fixed size chunk is optimized,
more data points can be used and convergence speed is increased. The
Sequential Minimization Optimization (SMO) takes decomposition to the extreme
and optimizes only two coefficients at a time and can thereby make sure that
the KKT condition, $\sum_{i=1}^{\ell}\alpha_{i}y_{i}=0$, is always true.
LibSVM uses a variant of this approach and it offer great
performance.CC01aLIBSVM ,
LIBSVMDimensionalityReductionViaSparseSupportVectorMachines

#### Beyond Gabor Filters

If modeling human brains is the objective, considering other approaches than
the Gabor wavelet would be interesting. A type of neurons in the first visual
cortex, called simple cells, have been recorded from monkey and cat. The
recordings and the elaborate analytical discussion in an article by Wallis
show that both difference of Gaussian$\times$Gaussian (DoGG) and Cauchy
functions model cortical cells better than Gabor wavelets for the measured
parameters.Wallis_Linear_Models_Of_Simple_Cells_Mammal_Vision_Model  In an
article by Ashour et al. three other types of transforms are suggested —
ridgelets, curvelets and
contourlets.Ashour2008SupervisedTextureClassificationUsingSeveralFeaturesExtractionTechniquesBasedOnAnnAndSvm
Perhaps they can show increased performance.

## References

  * [1] M. Amadasun and R. King.  ”Textural features corresponding to textural properties”.  Systems, Man and Cybernetics, IEEE Transactions on, 19(5):1264–1274, 1989.
  * [2] M.W. Ashour, M.F. Hussin, and K.M. Mahar.  ”Supervised texture classification using several features extraction techniques based on ANN and SVM”.  2008 IEEE/ACS International Conference on Computer Systems and Applications, pages 567–574, 2008.
  * [3] A. Barla, F. Odone, and A. Verri.  ”Old fashioned state-of-the-art image classification”.  Image Analysis and Processing, 2003.Proceedings. 12th International Conference on, pages 566–571, 2003.
  * [4] C. M. Bishop.  ”Pattern Recognition and Machine Learning”.  Springer, 2006.
  * [5] I. Buciu, C. Kotropoulos, and I. Pitas.  ”ICA and Gabor representation for facial expression recognition”.  Image Processing, 2003. ICIP 2003. Proceedings. 2003 International Conference on, 3:II–855–8 vol.3, 2003.
  * [6] C.-C. Chang and C.-J. Lin.  ”LIBSVM: a library for support vector machines”, 2001.
  * [7] O. Chapelle, P. Haffner, and V.N. Vapnik.  ”Support Vector Machines for Histogram-Based Image Classification”.  Neural Networks, IEEE Transactions on, 10(5):1055–1064, 1999.
  * [8] N. Christianini and J. Shawe-Taylor.  ”An Introduction to Support Vector Machines and other kernel-based learning methods”.  Cambridge, 2000.
  * [9] C.J. Du and D.W. Sun.  ”Multi-classification of pizza using computer vision and support vector machine”.  Journal of Food Engineering, 86(2):234–242, 2008.
  * [10] C.-W. Hsu and C.-J. Lin.  ”A Comparison of Methods for Multi-class Support Vector Machines”.  2002\.
  * [11] T. Lengauer, O. Sander, S. Sierra, A. Thielen, and R. Kaiser.  ”Bioinformatics prediction of HIV coreceptor usage”.  Nature biotechnology, 25(12):1407–1408, 2007.
  * [12] B.S. Manjunath and W.Y. Ma.  ”Texture features for browsing and retrieval of image data”.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, 18(8):837–842, 1996.
  * [13] B.S. Manjunath, J.-R. Ohm, V.V. Vasudevan, and A. Yamada.  ”Color and Texture Descriptors”.  Circuits and Systems for Video Technology, IEEE Transactions on, 11(6):703–715, 2001.
  * [14] J. M. Martínez.  ”MPEG-7 Overview”.
http://www.chiariglione.org/mpeg/standards/mpeg-7/mpeg-7.htm, 2004.  [Online;
accessed 26-August-2008].

  * [15] H.A. Park and K.R. Park.  ”Iris recognition based on score level fusion by using SVM”.  Pattern Recognition Letters, 28(15):2019–2028, 2007.
  * [16] V. Podlozhnyuk.  ”FFT-based 2D Convolution”.
http://developer.download.nvidia.com/compute/cuda/1_1/Website/projects/
/convolutionFFT2D/doc/convolutionFFT2D.pdf, June 2007.  [Online; accessed
21-July-2008].

  * [17] P.-H. Chen R.-E. Fan and C.-J. Lin.  ”Dimensionality Reduction via Sparse Support Vector Machines”.  Journal of Machine Learning Research, 6:1889–1918, 2005.
  * [18] S. Sergyan.  ”Color histogram features based image classification in content-based image retrieval systems”.  2008 6th International Symposium on Applied Machine Intelligence and Informatics, pages 221–224, 2008.
  * [19] E. Spyrou, G. Stamou, Y. Avrithis, and S. Kollias.  ”Fuzzy support vector machines for image classification fusing MPEG-7 visual descriptors”.  Integration of Knowledge, Semantics and Digital Media Technology, 2005. EWIMT 2005. The 2nd European Workshop on the (Ref. No. 2005/11099), pages 23–30, 2005.
  * [20] V.N. Vapnik.  Statistical Learning Theory.  Wiley, 1998.
  * [21] G. Wallis.  ”Linear models of simple cells: Correspondence to real cell responses and space spanning properties”.  Spatial vision, 14(3):237–260, 2001.
  * [22] Wikipedia.  ”Discrete Fourier Transform — Wikipedia, The Free Encyclopedia”.
http://en.wikipedia.org/w/index.php?title=Discrete_Fourier_transform&oldid=226961702
#Circular_convolution_theorem_and_cross-correlation_theorem, 2008.  [Online;
accessed 21-July-2008].

  * [23] Wikipedia.  ”Gaussian function — Wikipedia, The Free Encyclopedia”.
http://en.wikipedia.org/w/index.php?title=Gaussian_function&oldid=219261026,
2008.  [Online; accessed 24-June-2008].

  * [24] P. Wu, B.S. Manjunath, S. Newsam, and H.D. Shin.  ”A texture descriptor for browsing and similarity retrieval”.  Signal Processing: Image Communication, 16(1-2):33–43, 2000.
  * [25] J. Zhang, M. Marszalek, and S. Lazebnik.  ”Local Features and Kernels for Classification of Texture and Object Categories: A Comprehensive Study”.  International Journal of Computer Vision, 73(2):213–239, 2007.

##  6 Software Usage

The software produced in this project can be found at

  * •

`http://tobbe.nu/pub/2008/cell.morph.mpeg7.svm/`

The software has only been tested on an Ubuntu Linux system. However, the
software is written in portable C99 C++ and should work on all *nix platforms
that can supply the dependencies, perhaps even under cygwin under MS Windows.
The dependencies are

  * •

C99 compliant C++ compiler (GNU g++ tested)

  * •

Boost C++ Libraries, http://www.boost.org/

  * •

FFTW3 (Fastest Fourier Transform in the West 3), http://www.fftw.org/

  * •

GSL (GNU Scientific Library), http://www.gnu.org/software/gsl/

  * •

libjpeg

  * •

libpng

Below is a brief overview on how to use the most important programs in the
software package. There are other programs in the package but they are mostly
related to testing.

##  Appendix 6.A train – Train a model

This is the program where most processing is done. It can

  * •

train a model from a dataset

  * •

test a model with a dataset

  * •

load and/or save a model from/to a file

  * •

perform cross-validation

Here is the syntax of the program train

            MAIN ::= (MAIN_HELP | MAIN_DO)
       MAIN_HELP ::= ./train [-h]
         MAIN_DO ::= ./train MODE DATASET
                             MODEL_PARAMS SAVE_MODEL
            MODE ::= LOAD_MODEL XVALIDATION
      LOAD_MODEL ::= -l MODEL.model
     XVALIDATION ::= -f N_FOLDS
         N_FOLDS ::= 1 | INTEGER
         DATASET ::= -d INTEGER
    MODEL_PARAMS ::= -k KERN -p KERN_PARAM
                     -C DOUBLE -g GAP_TOL -m TERM
            KERN ::= KERN_LIST | KERN_TYPE
       KERN_LIST ::= 0
       KERN_TYPE ::= 1 | 2 | 3 | 4 | 5 | 6 | 7
      KERN_PARAM ::= DOUBLE
         GAP_TOL ::= DOUBLE
            TERM ::= BITMASK
         BITMASK ::= 1 | 2 | 3
      SAVE_MODEL ::= -o MODEL.model

Both cross-validation and saving of a model can be performed at the same time
if wanted. However, this will mean that train will create one model for each
fold but it is just the last one that will be saved. If cross-validation is
not wanted pass one (-f 1) fold. The double precision floating point number
passed with -C is a number used in the classifier, it is related to the KKT
conditions. The gap tolerance is also a double precision floating point number
which is used as a convergence criterion. It is the allowed gap between the
primal and dual objective function, the feasibility gap, which should be a
small number. The default gap is set to $10^{-3}$. The m terminator is a bit-
mask which control when a classifier is considered optimal, i.e. when training
will stop. The feasibility gap constraint is not used if -m 2 is passed, i.e.
when the first bit (1) is zero. The primary training terminator bit is 2 which
means that all KKT conditions must be satisfied to terminate training. The
default of 3 means that both these conditions must be satisfied.

##  Appendix 6.B cellfeatures – Generate examples from the cell database

To generate features from all pairs of (image,ground truth class) in the cell
database the program cellfeatures is used. The file cellfeatures.data is
backed up before writing the features generated to it. This file can be used
by the program train.

           MAIN ::= (MAIN_HELP | MAIN_DO)
      MAIN_HELP ::= ./cellfeatures
        MAIN_DO ::= ./cellfeatures DB

##  Appendix 6.C jpeg_genfeature – Feature generation from images

To generate a set of features from image(s) the program called jpeg_genfeature
is used. It generate a set of features that can be classified later with
predict.

         MAIN ::= (MAIN_HELP | MAIN_DO)
    MAIN_HELP ::= ./jpeg_genfeature -?
      MAIN_DO ::= ./jpeg_genfeature CROPIMAGE* -o FEATURESET.feat
    CROPIMAGE ::= -i IMAGE.jpeg [-x left -y top -w width -h height]

##  Appendix 6.D predict – Predicting a set of features

To predict a set of features, generated by jpeg_genfeature, the program called
predict is used. It needs a previously trained model generated by train.

         MAIN ::= (MAIN_HELP | MAIN_DO)
      MAIN_DO ::= ./predict -l MODEL.model -f FEATURESET.feat
    MAIN_HELP ::= ./predict -?

##  Appendix 6.E extractcelltype – Extract a class of images from the cell
database

To extract a specific class (as classified by CellaVision™ DM-96) from the
cell database, the program called extractcelltype is used.

           MAIN ::= (MAIN_HELP | MAIN_DO)
      MAIN_HELP ::= ./extractcelltype
        MAIN_DO ::= ./extractcelltype CLASS DB
          CLASS ::= INTEGER
            DB  ::= ALLXMLFILES | (XMLFILE ’ ’)*
    ALLXMLFILES ::= ’.’

##  Appendix 6.F extractcellid – Extract given instances from the cell
database

To extract given instances from a list of id numbers, the program called
extractcellid is used.

           MAIN ::= (MAIN_HELP | MAIN_DO)
      MAIN_HELP ::= ./extractcellid
        MAIN_DO ::= ./extractcellid IDLIST DB
         IDLIST ::= (INTEGER ’ ’)* ’x’

##  Appendix 6.G extractcellinfo – Extract statistics of instances from the
cell database

To extract statistics about size, resolution and number of instances of a
specific class or of all classes the program called extractcellinfo is used.

           MAIN ::= (MAIN_HELP | MAIN_DO)
      MAIN_HELP ::= ./extractcellinfo
        MAIN_DO ::= ./extractcellinfo CLASS DB
          CLASS ::= CLASS_ALL | CLASS
      CLASS_ALL ::= ’-1’

##  Appendix 6.H tolibsvm – Save cell features in libSVM format

This program load the features saved in cellfeatures.data and dump them in
libSVM format on standard output. It takes no parameters.

    ./tolibsvm > cellfeatures.data.libsvm


# Spatial-temporal switching estimators for imaging locally concentrated
dynamics

###### Abstract

The evolution of images with physics-based dynamics is often spatially
localized and nonlinear. A switching linear dynamic system (SLDS) is a natural
model under which to pose such problems when the system’s evolution randomly
switches over the observation interval. Because of the high parameter space
dimensionality, efficient and accurate recovery of the underlying state is
challenging. The work presented in this paper focuses on the common cases
where the dynamic evolution may be adequately modeled as a collection of
decoupled, locally concentrated dynamic operators. Patch-based hybrid
estimators are proposed for real-time reconstruction of images from noisy
measurements given perfect or partial information about the underlying system
dynamics. Numerical results demonstrate the effectiveness of the proposed
approach for denoising in a realistic data-driven simulation of remotely
sensed cloud dynamics.

Index Terms—  localized dynamics, Kalman filter, reconstruction, denoising.

##  1 introduction

The image formation of time-varying phenomena is commonplace in astronomy [1],
remote sensing [2], and many other disciplines, e.g., [3]. The challenge is to
recover spatial-temporal parameters of a physical phenomenon given noisy, ex
situ measurements. Such time-dependent inverse problems can be posed in the
framework of a linear dynamic system (LDS), where the state variables
represent the volumetric physical parameters of interest which may be
reconstructed through Bayesian inference. Oftentimes, the system dynamics
traverse or cycle through a sequence of distinct modes, e.g., a sudden,
transient shock and then recovery towards a quiescent background state or some
pattern caused by an external driver. In such cases, it is natural to augment
the LDS with a hidden random variable to represent the system dynamic mode and
formulate the problem under a switching LDS (SLDS). A switching Kalman filter
(SKF) [4] is then the optimal approach to jointly estimate the system state
and detect the sequence of dynamic modes, but the enormous dimensionality
typically involved precludes such an approach.

Physical phenomena often exhibit strong spatial correlations, providing a
priori knowledge which may be utilized to dramatically reduce the SLDS
parameter dimensionality. For example, cloud dynamics are (approximately)
localized in space where the underlying physics are governed by local
variables such as humidity and temperature [5]. In this work, we harness this
prior information to segment the state space as a collection of disjoint
convex regions and assume that time evolution at a given point is correlated
only with immediate neighbors. Under these conditions, we employ a “divide and
conquer” strategy and develop a patch-based technique to reduce the switching
mode cardinality as well as the parameter space dimensionality.

Patch-based algorithms are often used in image and video processing for such
tasks as denoising and inpainting [6, 7, 8]. For example in video denoising,
[9, 10, 11] spatial-temporal patches are used for real-time denoising. To
further reduce computation, recursive patch-based denoising algorithms [12,
13, 14, 15, 16] use online clustering to determine the most similar patches
per frame, estimate patch cluster dynamics, and track patches by using optical
flow.

Generally speaking, conventional patch-based video processing algorithms
assume little to no prior information regarding the measurement or the
evolution model. Instead, they make use of the inherent redundancy found in
natural video frames when collected at high temporal frequency to determine
those (possibly disjoint) regions exhibiting a common dynamic mode and track
feature motion. In contrast, the physics-based imaging applications considered
here must rely on prior information to compensate for measurement
deficiencies, e.g., severe noise and sparsity, and feature tracking is
challenging or impossible because the hidden object could be diffuse and state
dynamics are complicated, e.g., advective flow as opposed to rigid body
motion.

The remainder of the paper is organized as follows. Section 2 introduces the
SLDS model which captures locally concentrated dynamics. The patch-based
estimator is proposed in section 3 and its computational complexity is studied
in Section 4. Section 5 compares the performance of the proposed filter to a
SKF in terms of efficiency and accuracy.

##  2 Locally concentrated hybrid state-space model

A spatial-temporal multi-modal LDS may be represented as

$\displaystyle\begin{bmatrix}\bm{x}_{n}{(1)}\\\ \vdots\\\
\bm{x}_{n}{(d)}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\bm{A}_{s_{1,1}}(1,1)&\dots&\bm{A}_{s_{1,d}}(1,d)\\\
\bm{A}_{s_{2,1}}{(2,1)}&\dots&\bm{A}_{s_{2,d}}(2,d)\\\ \vdots&\ddots&\vdots\\\
\bm{A}_{s_{d,1}}(d,1)&\dots&\bm{A}_{s_{d,d}}(d,d)\\\
\end{bmatrix}\begin{bmatrix}\bm{x}_{n-1}{(1)}\\\ \vdots\\\
\bm{x}_{n-1}{(d)}\end{bmatrix}$ $\displaystyle\quad\quad+\bm{\nu}_{n},$
$\displaystyle\bm{y}_{n}$
$\displaystyle=\bm{H}_{n}\bm{x}_{n}+\bm{\omega}_{n}.$

In the above, the state vector is $\bm{x}_{n}{(i)}$, where subscript $n$ and
$i$ represent the time and the state index, respectively, $d$ is the dimension
of the state, $\bm{A}$ is the evolution operator, and $s_{i,j}$ refers to a
hidden switching random variable that determines the correlation between
dynamic evolution of pixel $i$ with pixel $j$. The given information are
$\bm{y}_{n}$, the measurement, $\bm{H}_{n}$, the measurement operator, and
$\bm{\nu}_{n}\sim N(0,\bm{Q}_{n})$, $\bm{\omega}_{n}\sim N(0,\bm{R}_{n})$
[$N(\bm{\mu},\bm{\Sigma})$ refers to the Gaussian distribution with mean and
covariance $\bm{\mu},\bm{\Sigma}$] are the evolution and measurement noise.

When the observed state variable is a function of variables that are locally
concentrated in space, the evolution model in the neighboring pixels depend
only on a single hidden random variable. It is assumed that each pixel in the
image is correlated with its neighboring pixels with distance of at most $r$
pixels (which means $\bm{A}_{s_{i,j}}{(i,j)}=0\text{ if }|i-j|>r$). The state
dynamic equation is then

$\displaystyle\begin{bmatrix}\bm{x}^{(1)}_{n}\\\ \vdots\\\
\bm{x}^{(B)}_{n}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\bm{A}_{s_{1}}^{(1)}&...&0\\\
\bm{A}_{s_{2}}^{(2,1)}&...&0\\\ \vdots&\ddots&\vdots\\\
0&\dots&\bm{A}_{s_{B}}^{(B)}\\\
\end{bmatrix}\begin{bmatrix}\bm{x}^{(1)}_{n-1}\\\ \vdots\\\
\bm{x}^{(B)}_{n-1}\end{bmatrix}+\bm{\nu}_{n},$

where $B$ is the total number of patches, $\bm{x}_{n}^{(j)}$ refers to the
pixels at time $n$ in the $j^{th}$ patch that follow the same dynamic
behavior, $\bm{A}_{s}^{(i)}$ is the evolution model for patch $i$, and
$\bm{A}_{s}^{(i,j)}$ is the correlating term between the dynamic model of
patch $i$ with its neighboring patch $j\in\mathcal{B}_{i}$ where
$\mathcal{B}_{i}$ is the set of patches neighboring of patch $i$.

Fig. 1: Patch-based filtering vs. full filtering for the simulation setup; the
dynamics of each quarter switch based on a hidden random variable between two
modes (blue and yellow). The variance of the two modes is known. The full SKF
with perfect knowledge has access to the spatial extent of the modes (quarter
regions in this example), while the full SKF and (s)wSKF are equally ignorant
of mode spatial content.

##  3 Patch-based estimation

For an LDS, the optimal Bayesian estimator is the Kalman filter (KF). In this
section, we first review the KF/SKF procedure, as they are the building blocks
of the proposed algorithm. Patch estimators using multiple low-dimensional
SKFs are then proposed for efficient estimation of a spatial-temporal, locally
concentrated evolving dynamic system.

###  3.1 Kalman filter

With $\bm{y}_{1}^{n}$ defined as the set of all samples
$\bm{y}_{1},\bm{y}_{2},...,\bm{y}_{n}$, the initial state is $\bm{x}_{0}\sim
N(\bm{x}_{0|0},\bm{P}_{0|0})$,
$\bm{x}_{n|n}=\operatorname{\mathbb{E}}[\bm{x}_{n}|\bm{y}_{1}^{n}]$ and
$\bm{P}_{n|n}=\operatorname{\mathbb{C}}(\bm{x}_{n}|\bm{y}_{1}^{n})$ are the
conditional mean and covariance, and $L_{n}=p(\bm{y}_{n}|\bm{y}_{1}^{n-1})$
the likelihood, one step of the KF is

$\displaystyle(\bm{x}_{n|n},\bm{P}_{n|n},L_{n})$
$\displaystyle=Filter(\bm{A}_{n},\bm{H}_{n},\bm{x}_{n-1|n-1},\bm{P}_{n-1|n-1},$
$\displaystyle\quad\quad\bm{Q}_{n},\bm{R}_{n},\bm{y}_{1}^{n}),$ (1)

which involves the following two recursive steps:
Step 1: Time update equations

$\displaystyle\bm{x}_{n|n-1}$ $\displaystyle=\bm{A}_{n}\bm{x}_{n-1|n-1},$
$\displaystyle\bm{P}_{n|n-1}$
$\displaystyle=\bm{A}_{n}\bm{P}_{n-1|n-1}\bm{A}_{n}^{T}+\bm{Q}_{n};$

Step 2: Measurement Update equations

$\displaystyle\bm{e}_{n}$ $\displaystyle=\bm{y}_{n}-\bm{H}_{n}\bm{x}_{n|n-1},$
$\displaystyle\bm{B}_{n}$
$\displaystyle:=\bm{H}_{n}\bm{P}_{n|n-1}\bm{H}_{n}^{T}+\bm{R}_{n},$
$\displaystyle\bm{K}_{n}$
$\displaystyle=\bm{P}_{n|n-1}\bm{H}_{n}^{T}\bm{B}_{n}^{-1},$ $\displaystyle
L_{n}$ $\displaystyle=N(\bm{e}_{n};\bm{0},\bm{B}_{n}),$
$\displaystyle\bm{x}_{n|n}$
$\displaystyle=\bm{x}_{n|n-1}+\bm{K}_{n}\bm{e}_{n},$
$\displaystyle\bm{P}_{n|n}$
$\displaystyle=(\bm{I}-\bm{K}_{n}\bm{H}_{n})\bm{P}_{n|n-1}.$

For an SLDS, $\bm{A}_{n}$ and $\bm{Q}_{n}$ can change with time and they must
be detected as a component of the SKF using the obtained likelihoods ${L}_{n}$
using a Bayesian approach [4].

###  3.2 Patch estimator

The state-space equations for patch $i$ may be written as

$\displaystyle\bm{x}_{n}^{(i)}$
$\displaystyle=\bm{A}^{(i)}\bm{x}_{n-1}^{(i)}+\sum_{j\in\mathcal{B}_{i}}\bm{A}^{(i,j)}\bm{x}_{n-1}^{(j)}+\bm{\nu}_{n}^{(i)},$
$\displaystyle\bm{\Gamma}_{i}\bm{y}_{n}$
$\displaystyle=\bm{\Gamma}_{i}\bm{H}\bm{x}_{n}+\bm{\Gamma}_{i}\bm{\omega}_{n},$

where $\bm{A}^{(i)}$, $\bm{A}^{(i,j)}$, and $\bm{Q}^{(i)}$ are functions of a
hidden switching random variable corresponding to patch $i$ and neighboring
patch $j$, and $\bm{\Gamma}_{i}$ is an operator applied to the measurements in
order to find the localized measurements for process $i$. For localization, it
is sufficient to have
$\bm{\Gamma}_{i}\bm{H}\in\mathcal{T}=\begin{bmatrix}0,...,0\\\
0,\mathbf{\Theta^{(i)}},0\\\ 0,...,0\end{bmatrix}$, if such
$\bm{\Gamma_{i}},\bm{\Theta}_{i}$ exist. Otherwise, one must solve the
constrained minimization problem
$\operatornamewithlimits{argmin}_{\bm{\Gamma}_{i}}{\|\bm{\Gamma}_{i}\bm{y}_{n}-\bm{\Gamma}_{i}\bm{H}\bm{x}_{n}\|_{2}^{2}}$
s.t. $\bm{\Gamma}_{i}H\in\mathcal{T}$. A special case is when local
measurements are calculated for each pixel. In this case, the minimization
problem to obtain local measurements is equivalent to solving an inverse
problem, and the patch estimator may then be applied to the inverse problem’s
solution for denoising.

The estimation of a locally concentrated SLDS when the modes’ spatial extent
is unknown can be formulated as the following optimization problem for the set
of possible dynamic models, patch sizes, shapes, and state variables:

$\displaystyle\operatornamewithlimits{argmin}_{\bm{x}_{n},s_{i},\bm{U}_{i},\bm{U}_{j|i},B}$
$\displaystyle L(\bm{x}_{n},s_{i},\bm{U}_{i},\bm{U}_{j|i},B),$ (2)
$\displaystyle L$
$\displaystyle=\sum_{i=1}^{B}\sum_{j\in\mathcal{B}_{i}}\left\|\bm{U}_{i}\bm{x}_{n}^{(i)}-\bm{F}(s_{i},\bm{U}_{i}{\bm{x}}_{n-1})\right.$
(3)
$\displaystyle\quad\quad\quad\left.-\bm{F}(s_{i},\bm{U}_{j|i}\hat{\bm{x}}_{n-1})\right\|_{(\bm{Q}_{n}^{(i)})^{-1}}^{2}+\left\|\bm{H}\bm{x}_{n}-\bm{y}_{n}\right\|_{\bm{R}_{n}^{-1}}^{2},$

such that $\bm{F}(s,\bm{b})={\bm{A}(s,\bm{b})}\bm{b}$ applies the linear
operator ${\bm{A}(s,\bm{b})}$ to pixels in patch $\bm{b}$, $\bm{U}_{i}$ refers
to the operator selecting the $i^{th}$ patch such that patches are disjoint
and their intersection forms the whole image, and $\bm{U}_{j|i}$ selects the
patches correlated with the $i^{th}$ patch. The loss function (3) may be
(approximately) decoupled when
$\left\|\bm{H}\bm{x}_{n}-\bm{y}_{n}\right\|_{\bm{R}_{n}^{-1}}^{2}$ is replaced
by
$\sum_{i=1}^{B}\left\|\bm{\Theta}_{i}(\bm{U}_{i})\bm{x}_{n}^{(i)}-\bm{\Gamma}_{i}(\bm{U}_{i})\bm{y}_{n}\right\|_{\Delta_{i}[(\bm{\Gamma}_{i}\bm{R}_{n}\bm{\Gamma}_{i}^{T})^{-1}]}^{2}$,
where the matrix $\bm{\Gamma}_{i}$ is full rank and operator $\Delta_{i}$
selects the rows and columns of the matrix that correspond to the pixels in
patch $i$. For now, we assume the matrix $\bm{H}$ to be invertible, so
$\bm{\Gamma}_{i}$ always exists and is full rank. Solving the optimization (2)
for a general measurement operator is left for future work.

Once the patches’ structures are known, the optimization problem decouples
into a minimization for each patch with
$\operatornamewithlimits{argmin}_{s_{i},\bm{x}_{n}^{(i)}}L_{i}$ such that

$\displaystyle L_{i}$
$\displaystyle=\sum_{j\in\mathcal{B}_{i}}\left\|\bm{U}_{i}\bm{x}_{n}^{(i)}-\bm{F}(s_{i},\bm{U}_{i}{\bm{x}}_{n-1})\right.$
(4)
$\displaystyle\left.\quad\quad-\bm{F}(s_{i},\bm{U}_{j|i}\hat{\bm{x}}_{n-1})\right\|_{(\bm{Q}_{n}^{(i)})^{-1}}^{2}$
$\displaystyle\quad\quad\quad\quad+\left\|\bm{\Theta}_{i}\bm{x}_{n}^{(i)}-\bm{\Gamma}_{i}\bm{y}_{n}\right\|_{\Delta_{i}[(\bm{\Gamma}_{i}\bm{R}_{n}\bm{\Gamma}_{i}^{T})^{-1}]}^{2},$

which is equivalent to running a SKF for the patch $i$. The additional terms
corresponding to the neighboring patches are considered as inputs to the local
state space equations.

Solving the optimization problem (2) is an NP/hard problem since the loss
function needs to be calculated for each possible patch size, shape, and
location. Restricting the class of possible patch shapes to specific shapes
and sizes reduces the computational burden. In practice, given the set of
dynamic evolution models and $r$ (the correlation length), the image may be
divided into windows of size $2r\times 2r$ and run the SKF for each window; we
refer to this approach as the windowed SKF (wSKF). To reduce the effect of
decoupling at window boundary points, it is also possible to use windows of
size $(\alpha+2r)\times(\alpha+2r)$, slide the window over the whole image,
and estimate each pixel using the centered window; we refer to this approach
as the sliding windowed SKF (swSKF). The swSKF helps to reduce the boundary
effects of windowing. It is notable that the window size must be greater than
$2r\times 2r$ and that $\alpha>0$ must be chosen such that the windows are not
too small (overfitting) or large (underfitting).

The decoupled estimation covariance will be larger compared the full SKF with
knowledge of the mode spatial extent, since some information propagates
through the image when the patches are dependent. Since running a SKF for
every possible mode spatial extent is intractable, patch-based processing will
have some information loss.

Result: $\hat{\bm{x}}_{1},\hat{\bm{x}}_{2},...,\hat{\bm{x}}_{n}$

Input: $\mathcal{F}=\\{\bm{F}_{1},\bm{F}_{2},...,\bm{F}_{l}\\}$ s.t.
$\bm{F}_{i}\in\mathcal{F},\bm{H},\bm{R},\hat{\bm{x}}_{0},\hat{\bm{P}}_{0},r$,

$\mathcal{Q}=\\{\bm{Q}_{1},...,\bm{Q}_{l}\\}$ s.t. $\bm{Q}_{i}\in\mathcal{Q}$

for  _$n=1:T$_ do

for  _$i\in\\{1,...,B\\}$_ do

$\hat{\bm{x}}_{n}^{(i)}=SKF(\mathcal{F},\mathcal{Q},\bm{\Gamma}_{i}\bm{y}_{n},\bm{R},\bm{x}_{0}^{(i)},\bm{P}_{0}^{(i)})$

end for

$\hat{\bm{x}}_{n}=[\hat{\bm{x}}_{n};\hat{\bm{x}}_{n}^{(i)}]$

end for

Algorithm 1 Patch-based estimation

(a)

(b)

Fig. 2: LABEL:sub@fig:cloud Realizations of the ground truth dynamical cloud
images and the measurements, as well as the estimates using different filters.
Visually, it is obvious that the (s)wSKF performs better than the full SKF
with the same amount of information and with far fewer computations.
LABEL:sub@fig:MSE The estimation MSE for the filters. The (s)wSKF has far less
MSE compared to the full SKF with the same amount of information.

##  4 Computational cost

Consider the estimation of a $\sqrt{d}\times\sqrt{d}$ image with $K$ non-
overlapping processes that switch between $l$ modes over time. Thus, an SKF
that has perfect knowledge of mode spatial extent must consider $K^{l}$
different modes at each time step with computational complexity
$\mathcal{O}(K^{l}d^{3})$. On the other hand, decoupling the image into $K$
windows where each window switches between $l$ modes only requires
$\mathcal{O}(Klr_{w}^{3})$ where $\sqrt{r_{w}}\times\sqrt{r_{w}}$ is the
window size. (For an illustration, see Fig.1). Similarly, running an swSKF
requires $\mathcal{O}(K^{\prime}lr_{w}^{3})$, where $K^{\prime}$ is the number
of sliding windows and $K^{\prime}>K$. Thus, the computational requirement of
the patch-based SKF is much smaller than that of the full SKF, even given
perfect knowledge of the spatial extent of the modes, and this difference
becomes more significant as $d$, $K$, and $l$ increase.

##  5 Simulation results

The performance of the proposed patch-based estimator for denoising
reconstructed video frames from tomographic measurements is studied in this
section. A sequence of $32\times 32$ images are reconstructed in time using a
nonlinear, locally switching dynamic model, and a set of measurement with an
average SNR of 11dB is generated accordingly; these images are meant to
represent the evolution of cloud density obtained by infrared imaging in which
clear local concentration is known to exist. The data are generated such that
the structures’ elements have random movements with slow/fast velocities
(movement velocities are $0.01,0.94$ pixels per time step) in each quarter of
the image (as in Fig. 1), and the velocity can switch randomly in each quarter
over time.

The goal is to denoise the set of tomographically reconstructed images.
Because the system dynamics are complicated, we use a purely stochastic model
with two evolution covariance matrices corresponding to slow and fast
velocities such that the variances are roughly set to the velocity of
movements in the two cases. Assuming knowledge of the variances only, we apply
a bi-modal full SKF, wSKF, and swSKF with window size $8\times 8$ to the
sequence of images. Fig.2a shows realizations of the ground truth, noisy
reconstructed image, and the estimates for three representative time steps.
The full SKF clearly cannot recover the ground truth details as well as the
wSKF/swSKF. This conclusion is even more pronounced in Fig.2b, where the error
for each filter is averaged over 100 realizations and the full SKF MSE using
the bi-modal model tends to diverge since the model used is mismatched. It is
also notable that the swSKF error is slightly smaller compared to the wSKF.
The (s)wSKF can find the estimates without the knowledge of the spatial extent
of the modes and by detecting it instead while requiring fewer computations
using the locality assumption. The run time of the filters are presented in
the table below.

| full SKF | wSKF | swSKF
---|---|---|---
Run time(sec) | 68.52 | 1.26 | 4.18

##  6 Conclusion

A patch-based filtering framework is proposed for the estimation of image
sequences governed by locally concentrated dynamics and shown to have superior
performance with respect to the full filter in terms of computation and
accuracy, when provided perfect/partial information about the evolution model.
Solving the reconstruction problem for a general measurement operator, fusing
the local estimates to obtain global estimates, and online learning of the
evolution models/statistics of the patches locally are future work.

## References

  * [1] M. D. Butala, F. Kamalabadi, R. A. Frazin, and Y. Chen,  “Dynamic tomographic imaging of the solar corona,”  IEEE Journal of Selected Topics in Signal Processing, vol. 2, no. 5, pp. 755–766, Oct 2008.
  * [2] D. A. Galvan, A. Komjathy, M. P. Hickey, P. Stephens, J. Snively, T. Song, M. Butala, and A. J. Mannucci,  “Ionospheric signatures of Tohoku-Oki tsunami of March 11, 2011: Model comparisons near the epicenter,”  vol. 47, pp. RS4003, 2012.
  * [3] A. K. George, M. D. Butala, R. A. Frazin, F. Kamalabadi, and Y. Bresler,  “Time-resolved CT reconstruction using the ensemble Kalman filter,”  in Proc. IEEE ISBI, Paris, France, 2008, pp. 1489–1492.
  * [4] Kevin P. Murphy,  “Switching Kalman filters,”  Tech. Rep., University of California Berkeley, 1998.
  * [5] B. Thurairajah and J. A. Shaw,  “Cloud statistics measured with the infrared cloud imager (ici),”  IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 9, pp. 2000–2007, 2005.
  * [6] P. Chatterjee and P. Milanfar,  “Patch-based near-optimal image denoising,”  IEEE Transactions on Image Processing, vol. 21, no. 4, pp. 1635–1649, 2012.
  * [7] P. Chatterjee and P. Milanfar,  “Practical bounds on image denoising: From estimation to information,”  IEEE Transactions on Image Processing, vol. 20, no. 5, pp. 1221–1233, 2011.
  * [8] A. Criminisi, P. Perez, and K. Toyama,  “Region filling and object removal by exemplar-based image inpainting,”  IEEE Transactions on Image Processing, vol. 13, no. 9, pp. 1200–1212, 2004.
  * [9] B. Wen, S. Ravishankar, and Y. Bresler,  “Video denoising by online 3d sparsifying transform learning,”  in 2015 IEEE International Conference on Image Processing (ICIP), 2015, pp. 118–122.
  * [10] M. Maggioni, G. Boracchi, A. Foi, and K. Egiazarian,  “Video denoising, deblocking, and enhancement through separable 4-d nonlocal spatiotemporal transforms,”  IEEE Transactions on Image Processing, vol. 21, no. 9, pp. 3952–3966, 2012.
  * [11] Thibaud Ehret., Pablo Arias., and Jean-Michel Morel.,  “Global patch search boosts video denoising,”  in Proceedings of the 12th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications - Volume 4: VISAPP, (VISIGRAPP 2017). INSTICC, 2017, pp. 124–134, SciTePress.
  * [12] A. Buades, J. Lisani, and M. Miladinović,  “Patch-based video denoising with optical flow estimation,”  IEEE Transactions on Image Processing, vol. 25, no. 6, pp. 2573–2586, 2016.
  * [13] A Buades and J L Lisani,  Patch-Based Methods for Video Denoising, pp. 175–205,  Springer International Publishing, Cham, 2018.
  * [14] A. Davy, T. Ehret, J. Morel, P. Arias, and G. Facciolo,  “A non-local cnn for video denoising,”  in 2019 IEEE International Conference on Image Processing (ICIP), 2019, pp. 2409–2413.
  * [15] T. Ehret, J. Morel, and P. Arias,  “Non-local kalman: A recursive video denoising algorithm,”  in 2018 25th IEEE International Conference on Image Processing (ICIP), 2018, pp. 3204–3208.
  * [16] P. Arias and J. Morel,  “Kalman filtering of patches for frame-recursive video denoising,”  in 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), 2019, pp. 1917–1926.