\chapter{Feature extraction}	\label{sec:feature}
\thispagestyle{empty}

  
In order to recognize an object in an image, it must be extracted some kind of measurements from it. These values should express the characteristics and
all the information contained in the image. Moreover, we have to look for the smallest set of values with these good properties in order to
speed up the recognition and to achieve more reliability. This combination of input data is called feature and it is based on some
understanding of the particular problem being tackled or by some automatic procedures, \cite{libro:Bishop}. It is the heart of the
classification system. If the features extracted are not representative of input images, then, no classifier will be able to do a good job.

In this chapter, we describe two different kinds of feature. The first one, is based on the
research of Schmid et al. \cite{paper:Schmid}. The latter one is derived from a work of Torralba et al. \cite{paper:ipcspc}, in
which they aim to discriminate between scenes containing man-made objects and scenes with natural objects. 
In the next chapter, we will show the very good results of classifiers using these two sets of
features.



  
  
\section{Local Jets}

In \cite{paper:Schmid}, the problem of matching an image in a large data set was addressed. In this method, from each interest point found in
a image (using Harris' detector) a set of local gray-value invariants were computed. Our features are based on these invariants, originally
studied by Koenderink et al. \cite{paper:invarianti}.

We want to extract from each image pixel a set of values, which are invariants with respect to shift and rotation. For example, we aim to
acquire the same measurements from two images with the same rotated and shifted bacterium. Furthermore, we want to get the same values when
the same object is at different distance from the camera, property referred as \emph{scale invariance}. Even though we are not interested in
specific object recognition, the scale, shift and rotation invariance is important in our classification task. With this kind of features, we
are able to cluster together the measurements of images belonging to a certain class because we acquire the same measurements, no matter the
position and the orientation of the cell and these values are robust against little variations in size. 

Let be $I(x,y)$ an image with spatial coordinates $x$ and $y$. We define \emph{local jet} of order $N$ the convolution of image $I$ with the 
derivatives with respect to $(\overbrace{x,\cdots,y,\cdots)}^{N \mbox{ times}}$ of a Gaussian kernel with standard deviation $\sigma$. We
indicate these local jets with $L_{i_1 \cdots i_n}$, $i_k \in \{x,y\}$. The subscripts refer to the kind of Gaussian kernel
derivation. Note that the convolution of $I$ with a derivative of Gaussian kernel is equal to the convolution of Gaussian kernel with the
same derivative of $I$.\\
In order to obtain invariance with respect to shifts and rotations, linear (or more precisely a differential) combinations of local jets are
computed. Our set of invariants is limited to 9, that is to a combination of third order derivatives. These invariants are\footnote{The first
invariant is the average luminance, the second one is the square of the gradient magnitude, the fourth one is the Laplacian. }:

\begin{eqnarray} \label{eq:invarianti}
I_0 & = & I  \nonumber \\
I_1 & = & L_x^2 + L_y^2  \nonumber \\
I_2 & = & 2 L_x L_{xy} L_y + L_x^2 L_{xx} + L_y^2 L_{yy} \nonumber \\
I_3 & = & L_{xx} + L_{yy}  \nonumber \\
I_4 & = & L_{xx}^2 + L_{yy}^2 + 2 L_{xy}^2 \nonumber \\
I_5 & = & L_{xxx} L_y^3 - L_{yyy} L_x^3 + 4 L_{xyy} L_x^2 L_y - 4 L_{xxy} L_x L_y^2 \nonumber \\
I_6 & = & L_{xxy} L_y^3 - L_{xyy} L_x^3 + L_{xxy} L_x^2 L_y - L_{xyy} L_x L_y^2 \nonumber \\
I_7 & = & L_{xyy} L_y^3 - L_{xxy} L_x^3 + L_{xxx} L_y L_x^2 - L_{yyy} L_x L_y^2 + 2 L_{xxy} L_x L_y^2 - 2 L_{xyy} L_x^2 L_y \nonumber \\
I_8 & = & L_{xxx} L_x^3 + L_{yyy} L_y^3 + 3 L_{xxy} L_x^2 L_y + 3 L_{xyy} L_x L_y^2
\end{eqnarray}

In order to deal with scale changes, these 9 invariants are computed at different scales. At each stage a blurred image (with a Gaussian
kernel) is considered as starting point to derive the invariant images. In our experiments, we considered 4 scales that is, we compute
invariants using the original image, its blurred version, the blurred version of this last version and so on for another stage. Totally,
given an image we build a stack with $9 \times 4 = 36$ images. \figref{esempioLJSintetico} shows an example of this images of invariants
using synthetic data: impulse (11x11 pixels) in first and second row, square and rhombus (100x100 pixels) in third, fourth and fifth, sixth row
respectively; the gray level values are (approximately) proportional to the invariants values. \figref{esempioLJ} shows all the 36 
invariants for a bacterium image.\\

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth,height=0.7\textwidth]{./figure/esempioLJSintetico.eps}
\caption{Invariants at scale 0 and 1 for synthetic data: impulse (11x11 pixel image), square and rhombus (100x100 pixel image). Note that some
asymmetries are due to the approximation used to represent a value in gray-level and the rotation of the square in the rhombus is
not perfect due to the poor resolution used.}
\label{figure:esempioLJSintetico}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth,height=0.7\textwidth]{./figure/esempioLJ.eps}
\caption{Example of invariants at scale $[0,\cdots,3]$ for an image of bacterium.}
\label{figure:esempioLJ}
\end{figure}

Till now, we have followed the same approach of \cite{paper:Schmid} to extract from pixels feature vectors. We now introduce a new stage
 that is, we apply some non-linearities. \\
Given an image of invariants, we compute the background mean from the
border (a ten pixel wide frame) and we split the image in three derived images: the positive and negative part and the absolute value with
respect to the background mean. Thus, each image gives 36 images of invariants that are splitted in 108 images thanks to this method.\\
Finally, we make an average of all the pixel values in each of the 108 planes and we get from each image a single feature vector with 108
components. 

Just to have an idea of how much these features are able to separate our particles, consider \figref{esempioClassificazione2}. 
Here, we have taken only 100 images from each of two pollen categories and we have extracted features with only two components (without
non-linearity application).
However, the separation between the data is already quite good.\\

\begin{figure}[hbtp]

  \centerline{\hbox{ \hspace{0.0in} 
  \includegraphics[width=0.45\textwidth,height=0.385\textwidth]{./figure/exampleFeatLJpoll.eps}
    \hspace{0.3in}
  \includegraphics[width=0.45\textwidth,height=0.385\textwidth]{./figure/pollenfeatspace.eps}
    }
  }

  \vspace{9pt}
  \hbox{\hspace{1.35in} (a) \hspace{2.10in} (b)} 
  \caption{(a) Left column: two kinds of pollen grains, pecan and rumex.
Central and right columns: processing on previous images in order to get from each image a feature vector with two components.
(b) Example of classification: points in feature space for pecan and rumex pollen, cross and circle respectively. A linear classifier is able to reach an error rate of about 15\%.
 }
  \label{figure:esempioClassificazione2}
\end{figure}

We found experimentally that the applications of these non-linearities boosts the performance of the classifier. For example, when we trained the system on the first 500 images of each class of urine database and we did not apply the nonlinearities, we achieved an error rate of 18\%. When we applied also the nonlinearities, the error rate decreased to 12\%.


\subsection{Decreasing the dimensionality}

In training, we aim to model the distribution of feature vectors belonging to all the images of a certain class with a mixture of Gaussians
(MoG) or with support vector machines (SVM) (see par. 
\ref{sec:mog}, \ref{sec:svm}). However, it is not possible to work in a 108 dimensional feature space because the number of parameters grows very
quickly with the space dimensions. For example, in \tabref{nrparamMoG} you can find these relations as a function of the choice of covariance matrix
structure and as a function of the number of components in MoG.

\begin{table}[htbp] 
\begin{center}
\begin{tabular}{||p{2.5cm}||*{3}{c|}|}
\hline
	 & \scshape{spherical} & \scshape{diagonal} & \scshape{full}\\
\hline 
\hline
& & & \\
\emph{nr. param.} & $N(D+1)+N-1$ & $2ND + N -1$ & $N(D + \frac{D(D+1)}{2}) + N - 1$ \\
& & & \\

\hline
\hline
\end{tabular}
\end{center}
\caption{Number of parameters in MoG using $N$ components in a $D$ dimensional feature space as a function of the covariance matrix structure.}
\label{tab:nrparamMoG}
\end{table}

For instance, if $D = 108$ and $N = 2$ then, we should estimate $220,433,11989$ parameters in accordance to the chosen covariance matrix
structure. On the other hand, a common rule of thumb says that you need at least
3 or 4 points to estimate a single parameter. Since our database has not a so huge number of images we have to reduce the
dimension of feature space.

In order to reduce the feature space dimensionality we could apply PCA \cite{libro:pca} \cite{libro:Bishop}. Experimentally, we found that
 linear discriminants analysis works better than PCA. Particularly, given the whole set of training feature vectors we compute Fisher's Linear
 Discriminants (FLD) \cite{libro:Bishop} \cite{paper:fld}.
 
 We now introduce this technique. Suppose for simplicity to have data ${\xv} \in \Bbb{R}^D$ belonging only to two classes $C_1$ and $C_2$. We
 want to find the vector ${\wv}$ such that the projected data
 
 \begin{equation} \label{eq:proj}
 y_n = {\wv}^T {\xv}_n
 \end{equation}
 
 are optimally separated. The mean vectors of the two classes are
 
 \begin{equation}
 {\mv}_1 = \frac{1}{N_1} \sum_{n \in C_1} {\xv}_n \hspace{1cm} {\mv}_2 = \frac{1}{N_2} \sum_{n \in C_2} {\xv}_n
 \end{equation}
 
 where $N_i$ is the number of items in class $i$-th.\\
Then, we define

\begin{eqnarray}
m_k = {\wv}^T{\mv}_k \nonumber \\
s_k^2 = \sum_{n \in C_k} (y_n - m_k)^2
\end{eqnarray} 

Fisher's discriminants maximize a function which is in the two-class problem,

\begin{equation}
\frac{(m_2-m_1)^2}{s_1^2 + s_2^2}
\end{equation}

that is, we are trying to maximize the distance between the means of the projected clusters and minimize the spread of the projected data. While
PCA seeks the sub-space that best \emph{represents} the data, FLD seeks the sub-space that maximizes the \emph{separability} between the classes. \\
The generalization of FLD to several classes follows the same reasoning. It turns out that if you have $K$ classes you can have no more than
$K-1$ directions where to project the data. For example, if we think to work with the 12 categories of urine database, we should have the
constraint to reduce the feature space dimension from 108 to less than 12. In order to reduce the feature space dimensions without any
constraint, we split randomly the data of each class in a proper number of sub-categories \footnote{We tried also to divide each cluster applying
\emph{k-means} algorithm \cite{libro:Bishop} and estimating a MoG but without any meaningful improvement.}.



\section{Image and Power Spectrum Principal Components}	\label{sec:ipcspc}

We describe a feature extraction based on principal component analysis (PCA), \cite{libro:pca} \cite{libro:Bishop}. \\
Given data in a $N$ dimensional space, the goal of PCA is to find a $D$ dimensional subspace such that
the projected data are the closest in $L_2$ norm (mean square error) to the original data. This means that we are
looking for the subspace which gives the best \emph{representation} of the original data. This subspace is spanned
by the eigenvectors of the data covariance matrix having the highest corresponding eigenvalues. Often the full covariance 
matrix can not be reliably estimated from the number of example available, but the approximate highest eigenvalue basis
vectors can be computed using singular value decomposition (SVD). Thus, if $\Sigma$ is the covariance matrix of the data, 
its SVD decomposition is

\begin{equation} \label{eq:svd}
\Sigma \stackrel{svd}{=} USV^T
\end{equation}
  
where if $\Sigma \in \Bbb{R}^{N \times N}$, then $U \in \Bbb{R}^{N \times N}$, $V \in \Bbb{R}^{N \times N}$ and $U U^T = I$, 
$V V^T = I$. $S$ is a diagonal matrix with the elements on the diagonal (the singular values) in descending order. The first columns
of $U$ are the eigenvectors with the highest eigenvalues and can be used as basis for the original data.\\
In our specific case, let be $i(x,y)$ the intensity distribution of the image along spatial variables $x$ and $y$ with $x,y \in [1,N]$.
Let be ${\iv}$ the values of image pixels rearranged in a column vector, then the covariance matrix is $ \Sigma =
E[({\iv}-{\mv})({\iv}-{\mv})^T] $ with $ {\mv} = E[{\mv}] $. If we call image principal components $IPC_n$ the eigenvectors 
of $\Sigma$ (computed by singular value decomposition) and $IPC_n(x,y)$ the same eigenvectors reshaped into a matrix $N \times N$, then 
we can write the following decomposition,

\begin{equation} \label{eq:ipc}
i(x,y) = \sum_{n=1}^P v_n IPC_n(x,y)
\end{equation}

with $P \leq N^2$ and $v_n$ the coefficients used for describing the image $i(x,y)$ in this new basis. Obviously, $v_n$ is the projection of
the image along the direction given by $IPC_n$. \\
\figref{esempioIPC} shows 36 IPCs (with the highest corresponding eigenvalues) computed from the first 500 images of each class belonging to
the urine database. In our experiments we used always less than 10 IPCs and you can note from \figref{andamentoSingValueIPC} that the first
components are the most important in the data representation. Before any computation, we subtract the background mean estimated from the border of the image
(a frame 10 pixels wide) in order to normalize the data with respect to the background brightness.\\

\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{./figure/esempioIPC.eps}
\caption{First IPCs derived collecting 500 images from each class in urine database. Each image was normalized to a common resolution of 40x40
pixels, thus the space dimension is 1600 and there are 1600 IPCs and singular values.}
\label{figure:esempioIPC}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth,height=0.5\textwidth]{./figure/andamentoSingValueIPC.eps}
\caption{The singular values corresponding to IPCs indicate the importance of each component in representing the data.}
\label{figure:andamentoSingValueIPC}
\end{figure}


If we compute the discrete Fourier transform (DFT) $I(k_x,k_y)$ of an image, we can approximate its power spectrum with the squared
magnitude of $I$:

\begin{equation}
R(k_x,k_y) = \left| I(k_x,k_y) \right|^2 
\end{equation}

with

\begin{equation}
I(k_x,k_y) = \frac{1}{N^2} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} i(x,y) e^{ -j\frac{2\pi}{N}(x k_x + y k_y) }
\end{equation}

and $f_x = k_x/N$, $f_y = k_y/N$ spatial frequencies. The power spectrum encodes the energy density for each spatial frequency and
orientation over the whole image.
PCA applied to power spectra gives the main components that take into account the structural variability between images. \\
First, we normalize the power spectrum with respect to its variance for each spatial frequency: $R^*(k_x,k_y) = R(k_x,k_y) / std[ R(k_x,k_y) ]$,
with $ std[ R(k_x,k_y) ] = \sqrt{ E[ (R(k_x,k_y) - E[R(k_x,k_y)])^2 ] } $. Then, we apply PCA as done before for the image, we find the
spectral principal components (SPCs) and we decompose the normalized power spectrum in

\begin{equation} \label{eq:spc}
R^*(k_x,k_y) = \sum_{n=1}^P u_n SPC_n(k_x,k_y)
\end{equation}


\figref{esempioSPC} shows the SPCs (with the highest corresponding eigenvalues) computed from the first 500 images of each class belonging to
the urine database.\\

\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{./figure/esempioSPC.eps}
\caption{First SPCs derived collecting 500 images from each class in urine database. Each image was normalized to a common resolution of 40x40
pixels.}
\label{figure:esempioSPC}
\end{figure}

Till now, we have described the basic approach as presented in \cite{paper:ipcspc}. In our research, we added a new step.\\
After the normalization of the image with respect to the background mean that makes this mean equal to zero,
 we split the image in three derived images: the first one is its positive part $I^+$, the second one its negative part $I^-$and the third one 
 the absolute value $I^a$. This operation is equivalent to the application 
 of three non-linear functions that are selective of high, low and mid-range responses.
 Then, we apply for each of these derived images the PCA and we get $IPC_+,IPC_-,IPC_a$. A set of these principal components are shown in
 \figref{esempioIPCNL} and are computed from the same data used to calculate the basis shown in \figref{esempioIPC}.\\
 It is experimentally proved that the classification is improved when this splitting is done; evidently, particles belonging to different
 categories behave differently with respect to this operation. Note that the split in positive, negative and absolute value parts can be done
 very efficiently. We are adding complexity only in training stage because we double the number of PCA required (from two to four).\\
With this trick, we decrease the error rate from nearly 19\% to 16\% with reference to the classification of urine database.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{./figure/esempioIPCNL.eps}
\caption{First 12 $IPC_+,IPC_-,IPC_a$ derived collecting 500 images from each class in urine database. Each image was normalized to a common resolution of 40x40
pixels.}
\label{figure:esempioIPCNL}
\end{figure}

If we consider the first $N_1^i$,$N_2^i$,$N_3^i$,$N^s$ coefficients computed by the projection of $I^+$ along $IPC_+$, $I^-$ along $IPC_-$, 
$I^a$ along $IPC_a$ and $R^*$ along $SPC$ we get a feature vector from each image

\begin{equation} \label{eq:featureIPCSPC}
{\fv} = [v_{11}, \cdots v_{1 N_1^i}, v_{21}, \cdots v_{2 N_2^i},v_{31}, \cdots v_{3 N_3^i},u_1, \cdots u_{N^s}]
 \end{equation}
 
where $v_{ij}$ is the projection of the image $I^i$ along j-th $IPC_i$ for $i\in\{+,-,a\}$.\\

How do we choose values for $N_1^i,N_2^i,N_3^i,N^s$? In preliminary experiments, we found that the nearly same performance can be achieved
considering 

\begin{eqnarray}
N_1^i = N_2^i = N_3^i \nonumber \\
N^s = 3 N_j^i, \hspace{.5cm} \mbox{ for } j=1,2,3
\end{eqnarray}

and applying Fisher's linear discriminants to the whole set of principal components in order to reduce the dimension of
feature space to $N=N_1^i+N_2^i+N_3^i+N^s$. In the next chapter, we will only consider the first option because of its simplicity. 
  


\section{Interpretation} 

The two sets of features described in this chapter are similar, in spite of the different techniques used to compute them. 
Features based on PCA give the average information content of the signal at frequencies higher and higher when it is considered the 
projections on eigenvectors with smaller and smaller eigenvalues. 
Similarly, features based on invariants give the global amount of information 
at different bands when it is taken the average of these invariants at different scales. This is proved experimentally. The implementation of 
classifiers based on these features skipping the non-linearity application gave nearly the same performance, with an error rate of 20\% on
the twelve categories of urine database.\\
The application of non-linearities boosts the performance; in particular, it nearly halves the error rate of classifier using the features 
based on local jets.



  
  
\clearemptydoublepage



