\section{Introduction}

Developing chemo-sensing solutions and standards for early warning against chemical and biological hazards   has been an active research area~\cite{chemHazards,guideHazards}. To construct an accurate and reliable chemical warning system, the information from several sensors must be integrated to provide a clear indication for the composition of the chemical substances as well as their propagation profile. We propose a machine learning approach, that is based on recent advances, for classification and regression problems in machine olfaction. Machine olfaction problems include odor classification problems, gas consecration detection, and chemical source localization.

Our algorithm consists two steps of feature generation and classification. To transform the raw data of the sensor array platform into discriminative features, we propose the Scattering Time Series on Graphs(STSG) transform, which is an hierarchical feature extraction method from multiple time series. This transform is an extension of the recently proposed scattering transform for time series~\cite{mallat2012group,bruna2013invariant,bruna4104multiscale,bruna2013audio} and graph signals~\cite{chen2014unsupervisedHaarScatGraphs}, to multivariate time series defined on a graph. The resulting features are then classified using a random forest (RF) based classifier~\cite{breiman2001random, breiman1984classification}.

We demonstrate and evaluate our algorithm by means of the “Dataset from chemical gas sensor array in turbulent wind tunnel”~\cite{fonollosa2015dataset}, available through the UCI Machine Learning Repository~\cite{webUCIgasSensors}.
This dataset has been used to validate existing algorithms~\cite{vergara2013performance, vembu2012time}.

The structure of our paper is as follows. In Section~\ref{ch:dataset} we describe a set of machine olfaction problems we wish to solve and the dataset that are used. In Section~\ref{sec:sateOFart} we review prior art methods.

In Section~\ref{ch:Math}, we present the theoretical building blocks of our method: a redundant Haar wavelet decomposition over a (possibly irregular) graph, the scattering convolution network and Random Forests. Equipped with these building blocks, we present in Section~\ref{ch:MainAlgo} the main contribution of this paper, the \emph{\textbf{STSG}} (Scattering of a Time Series on a Graphs) algorithm.

Finally, in Section \ref{ch:Exper}, we show experimental results of our methodology applied in the machine olfaction setting: odor discrimination, odor consecration, and odor localization. We compare the performance of our method with prior techniques of machine olfaction~\cite{vergara2013performance,vembu2012time}.

\section{Olfaction datasets} \label{ch:dataset}
In this paper we use the dataset ``Gas sensor arrays in open sampling settings" from the UCI archive (Machine learning Repository)~\cite{webUCIgasSensors}. This dataset is a collection of multidimensional time series data along with some static environmental parameters that include the responses of a chemical detection platform to different gases at different levels of concentration. The challenge is to develop machine learning techniques for gas classification and source prediction models. In this section, we review the dataset in details.

\begin{figure}
  \centering
  \includegraphics[scale=0.5]{Wind_tunnel.eps}\\
  \caption{Wind tunnel used to collect time series data from sensor arrays\cite{fonollosa2015dataset}.}\label{fg:Wind_tunnel}
\end{figure}

The data in~\cite{webUCIgasSensors} was collected in a  $2.5m \times 1.2m \times 0.4m$ wind tunnel test-bed facility (see Figure~\ref{fg:Wind_tunnel}), into which the gaseous substances of interest were released.

The wind tunnel operates in a propulsion open-cycle mode, by continuously drawing external turbulent air throughout the tunnel and exhausting it back to the outside, creating a relatively less turbulent airflow moving downstream towards the end of the test field. In order to construct various distinct artificial airflows in the wind tunnel, the wind tunnel contains a motor-driven exhaust fan at the outlet of the test section. The motor can be set to rotate at three different constant speeds. The wind tunnel measures the ambient temperature and relative humidity during the entire experiment.

The chemical detection platform inside the wind tunnel consists of boards. Each board has eight commercialized metal-oxide gas sensors (MOX)~\cite{FigaroIncSite}, which are sensitive to rapid changes in the analytes concentration. Thus, the output of each board is an $8$-dimensional time series. The chemical detection platform consists of columns of nine boards each located at six equally-spaced positions along the wind tunnel, that is, a total of $72$ sensors per location (see Figure~\ref{fg:Wind_tunnel}). Figure~\ref{fg:response} depicts a typical time-series response of one board. The sensor responses are affected by the air turbulence in the wind tunnel and depend on the concentration of the gas substance. As the operating temperature of the sensors affects their performance, it is adjustable by setting the voltage of the built-in heater of each sensor to one of five different levels.

\begin{figure}
  \centering
  \includegraphics[scale=0.75]{Multivariate_response.eps}\\
  \caption{Multivariate response of a 8-sensor array when methane is released in the wind tunnel~\cite{fonollosa2015dataset}.}\label{fg:response}
\end{figure}

The dataset~\cite{webUCIgasSensors} was generated by releasing ten different types of gas into the wide tunnel: acetone, acetaldehyde, ammonia, butanol (butyl-alcohol), ethylene, methane, methanol, carbon monoxide, benzene, toluene, and carbon monoxide(CO).

Each kind of chemical substance is released at same  nominal  concentration  values  at  the  outlet  of  the gas  source  in  parts-per-million  by  volume(ppmv).The CO  was released in two different nominal concentrations. The value of  concentration at gas  source is from $100$ ppm to $400$ ppm. Note that the actual concentration in the wind tunnel decreases as the generated gas plume spreads out along the tunnel.

For each gas released, the motor of the exhaust fan was set to one of three rotation speeds. A  turbulent airflow was thereby generated within the wind tunnel.
The outputs of the sensors at each of the six locations in the wind tunnel were measured separately, resulting in six $72$-dimensional time series datasets capturing the chemical analyte circulated throughout the wind tunnel. The temperature of all sensors remained fixed during the measurement. After each test, the wind tunnel was thoroughly ventilated. The measurement was repeated for three different air velocities and five operating temperatures.
 For each combination of the type of gas, airflow velocity, and temperature, the measurement (generation of time-series) was repeated $20$ times. The time-series measure at each location was sampled for approximately $250$ seconds, with $10$ samples per second, that is, a total of about $2500$ samples per location per experiment. The time series measured at all locations was fully synchronized, although each was of a slightly different length. The data for each of the ten gas classes was collected $20$ times for each three speeds, five temperatures modes and six locations. The resulting dataset thus consists of 18,000 72-dimensional time series.

   In addition, each experiment recorded the ambient temperature and relative humidity during the experiment. Although these parameters are also represented as time series, they have very low variance and can thus be represented by their average values. For a more detailed description of the experimental protocol see~\cite{fonollosa2015dataset,vergara2013performance}.

\section{Related studies on odor classification}\label{sec:sateOFart}
Several previous studies~\cite{nimsuk2007improvement, nimsuk2008study, vembu2012time, vergara2013performance} have considered the problem of odor classification. These studies all follow the same approach of extracting features from the data followed by some classification scheme.

Earlier approaches of odor classification problems~\cite{nimsuk2007improvement, nimsuk2008study} use features extracted by applying stepwise discriminant analysis (DA) to the short-time Fourier transform (STFT), followed by classification of these features with an LVQ (learning vector quantization) neural network. The Fourier power spectrum of time series $X(t)$ is defined by
$$
	\mathcal{F}_{w,\phi} X(t,\xi):=\left|\cap X(\tau)w(t-\tau) \Exp{-i\xi\tau} d\tau \right|^2\ast\phi(t),
$$
where $w$ is a short-time low pass filter, $\phi(t)$ is a normalized smoothing window, and  $\xi$ is a frequency value. The Fourier power spectrum is defined as the expected value of square of the modulus of Short Time Fourier Transform(STFT $L_2$-moments)

\begin{equation}\label{fr:Fourier}
\bar{\mathcal{F}_{w}}(\xi) = {\mathbb{E}}_{t} \left[\mathcal{F}_{w,\phi} X(t,\xi)\right].
\end{equation}

Then, the STFT-based feature mapping of time series  $X(t)$ is
$$
\Phi: X \to \bar{\mathcal{F}_{\omega}}(\xi), \qquad \xi \in \Xi,
$$
where $\Xi$ is a given set of frequencies.

As the STFT feature extraction method constitutes prior art for our scattering approach, we applied it to our data as a performance benchmark.

Recent approaches~\cite{vergara2013performance,vembu2012time} use more elaborate statistical modeling of the data, as described below. Moreover, they have been applied to the dataset described in Section~\ref{ch:dataset} and are thus of greater interest to us. In this section we briefly describe these algorithms~\cite{vergara2013performance,vembu2012time}, which are used in Section~\ref{ch:Exper} as performance benchmarks.

The study~\cite{vembu2012time} models the time series at the output of the sensors, denoted by $X(t)$, using an auto-regressive linear model of order $p$ ~\cite{cuturi2011autoregressive,lutkepohl2007newTimeSeries}

\begin{equation*}
X(t) = C + \sum \limits_{i=1}^{p} A_i X(t-i) + e_t
\end{equation*}
where $e_{t}$ is a stochastic noise term, and $C,A_{1},\ldots,A_{p}$ are parameters determined from the given observations.

 The parameters represent the time series $C,A_{1},\ldots,A_{p}$  in the sense that the parameters can predict its value at time $t$ from past observations, allowing for an error $e_{t}$. Once these parameters have been estimated, we map each time series $X(t)$ to its feature vector by  $(C,A_{1},\ldots,A_{p})$ and use these vectors as the input for the classification step. The classification step in~\cite{vembu2012time} uses a kernel SVM with a Gaussian kernel function. The idea behind kernel SVMs is to use a function $k(x,x')$ that measures the similarity between the features corresponding to each pair of instances in the dataset. The most commonly used function is the Gaussian $k(x,x')=\exp(-\gamma\|x-x'\|^2),$  where  $\gamma$ hyper-parameter, which is usually learned by cross-validation. It can be shown that the kernel SVM algorithm is equivalent to mapping each feature vector to some high (possibly infinite) dimensional space, followed by linear partitioning of that space. Passing the points through a kernel function gives rise to non-linear decision boundaries, which cannot exist in linear classification. The results reported in~\cite{vembu2012time} pertain to a reduced dataset and simplified clarification/prediction, which used only four gases and source location only is classified to be either the right or left sides of wind tunnel.

The study~\cite{vergara2013performance} is of particular interest to us since the conditions were exactly the same as in our study, namely classification of $10$ gases with similar training scenarios. The features space is the maximum of the normalized response of the sensors of each board.
The classification scheme of~\cite{vergara2013performance} is based on an inhibitory SVM classifier with a Gaussian kernel, whose main concept is to train a classier $f_j$ for each possible label $j=1,\...,L$ where function  $f_j$  constitutes the distance between the correct label and the most offending incorrect answer. For a more detailed description, see~\cite{huerta2012inhibition}.

\section{Mathematical background}\label{ch:Math}

In this section we present the theoretical background for our approach. In Section ~\ref{ch:Wav_graph} we present a generalization of the critically sampled Haar wavelet transform on graphs~\cite{gavish2010multiscale}  to an over-complete representation. This transform allows us to exploit correlations between different sensing boards and plays a critical role in our feature extraction process. In Section~\ref{ch:scat}, we review the scattering convolution network, introduced by Mallat, which is one of building blocks of our method. Finally, in Section~\ref{sec:RF}, we provide some details on the Random Forest algorithm.

\subsection{Redundant wavelet decomposition on graph}\label{ch:Wav_graph}

Let $G=(V,E)$ be an undirected and unweighed graph, where $V=\{v_i\}_{i=1}^N$  is the vertex set and $E \subset V \times V$ is the edge set. The multiscale folder-decomposition ~\cite{gavish2010multiscale} $\mathscr{V}=\{\mathcal{V}_j\}_{j=0}^{J}$ is collection of vertex set partitions where $\mathcal{V}^0=V$ is at resolution zero and for $j>0$, $\mathcal{V}^j=\{\upsilon^j_i\}_{i=1}^{n_j}, n_j=2^{-j}N$,  with

\begin{equation*}
\upsilon^j_i:=\upsilon^{j-1}_{\alpha_i} \cup \upsilon^{j-1}_{\beta_i}.
\end{equation*}
That is, the $i$-folder $\upsilon^j_i$ at level $j$ is obtained by grouping two folders $\upsilon^{j-1}_{\alpha_i}$  and $\upsilon^{j-1}_{\beta_i}$ from previous scale $j-1$.

A function $X:G \to \mathbb{R}^{n}$ on a graph $G$ is defined by mapping each vertex $v \in V$ to $X(v) \in \mathbb{R}^{n}.$

The dot-product of two functions $X_1,X_2: G \to \mathbb{R}$, is defined by

\begin{equation}\label{fr:dotProdG}
\langle X_1,X_2 \rangle_{G} := \sum \limits_{v \in G} X_1(v) \cdot X_2(v).
\end{equation}

We now define an Haar wavelet ortho-basis transform over a folder decomposition $\mathscr{V}$. The Haar wavelet is a function on the graph $G$ which is defined for each folder $\upsilon^j_i$ by
\begin{equation}\label{fr:HaarOrthoBasis}
\varphi_{j,i}:=\mathbbm{1}_{v^{j-1}_{\alpha_i}} - \mathbbm{1}_{v^{j-1}_{\beta_i}},
\end{equation}
where  $\mathbbm{1}_{W}$ is the indicator function on the set $W\subset V$ of the graph $G$, given by

$$
\mathbbm{1}_{W}(v)=
\left\{
  \begin{array}{l l}
    1 \text{ if } v \in W\\
    0 \text{ otherwise}
  \end{array} \right.
$$
According to~\cite{gavish2010multiscale}, the set of functions $\{\varphi_{j,i}\}$
is an orthogonal system defined on the graph $G=(V,E)$.
The set
\begin{equation}\label{fr:HaarBasis}
\Phi_{\mathscr{V}}:=\{\varphi_{j,i}\cup \mathbbm{1}_{v^j_{i}}\}_{i=0}^{n_j}  \text{, such that } n_j=2^{-j}N \text{ and }0\leq j \leq J
\end{equation
defines  redundant wavelets on the graph $G.$
Note that on each level $j$, there are $2^{-j}N$  folders and that the maximum scale of a folder-decomposition is  $K_{\max}=\lfloor \log_2 N\rfloor$. The application of the Haar transform to a signal $X$ defined on a graph $G$ proceeds as follows. For first scale $k=$1, the Haar coefficients are
\begin{equation}
 \left\{
  \begin{array}{l l}
    X^1(i,0) = \langle \mathbbm{1}_{v_i^1} , X \rangle_{G}, \\
    X^1(i,1)= \langle \varphi_{i,1}  , X \rangle_{G}.  \
  \end{array} \right.
 \end{equation}\label{fr:HaarZerroCoeef}

Then, for $k=2,...,K_{\max}$,
\begin{equation}\label{fr:HaarCoeef} \left\{
  \begin{array}{l l}
    X^k(i,2j) = \langle \mathbbm{1}_{v_i^k}, X^{k-1}(i,j) \rangle_G, \\
    X^k(i,2j+1)= \langle \varphi_{i,k} , X^{k-1}(i,j) \rangle_G.
  \end{array} \right.
  \end{equation}
where $i=1\...2^{-k}d.$
Observe that when $X(t)$, $t \in {\mathbb{T}}$, is in fact a time series signal over the graph, then we use the notation $X^k(t,i,j)$ for the time dependent multiscale Haar coefficients.

Next, we introduce a generalization of the redundant wavelet transform~\cite{beylkin1992representation,shensa1992discrete} to functions defined on a graph. This approach is effective for signal processing and pattern recognition problems such as image denoising~\cite{ram2012redundant} and speech recognition~\cite{tohidypour2012new}. Our approach is to build a redundant scheme of wavelet decomposition by overlapping the folders at each level. Figure~\ref{fg:RedundantHaar} illustrates a multiscale decomposition with folder overlapping, which leads to an over-complete wavelet representation of a signal over a graph domain. It is easy to demonstrate that the system of redundant Haar-functions $\Phi_{\mathscr{V}}$ defined in (\ref{fr:HaarBasis}), where ${\mathscr{V}}$ is a folder decomposition with overlapping, is not necessarily an orthogonal basis. However, redundant Haar functions provide more mutual information between vertices.

\begin{figure}
  \centering
  \includegraphics[scale=0.5]{RedundantHaar.eps}\\
  \caption{Multiscale folder-decomposition with overlapping.}\label{fg:RedundantHaar}
\end{figure}

For our olfaction problem, an overlapped folder decomposition using the neighborhood relationships of the boards. In Figure~\ref{fg:Lift_boards}, we can see the Haar wavelet lifting scheme, where overlapping is applied on center boards of the line position. Due to the wind direction in the tunnel, the air wind tunnel generates a diffusion chemical analyte. Therefore, centralized sensing boards contain more informative time series data~\cite{vergara2013performance}.
 \begin{figure}
  \centering
  \includegraphics[scale=0.125]{BoardLifting.eps}\\
  \caption{Multiscale folder-decomposition with overlapping Board Position}\label{fg:Lift_boards}
\end{figure}

\subsection{Scattering convolution network on the graph time series}\label{ch:scat}
In this section we review scattering convolution networks and their adaptations for a graph time series.
In~\ref{subch:DeepLearn}, we review fundamentals of convolutional networks and deep learning. In \ref{subch:Scat}, we review the wavelet-based scattering network for time series. In \ref{subch:ScatStabl}, we emphasize the stability properties of the wavelet approach, which constitute an advantage over the Fourier power spectrum. In \ref{subch:ScatMoments}, we define feature extraction and classification methods based on scattering.

\subsubsection{Fundamentals of Convolution Networks and Deep Learning} \label{subch:DeepLearn}
Deep learning (DL) architectures~\cite{bengio2009learning}  are neural network-based algorithms modeled to mimic the functionality of the human nervous system. The DL methods have been successfully applied in a variety of pattern recognition problems such as computer vision (face recognition, image classification, and annotation), natural language processing, speech recognition, and audio representations signals.
Learning is achieved through hierarchical feature extraction of the observed data. The main idea is to apply nonlinear processing on the data that flows between the layers. In fact, hierarchical features are fitted weight parameters of neural networks (free parameters of each unit), which are calculated as the results of one complex optimization process.

The convolutional neural network (CNN)~\cite{lecun1998gradient} is one the most popular deep learning architectures for pattern recognitions problems for grid-based signals such as time series, images, and videos, such as ImageNet~\cite{krizhevsky2012imagenet}. The CNN consists of units that use overlapping patches of input signals to apply neurons~\cite{lecun2010convolutional}. In other words, the CNN learns a hierarchical net of convolutions (filters) of signals.

Mallat introduced a mathematical class of deep convolution networks \cite{mallat2012group}, which is called \emph{`Scattering Convolution Networks'}. The scattering convolution networks is unsupervised CNN architecture for grid-based signals that are obtained by cascading wavelet transforms and modulus pooling operators with the average of the amplitude of iterated wavelet coefficients. The scattering-based feature extraction is translation invariant and Lipschitz continuous to deformations~\cite{bruna2013invariant}.

When trying to apply DL/CNN techniques in problems such as ours, a difficulty arises since the geometric configuration of the sensors is not necessarily regular, as in image processing where the pixels are well aligned on a uniform grid. Thus, in this work, the concept of the uniform grid is replaced by a graph structure.

\subsubsection{Scattering convolution network} \label{subch:Scat}
In this work, the Scattering network is applied separately to each time series of the type $X(t)=X^k(t,i,j)$ computed by the graph analysis of Section \ref{ch:Wav_graph}. A scattering convolution network is obtained based on a cascade of wavelet convolution and modulus operators with smoothing operator (low-pass filter)  \cite{mallat2012group,bruna4104multiscale,bruna2013intermittent}. Let $\psi(t)$ be a complex wavelet, whose real and imaginary parts are orthogonal and have the same  $L_2$-norm, and $\cap_{{\mathbb{R}}} \psi(t)dt=0$ with $|\psi(t)|=O\left( (1+t^2)^{-1} \right)$  with  dyadic dilations:
$$
\psi_j(t)=2^{-j}\psi(2^{-j}t), \quad \forall j \in {\mathbb{Z}}.
$$
The wavelet transform of time series $X(t)$  at scale $2^j$

$$
X \star \psi_j(t)=\cap_{{\mathbb{R}}} X(u)\psi_j(t-u)du.
$$
We calculate the absolute value of the complex value coefficient
$$
U_1[j]X(t)=|X\star \psi_j(u)|.
$$
For each sequence of indices  $\bar{p}=(j_1,\...,j_m):j_1<j_2<\...j_m$ the \emph{order-$m$ scattering propagator} $U[p]$ is defined by:
$$
U_m[\bar{p}]X(t)=U[j_m]\...U[j_1]=||\...||X\star\psi_{j_1}|\star \psi_{j_2}|\...|\star \psi_{j_m}(t)|
$$

The \emph{windowed scattering} $S_{m,J}[\bar{p}]$ (scattering coefficients of order-$m$) is defined by
\begin{equation}\label{fr:ScatCoeff}
S_{m,J}[\bar{p}]X(t)=U[\bar{p}]X\star \phi_J=||\...||X\star\psi_{j_1}|\star \psi_{j_2}|\...|\star \psi_{j_m}(t)|\star \phi_J
\end{equation}
where $j_1<\...<j_m<J$ and $\phi_J(t)=2^{-J}\phi(2^{-J}t)$ is low-pass filter
with $2^J$ scale and $\cap \phi(t)dt=1$.

The  scattering operator $S_J X(t)$ aggregates all scattering coefficients with order until layer $M$
\begin{equation}\label{fr:ScatAggregate}
S_J X(t)=\left( S_{m,J}X(t) \right)_{0\leq m\leq M}
\end{equation}
where $S_{0,J}X(t)= X \star \phi_J.$

The iterated procedure of scattering convolution network for a time series $X(t)$ is illustrated in Figure \ref{fg:Scat_timeSer}.
\begin{figure}
  \centering
 \includegraphics[scale=0.125]{Scat_timeSer.eps}\\
  \caption{Two levels scattering convolution network of time series $X(t)$ with scattering propagators $U[j_1]$ and $U[j_1,j_2]$ } \label{fg:Scat_timeSer}
\end{figure}

For most types of the signals, such as audio, images, biomedical signals, and finance time series, it is sufficient to compute the scattering coefficient of layers 0,1 and 2 (M=2).
$$
S_{2,J}X=
\begin{pmatrix}
S_{0,J}X\\
S_{1,J}\left[j_1\right]X\\
S_{2,J}\left[j_1,j_2\right]X
\end{pmatrix} _{j_1,j_2 \in \Lambda}
$$

\subsubsection{Scattering Deformation stability}\label{subch:ScatStabl}

The efficiency of a scattering representation comes from its invariance to local translations due to convolutions with $\phi_j$ and from its ability to linearize deformations, that is, its stability to time-warping.

The Fourier transform is unstable to deformation because dilating a sinusoidal wave yields a new sinusoidal wave of different frequency that is orthogonal to the original one~\cite{talmon2014manifold}. Let us define the deformation operator $D_{\tau}$  of signal $X(t)$

$$
D_{\tau}X(t)=X(t-\tau(t))
$$
where $\tau(t)$ non-constant deformation term.
As proven in~\cite[Theorem 2.12]{mallat2012group}, the scattering transform $S_J$ of a signal $X$ with compact support is Lipschitz continuous under action of deformation operator $D_{\tau}$.

\begin{equation}\label{fr:Scat_stable}
\|S_J(D_{\tau}X) -S_J(X)\|\leq C M\|X\| \left(2^{-J}|\tau|_{\infty} +|\nabla\tau|_{\infty}  \right)
\end{equation}
where  $|\tau|_{\infty}=\sup_t|\tau(t)|$ and $|\nabla\tau|_{\infty}=\sup_t|\nabla\tau(t)|<1.$

This property guarantees stability of signals $X(t)$. It is clear that the deformation error is small if the scaling factor $J$ is $2^J\gg|\tau|_{\infty}$
and the signal $X(t)$ is smooth in $L_1$ meaning. In other words, the scattering metric satisfies invariance to local transformations and deformations.

\subsubsection{Scattering moments}\label{subch:ScatMoments}
As noted in Section \ref{subch:ScatStabl} above, the scattering network is stable under small deformation. It can therefore  be used as an effective feature space for many kinds of classification and regression problems.

State-of-the-art results of the scattering approach have been obtained for handwritten digit recognition and texture classification \cite{bruna2010classification} compared to  convolution neural networks (CNN)~\cite{lecun2010convolutional,ranzato2007unsupervised} and dictionary learning (DL)~\cite{mairal2012task}.

Scattering moments are defined as expected values over time of scattering coefficients, for each path  $\bar{p}=(j_1,\...,j_m):j_1<j_2<\...j_m$
$$
\bar{S}[\bar{p}]X = {\mathbb{E}} \left( S_{m,J}[\bar{p}]X \right) = {\mathbb{E}} \left(U_m[\bar{p}] X \right)
$$
For finite time series signal $|X(t)| = N$
$$
\bar{S}[\bar{p}]X = \frac{1}{N}\sum \limits_{t=1}^{N} ||X\star\psi_{j_1}|\star \psi_{j_2}|\...|\star \psi_{j_m}(t)|
$$
According to \cite{anden2013deep,bruna2010classification}, the standard way to build feature space based on scattering moments for finite time series signal $|X| = N$ is
\begin{equation} \label{fr:scat_moments}
\Phi: X \to \begin{pmatrix}
\bar{S}[j_1] X\\
\bar{S}[j_1,j_2] X
\end{pmatrix} _{j_1,j_2 \in \Lambda}
\end{equation}
where a scaling set
$ \Lambda = \{ (j_1,j_2) : 1\leq j_1 = 2^{z_1/Q_1}  \leq  N  \text{ and } 1 \leq j_2=2^{z_2/Q_2} \leq j_1 \: z_1, z_2 \in {\mathbb{Z}} \}$
defines a filter bank of scattering transform, such that  are number wavelets per octave of the first and the second layer.

Scattering moments have been used as features space\cite{anden2013deep} for time series classification problems: musical genre classification (GITZAN) and phone segment classification. In these kinds of signals, the best state-of-the-art results were obtained by an SVM classifier with Gaussian kernel.

\subsection{Random Forest}\label{sec:RF}
Random forest (RF)~\cite{breiman2001random} is a popular ensemble learning method for classification and regression. At training time, a diverse set of decision trees is constructed using randomization techniques. Their output is then averaged to overcome the potential bias of each tree. In some implementations, the decision trees are pruned to reduce variance. RF easily allows a  parallel architecture to be implemented in applications for testing and training scenarios. RF has been implemented in many recent applications for classification and regression problems~\cite{verikas2011mining}. In the current study, we found RF to be an effective classifier for our machine olfaction problem. Due to the fact that the RF approach is based on an ensemble of decision trees, the feature set consists of variables from different domains, including  categorical and continuous variables. This allows us to add a set of static parameters, such as airflow velocity and operating temperature, to the feature space of time series.

Let a set of vectors $X^{tr}=\{x_1\...x_n : x_i \in {\mathbb{R}}^d\}$
be a training set consisting of $n$ samples, each a $d$-dimension feature vector with associated response values $Y^{tr}=\{y_1,y_2,\...,y_n\}$ that can be categorical variables in a classification problem or continuous variables in a regression problem. The training data can be interpreted as samples governed by an unknown distribution from $X \times Y$.  The goal of the RF classifier is to learn (approximate) the oracle function $F:X \to Y$, such that learned function ̂$\hat{f}$ is the  best approximation of $F$ with respect to an error metric (loss function). In the tree bagging approach, for each $k=1\...N$, the algorithm randomly selects samples with replacement from the training set  $\{X^{tr}_k,Y^{tr}_k\}$ in order to build decision trees $\{T_k\}$. The final prediction function $\hat{f}$, is an average of the predictions of the decision tree $\{T_k\}$, or by weighted voting of the ensemble.

\section{Main algorithm}\label{ch:MainAlgo}
In this section, we present the main algorithm developed in this paper, in which all the elements introduced in the previous sections are brought together. The feature construction is performed using the {\bf Scattering Time Series on Graph (STSG)}.The computed feature vectors, along with the response variables, serve as the training set for an RF algorithm.

The STSG algorithm combines the Haar scattering transform on the graph and standard wavelet-based scattering net described previously.
The Haar scattering network architecture for the signal defined on graphs (not time series) was described in~\cite{chen2014unsupervisedHaarScatGraphs}. The scattering architecture is obtained by cascading multiscale Haar wavelet transform defined on an embedded subset (folders construction~\cite{gavish2010multiscale}), but without redundant wavelet decomposition.

$X$ denotes a multivariate time series defined on an unweighted graph domain $G=(V,E)$, with $\dim(V)=d$ and finite time domain ${\mathbb{T}}$

$$
X: G \times {\mathbb{T}}  \to {\mathbb{R}}.
$$

The scattering time series on graphs $X(n,t)$ defined by:

\begin{equation}
	S_{J,\ell,k}[p_{\ell},\hat{j}]X=S_{J,\ell}[p_{\ell}]X^k(t,i,\hat{j}).
\end{equation}\label{fr:Scat_TimeGraph}
$S_{J,\ell}[p_{\ell}]$ consists of two actions. The first of these two actions is the calculation of $k$-levels at redundant Haar wavelet coefficients (\ref{fr:HaarCoeef}) with respect to the selected folder decomposition $\mathscr{V}$. The architecture of folder decomposition $\mathscr{V}$ can be derived from our understanding of the geometry of the sensing platform. The second action is the calculation of $\ell$-levels scattering coefficients of time series such that $p_{\ell}=j_1,\...j_l$, are their scaling paths (\ref{fr:ScatCoeff}).

The feature space of signal $X(n,t)$ consists of STSG moments, that is, the expected value scattering time series on graphs transform $S_{J,\ell,k}[p_{\ell},\hat{j}]$ over time.

 \begin{equation}\label{fr:Scat_TimeGraphFr}
 \Phi: X \to \bar{S}_{\mathcal{V}} =
  \mathbb{E}_{t} \left[S_{J,\ell,k}[p_{\ell},\hat{j}]X(t,.)\right].
 \end{equation}

Our machine learning task is as follows: Given a set of training signals $X^{tr},Y^{tr}$, such that each instance is a time series defined on the same graph $G$, the learning algorithm must seek a prediction function $f: X \to Y$.

Equipped with our training set, we calculate STSG moments  $\bar{S}_{\mathcal{V}}(x_i)$, $x_i \in X^{tr}$ (see \ref{fr:Scat_TimeGraphFr}). The second step is applying dimensionality reduction of the STSG moments to $d$ principal components
$$
\Phi: x_i \to \Phi^{d}(\bar{S}_{\mathcal{V}}(x_i)), \forall x_i \in X^{tr}.
$$
The dimensionality reduction increases stability. The scattering domain lies on a low dimensional manifold~\cite{bruna2010classification}.

In some of the classification scenarios, we have additional information about each instance $x_i$, which it is added to the feature space.  Since our learning algorithm is based on an RF classifier, with trees constructed over bagging of the training set, each iteration $k$ begins with random sampling $\{X_k,Y_k\}$ from the full training set, followed by feature mapping, $\bar{S}_{\mathcal{V}}(x_i)$, $x_i \in X_k$.

\begin{algorithm}
\caption{Training Random Forest of Scattering Graph Net}\label{algo:clasif}
\begin{algorithmic}[1]
\Procedure{TrainRF}{$X^{tr},Y^{tr}$}

\For{$k\gets 1, T$} \Comment{Loop for ensemble trees}

\State{$(X_k,Y_k) \gets \textbf{Bagging}(X^{tr},Y^{tr})$} \LONGCOMMENT{\textbf{Bagging}: Random sampling with replacement}

\State $\bar{S}_{\mathcal{V}}(x_i) \gets \textbf{STSG}(x_i), x_i \in X_k$ \LONGCOMMENT{Computing STSG moments}

\State $(\Phi^d_k,P_k^d) \gets PC_d\{\bar{S}_{\mathcal{V}}(x_i), x_i \in X_k \})$ \LONGCOMMENT{\textbf{Features space}: Dimensionality reduction to $d$ principal components}

\State $f_k \gets \textbf{TreeGrow}(\Phi^d_k,Y_k)$ \LONGCOMMENT{Growing a decision tree with feature bagging}

\State $F \gets F +f_k$ \LONGCOMMENT{Ensemble tree building}
\EndFor
\State $F \gets \frac{1}{N}F$\LONGCOMMENT{Final classifier uniform normalization}

\State \textbf{return} $F$\LONGCOMMENT{Output: Final tree ensemble and learned PCA transform}
\EndProcedure
\end{algorithmic}
\end{algorithm}

\section{Results} \label{ch:Exper}
In this section we demonstrate the application of our approach to machine olfaction problems.
In Section \ref{subch:features} , we  review the raw feature space. In Section \ref{subch:scenarios} we define classification scenarios.
 Finally, in Sections \ref{ch:Clasif}, \ref{ch:CO_consentr} and \ref{ch:SourceLoc}, we compare our results with prior-art scattering techniques and state-of-the-art machine olfaction techniques described in Section ~\ref{sec:sateOFart}.

\subsection{Feature Space}\label{subch:features}

To the STSG features (see Section~\ref{ch:MainAlgo}), we add some of the pre-established conditioning parameters described in Section~\ref{ch:dataset}

\begin{itemize}
  \item \emph{Heater voltage} $V_H \in \{4.0V, 4.5V, 5.0 V, 5.5 V, 6 V \}$,
  \item \emph{Airflow velocity} measured in rotation per minute: $rmp \in \{1500, 3900, 550\}$,
  \item \emph{Nominal  concentration} measured  in  parts-per-million  by  volume(ppmv).
\end{itemize}

The features of ambient temperature and relative humidity are not included since they have very low variance.  Our experiments clearly show that adding these static features significantly improves the performance of the prediction algorithm.

In addition, the location of each time series is known, namely, the position and board number with respect to the chemical source. The location of the sensing board is

\begin{equation} \label{fr:loc}
F_{loc}=\{X_{pos} \times X_{bord}\},
\end{equation}
where
\begin{itemize}
  \item $X_{pos} \in \{0.25, 0.5, 0.98, 1.18, 1.40, 1.45 \}$,
  \item $X_{board}= 0.13:0.13:1.2$.
\end{itemize}

\subsection{Classification scenarios}\label{subch:scenarios}
Our results are for complex scenarios of gas classification and source detection. These scenarios are motivated by proposals of the Department of Defense (DoD) for the development of chemo-sensing solutions and standards for early warning and protection of military forces against potential chemical and biological attacks (see \cite{guideHazards},\cite{chemHazards}). Specifically, the scenarios are as follows:

\begin{itemize}
  \item Gas classification problem (10 labels),
  \item Gas concentration prediction for $CO$ only (binary classification problem: $1,000$ ppm and $4,00$ ppm only),
  \item Source localization problem – prediction of source location with respect to local coordinates in the wind tunnel.
\end{itemize}

In all the above scenarios, we applied learning with/without static features and with/without source location information.

Based on prior studies, we constructed two ways to aggregate raw features:
\begin{enumerate}
  \item \textbf{Board Column}\cite{vergara2013performance}  aggregating all nine boards from than same position to a $72$-dimension time series,
  \item \textbf{Single Board}\cite{vembu2012time} using only the $8$-dimension time series from a single board.
\end{enumerate}
The `Single Board' scenario is significantly more difficult because fewer features are used.
We compared our results with frequency-domain features based on short-time Fourier transformation (\ref{fr:Fourier}) and simple features based first two statistical moments of each time series of sensor response.\\
In other words, we  compared three groups of features:
\begin{enumerate}
  \item STSG features (\ref{fr:Scat_TimeGraphFr}),
  \item Statistical moments,

  \item Fourier Power spectrum(STFT $L_2$-moments) (\ref{fr:Fourier}).
\end{enumerate}
For these three groups of features, we applied the same RF classifier with $200$ grown trees. The error rate was calculated as expected value and variance of 5-folder cross validation.

\subsection{Gas classification results}\label{ch:Clasif}

 Table \ref{tbl:Classif_byPos} shows the classification performance of Board Column suitable for the conditions  \cite{vergara2013performance}. The performance of the random forest classifier for all three groups of features is significantly better than the mean results of a previous study \cite{vergara2013performance}. As shown, using our STSG features provides significantely better performance in the more complex scenarios.

\begin{table}[h]
	\centering
\caption{\textbf{Classification performance of \emph{`Board Column'} scenario}.Testing error rate RF models with proposed STSG-algorithm, Short Time Fourier Transform(STFT) and statistical moments features, with/withput static features and with/withput location information  }
\begin{tabular}{c|c||c|c|c|c}
\hline\hline
\begin{tabular}[c]{@{}c@{}}Static\\ feature\end{tabular} & Location & STFT           & \begin{tabular}[c]{@{}c@{}}Statistical\\ moments\end{tabular} & STSG                 & \begin{tabular}[c]{@{}c@{}}\cite{vergara2013performance}\\ \end{tabular}\\[0.5ex]  \hline\hline
\multirow{2}{*}{FALSE}                                   & FALSE    & 2.52\%($\pm$0.09\%) & 0.574\%($\pm$0.07\%)                                             & {\bf 0.47\%($\pm$0.07\%)} & \multirow{4}{*}{7.89\%}                                           \\
                                                         & TRUE     & 2.41\%($\pm$0.16\%) & 0.58\%($\pm$0.06\%)                                              & {\bf 0.48\%($\pm$0.03\%)} &                                                                   \\
\multirow{2}{*}{TRUE}                                    & FALSE    & 1.26\%($\pm$0.15\%) & {\bf 0.09\%($\pm$0.07\%)}                                        & 0.26\%($\pm$0.08\%)       &                                                                   \\
                                                         & TRUE     & 1.31\%($\pm$0.25\%) & {\bf 0.09\%($\pm$0.02\%)}                                        & 0.24\%($\pm$0.07\%)       &                                                                   \\ \hline\hline
\end{tabular}
\label{tbl:Classif_byPos}
\end{table}

A comparison in the \emph{'Single Board'} scenario is presented in Table \ref{tbl:Classif}.

\begin{table}[h]
\centering
\caption{\textbf{Classification performance of \emph{`Single Board'} scenario}. Testing error rate RF models with STSG, STFT and statistic moments features with/without static and location information }
\begin{tabular}{c|c||c|c|c}
\hline\hline
\begin{tabular}[c]{@{}c@{}}Static \\ features\end{tabular} & Location & STFT            & \begin{tabular}[c]{@{}c@{}}Statistic\\  moments\end{tabular} & STSG            \\[1ex]  \hline\hline
\multirow{2}{*}{FALSE}                                     & FALSE    & 19.68\%($\pm$0.37\%) & 24.96\%($\pm$0.15\%)                                              & \textbf{18.11\%($\pm$0.11\%)} \\[0.5ex]
                                                           & TRUE     & 14.74\%($\pm$0.20\%) & 18.20\%($\pm$0.11\%)                                              & \textbf{13.63\%($\pm$00.20\%)} \\[0.5ex]
\multirow{2}{*}{TRUE}                                      & FALSE    & 3.15\%($\pm$0.12\%)  & 4.78\%($\pm$0.08\%)                                               & \textbf{2.80\%($\pm$00.08\%)}  \\[0.5ex]
                                                           & TRUE     & 2.30\%($\pm$0.07\%)  & 3.15\%($\pm$0.05\%)                                               & \textbf{2.14\%($\pm$0.03\%) } \\[0.5ex] \hline\hline
\end{tabular}
\label{tbl:Classif}
\end{table}
Figure \ref{fg:Learn_rate_classif} shows error bars on learning curves of validation error and out-of-bag error over the number of grown classification trees in the random forest ensemble for each kind of features.

 \begin{figure}
  \centering
  \includegraphics[scale=0.5]{Learing_clasif_rate.eps}\\
  \caption{Learning curves of the validation error and out-of-bag errors of classification performance of Single Board scenario scenario.}\label{fg:Learn_rate_classif}
\end{figure}
 The STSG learning curves have the fastest decay rate. This proves the learnability of the proposed method.

\subsection{Detection of $CO$-concentration }\label{ch:CO_consentr}
We now address a binary classification problem, in which the goal is to determine concentration of the $CO$ substance (carbon monoxide). Note that the given dataset includes every chemical substance in the experiment in only one nominal concentration, except for the carbon monoxide, which was collected in two different concentrations of $1,000$ ppm and $4,000$ ppm. We built a subset with binary labeling from a given dataset that contained only experiments with $CO.$ We then applied our RF classifier for compared feature space.
Table~\ref{tbl:ClassifCO} presents the classification performance of the Single Board scenario of the binary classification problem. Note that,  in this scenario, we cannot compare results with static features because these features have been simulated with different unique static parameters.

\begin{table}[h]
\centering
\caption{\textbf{Classification performance of  \emph{'Single Board} scenario for $CO$ concentration}. Testing error rate of RF classifier of proposed STSG, STFT, and statistical moments features space with/without location information .}
\label{tbl:ClassifCO}
\begin{tabular}{c||c|c|c}
\hline\hline
Location & STFT & \begin{tabular}[c]{@{}c@{}}Statistic \\ moments\end{tabular} & STSG \\ \hline\hline
FALSE & 0.30\%($\pm$0.05\%) & 1.00\%($\pm$0.16\%) & \textbf{0.24\%($\pm$0.07\%)} \\% \hline
TRUE & 0.26\%($\pm$0.06\%) & 0.61\%($\pm$0.10\%) & \textbf{0.20\%($\pm$0.05\%)} \\ \hline\hline
\end{tabular}
\end{table}
The proposed STSG feature space clearly demonstrates performance that is superior to that of other methods.. We did not apply this problem ’Board Column’, because ’Board Column’ scenario for all gases has perfect performance yet.

\subsection{Source localization }\label{ch:SourceLoc}
The source localization scenario is a learning of regression model for high accuracy detection of MOX-sensors location with respect to the gas substance. Tables \ref{tbl:LocPred_pos} and \ref{tbl:LocPred} show the rate error of the compared features space. The two types of features aggregation, ’Board Column’ and ’Single Board’ scenarios, were used respectively.
The error rate is the  $L_2$-norm error in meters between sensor board location and predicted location.
Note that if we use feature aggregation by Single Board, the classifier predicts the 2D location of the sensing board $F_{loc}$ (\ref{fr:loc}). When feature aggregation by Board Column is used, the classifier predicts distance of line position only $X_{pos}$.

\begin{table}[h]
\centering
\caption{\textbf{Location prediction performance of \emph{'Board Column'}}: $L_2$-norm error in meters with respect to source location}
\label{tbl:LocPred_pos}
\begin{tabular}{c||c|c|c}
\hline\hline
\begin{tabular}[c]{@{}c@{}}Static\\ Features\end{tabular} & STFT & \begin{tabular}[c]{@{}c@{}}Statistic\\ moments\end{tabular} & STSG \\ \hline\hline
FALSE & 0.02155($\pm$0.00063) & 0.00745($\pm$0.00028) & \textbf{0.00477($\pm$0.00060)} \\\hline
TRUE & 0.02025($\pm$0.00076) & 0.00773($\pm$0.00032) & \textbf{0.00481($\pm$0.00043)}\\ \hline\hline
\end{tabular}
\end{table}

\begin{table}[h]
		\centering
\caption{\textbf{Location prediction performance of \emph{'Single Board} scenarios}: $L_2$-Norm error in meters with respect to source location}
\begin{tabular}{c||c|c|c}
\hline\hline
\multicolumn{1}{c||}{\begin{tabular}[c]{@{}c@{}}Static \\ features\end{tabular}} & \multicolumn{1}{c}{STFT} & \multicolumn{1}{|c|}{\begin{tabular}[c]{@{}c@{}}Statistic \\ moments\end{tabular}} & \multicolumn{1}{c}{STSG} \\[1ex]  \hline\hline
FALSE                                                                          & 0.174($\pm$0.001)             & 0.288($\pm$0.001)                                                                     & {\bf 0.162($\pm$0.001)}       \\ \hline
TRUE                                                                           & 0.161($\pm$0.010)             & 0.280($\pm$0.006)                                                                     & {\bf 0.150($\pm$0.008)}
\\ \hline\hline
\end{tabular}
\label{tbl:LocPred}
\end{table}

\section{Conclusions}
In the current study we presented a novel methodology that can be used for pattern recognition problems in which the signals are collected from an array of possibly non-uniform ensemble of sensors. We defined a novel scattering transform for multidimensional time series on the graph (STSG), which used redundant wavelet graph decomposition. In this study, we focused on machine olfaction problems in which the dataset was obtained from a chemical gas sensor array in a turbulent wind tunnel. We applied our methodology to three machine learning problems: classification of $10$ different gases at various concentrations, detection of $CO$ consecration, and chemical substance localization.\\
The next step of this research is to develop the proposed method formachine olfaction in cases involving  turbulent gas mixtures~\cite{fonollosa2014chemical}.
Another interesting research direction would be to add \textit{soft supervised learning to scattering net}. Recent advances in classical deep learning, such as CNN, have demonstrated high accuracy performance. These methods are based on the learning of convolutional filters. However, it must be noted that deep learning architecture typically requires a large training dataset and the learning process is computationally intensive. Recall that the learning process when using scattering networks is faster since we use pre-designed filters. One could then consider `soft' learning variant of the scattering network, where standard wavelet filters are used as initial filters and then optimization is applied to only a few wavelet parameters.


\section{Introduction}
\label{INTRODUCTION}

Fission tracks are dislocated zones caused by nuclear fragments released
in spontaneous fission of uranium-238. Information about fission tracks
can be related to geologic events, as mineral crystallization age, geologic
fault zones and thermal events\cite{WAGNER1992}.

Tracks crossing a polished mineral surface can be etched and visualized
under an optical microscope, and its selection is based on the following
relatively simple criteria\cite{FLEISCHER1964,PRICE2005}:

\begin{itemize}
	\item Fission tracks form straight line defects of a limited length ($< 20 \mu m$);
	\item They exhibit no preferred orientation and disappear after suitable heating.
\end{itemize}

Manual identification of fission tracks has some practical problems, such as
variation due to observer-observation efficiency. Also, Gleadow
\textit{et al.}\cite{GLEADOW2009} list some problems in discrimination of
fission tracks from non-track defects as polishing scratches and resolving
multiple track overlaps and small tracks amongst a similarly sized
background of surface defects.

An automatic identification of fission track could solve this problem and
improve the speed of track counting. Image processing can be used to
automatize such task; however, separation of nontrivial images is one
of the most difficult tasks in image processing\cite{GONZALEZ2008}. Several
commercial and free softwares are available for this purpose. Nonetheless,
these softwares are meant to be used in specific images\cite{USAJ2011}.

\subsection{Proposed methodology}
\label{PROPMET}

In this paper we propose an automatic method based on starlet wavelets,
in order to segment fission tracks in images of natural minerals obtained
by optical microscopy. Commonly used objective lenses (dry or oil immersion
type) have total magnification up to 1500 times. In combination with a
reflected-transmitted light system, it is possible to analyze fission
tracks\cite{WAGNER1992,GLEADOW1986}.

The proposed approach consists of applying starlet wavelets in a sample
image to obtain its detail decomposition levels. Based on information
retrieval (precision, recall and accuracy) and Matthews correlation
coefficient (MCC), the segmentation level that better represents fission
tracks of the original image is automatically chosen.

This technique is an improvement of a recent study aimed at segmentation
of scanning electron microscopy images\cite{SIQUEIRA2014}. An application
of this method is the separation of fission tracks in images of natural
minerals, such as volcanic glasses, apatite, zircon, muscovite, epidote,
among others.

In this paper, the proposed methodology is applied to segment fission
tracks in photomicrographs of epidote crystals. Results presented in this
study will be used as a basis to develop an open source software capable
of extracting fission tracks from images of natural mineral samples in
order to establish the age of the material using the fission track dating
method\cite{WAGNER1992}. A prototype of this software, containing the
algorithms used in this study, is available for download on this journal
website (see \ref{SUPPLMAT}).

The remainder of this paper follows. Section \ref{MATANDMET} introduces
the material used in this study and starlet wavelets, as well as an
overview of evaluation and automatization methods. Next, Section \ref{RESULTS}
presents the results from this method application in test photomicrographs.
Moreover, the method performance is discussed. In the following, Section
\ref{CONCLUSION} presents the final considerations about this study.
Finally, \ref{SUPPLMAT} explains where to obtain the cited
algorithms and how to use them.

\section{Material and Methods}
\label{MATANDMET}

\subsection{Epidote crystals}

Epidote is a mineral with monoclinic crystal structure and general formula
$Ca_{2}(Al,Fe)_{3}Si_{3}O_{12}(OH)$\cite{ITO1950}. According to Poli and
Schmidt\cite{POLI2004}, it is possible to have epidote formation at
temperatures of $500 \sim 700^{\circ}C$ (pressure range of 0.2 to 0.6 GPa),
and also at $720 \sim 760 ^{\circ}C$ (pressure range of 1.6 to 3 GPa). Their
formation is given by different means. One of them is deuteric action,
during the late phase of magmatic crystallization stage, by regional
metamorphism and hydrothermal activity, \textit{i.e.} percolation of solutions
which chemically react with the rock through fractures, often in
temperatures between $300$ and $500^{\circ}C$\cite{BAR1974}.

In order to evaluate the proposed methods, we used a data set consisting
of 45 images. These images were obtained from epidote phenocrystals using
a Carl ZEISS optical microscope with Axiocam Imager.M1m system, nominal
magnification factor of $1000X$ (dry) and transmitted light.

\subsection{Starlet transform}

Starlet wavelet transform is an isotropic redundant wavelet based on the
algorithm ``\textit{à trous}'' (with holes)\cite{HOLSCHNEIDER1990,SHENSA1992}.
The construction of this wavelet is given by its scale and wavelet functions,
respectively $\phi_{1D}$ and $\psi_{1D}$ (Eqs. \ref{STARSCALE} and \ref{STARWAVE},
\cite{STARCK2006,STARCK2010}), where $\phi_{1D}$ is the third order B-spline
($B_{3}$-spline).

\begin{eqnarray}
\phi_{1D}(t) = \frac{1}{12}\left(|t-2|^3-4|t-1|^3+6|t|^3-4|t+1|^3+|t+2|^3\right) \label{STARSCALE} \\
\frac{1}{2}\psi_{1D}\left(\frac{t}{2}\right) = \phi(t)-\frac{1}{2}\phi\left(\frac{t}{2}\right) \label{STARWAVE}
\end{eqnarray}

An extension to two dimensions is achieved by a tensor product (Eq. \ref{STAR2D}),
\begin{eqnarray}
\phi(t_1,t_2) = \phi_{1D}(t_1)\phi_{1D}(t_2) \nonumber\\
\frac{1}{4}\psi\left(\frac{t_1}{2},\frac{t_2}{2}\right) = \phi(t_1,t_2)-\frac{1}{4}\phi\left(\frac{t_1}{2},\frac{t_2}{2}\right) \label{STAR2D}
\end{eqnarray}

These wavelets were successfully employed in analysis of astronomical
\cite{STARCK2006,STARCK2010,STARCK2011} and biological \cite{GENOVESIO2003}
images, being suitable to evaluate images that contains isotropic objects.
Isotropic transforms retrieve only one detail set per level instead of
several detail sets (\textit{e.g.} Daubechies wavelets have horizontal, vertical
and diagonal detail levels), facilitating the interpretation of results. The
properties of this wavelet (isotropy, redundancy, translation-invariance)
make it a good alternative in image processing and pattern recognition.

Similarly to Eq. \ref{STAR2D}, the pair of filters $(h, g)$ related to
this wavelet is (Eq. \ref{FILTERPAIR}, \cite{STARCK2010})
\begin{eqnarray}
h_{1D}[k] = [\begin{array}{ccccc} 1 & 4 & 6 & 4 & 1\end{array}]/16, k = -2,...,2 \nonumber\label{H1D}\\
h[k,l] = h_{1D}[k]h_{1D}[l] \nonumber\\
g[k,l] = \delta[k,l]-h[k,l] \label{FILTERPAIR}
\end{eqnarray}
where $\delta$ is defined as $\delta[0,0] = 1$, $\delta[k,l] = 0$ for $[k,l] \neq 0$.
From Eqs. \ref{STAR2D} and \ref{FILTERPAIR}, detail wavelet coefficients
are obtained from the difference between the current and previous resolutions.

Starlet transform application is given by a convolution between an input
image $c_{0}$ and the finite impulse response (FIR) filter derived from
$\phi$ (Eq. \ref{H2D} \cite{STARCK2010}),

\begin{center}
\begin{equation}
h = \left[\begin{array}{ccccc} \frac{1}{256} & \frac{1}{64} & \frac{3}{128} & \frac{1}{64} & \frac{1}{256} \\
\frac{1}{64} & \frac{1}{16} & \frac{3}{32} & \frac{1}{16} & \frac{1}{64} \\
\frac{3}{128} & \frac{3}{32} & \frac{9}{64} & \frac{3}{32} & \frac{3}{128} \\
\frac{1}{64} & \frac{1}{16} & \frac{3}{32} & \frac{1}{16} & \frac{1}{64} \\
\frac{1}{256} & \frac{1}{64} & \frac{3}{128} & \frac{1}{64} & \frac{1}{256} \end{array}\right]
\label{H2D}
\end{equation}
\end{center}

This convolution results in a set of smoothing coefficients which correspond
to the first decomposition level, $c_{1}$. Detail wavelet coefficients of
the first decomposition level are obtained from the difference
$w_{1} = c_{0} - c_{1}$.

Let $L$ be the last resolution level. Therefore, resolution levels can be
calculated by:
\begin{eqnarray*}
c_{j} = c_{j-1} * h, \\
w_{j} = c_{j-1} - c_{j},
\end{eqnarray*}
\noindent with $j = 0, \ldots, L$, and $*$ the convolution operation. The set
$W = \{w_{1}, \ldots, w_{L}, c_{L}\}$
obtained by these operations is the starlet transform of the input image.

\subsection{Evaluation of the results}
\label{SECEVAL}

In order to evaluate the proposed methodology, we employed precision,
recall and accuracy\cite{OLSON2008,WANG2013}. These values are based on
the concepts of true positives (TP), true negatives (TN), false positives
(FP) and false negatives (FN).

Fission tracks in an image sample are represented in a ground truth (GT)
image: black represents the background, whereas white represents
fission tracks in this image. Comparing an input image and its ground
truth, TP, TN, FP and FN values could be established as:

\begin{itemize}
\item TP: pixels correctly labeled as fission tracks.
\item FP: pixels incorrectly labeled as fission tracks.
\item FN: pixels incorrectly labeled as background.
\item TN: pixels correctly labeled as background.
\end{itemize}

Based on these considerations, precision (retrieved pixels that are
relevant), recall (relevant pixels that were retrieved) and accuracy
(proportion of true retrieved results) are defined:

\begin{eqnarray}
precision &=& \frac{TP}{TP+FP} \nonumber \\
recall &=& \frac{TP}{TP+FN} \nonumber \\
accuracy &=& \frac{TP+TN}{TP+TN+FP+FN} \nonumber
\end{eqnarray}

\subsubsection{Method automatization}

In order to establish the optimal level for method application, Matthews
correlation coefficient (MCC, \cite{MATTHEWS1975}) is used. MCC uses TP,
TN, FP and FN, and may offer an evaluation of the segmentation correctness:

\begin{center}
\begin{equation}
MCC = \frac{TP*TN-FP*FN}{\sqrt{(TP+FN)(TP+FP)(TN+FP)(TN+FN)}} \label{MCC}
\end{equation}
\end{center}

This coefficient measures how variables tend to have the same sign and
magnitude, where $1$, zero and $-1$ indicates perfect, random and imperfect
predictions, respectively\cite{BALDI2000}.

Thereby, automatic retrieval of the optimal segmentation level is achieved by:
\begin{itemize}
\item applying the method for $L$ desired starlet decomposition levels;
\item comparing method segmentation results with the image GT and obtaining
TP, TN, FP and FN;
\item calculating MCC (Eq. \ref{MCC}) for each $L$.
\end{itemize}
As the optimal segmentation level is reached, values yielded by MCC become
higher. Then, the best segmentation level obtained by this method is the
one that returns the highest MCC value.

\subsection{Method overview}

The proposed automatic segmentation method is defined as follows:
\begin{itemize}
\item Starlet transform is applied in an input image $c_{0}$, resulting
in $L$ detail levels: $D_{1}, \cdots, D_{L}$, where $L$ is the last
desired resolution level.
\item First and second detail levels, $D_{1}$ and $D_{2}$, are ignored
due to the large amount of noise; third to $i$ detail levels are summed
($Ri = D_{3} + \cdots + D_{i}$), where $3 \leq i \leq L$. This is the
result of the method application related to starlet level $i$.
\item TP, TN, FP and FN are obtained comparing Ri with its GT (see
Section \ref{SECEVAL}). MCC (Eq. \ref{MCC}) is calculated using these
values.
\item Therefore, the optimal segmentation level is the one that has a
higher MCC between the $R_{L}$ levels of method application.
\end{itemize}

To apply this method, one can use the pseudocodes given in Algorithm
\ref{ALGORITHM1}. Also, the source code of these algorithms in
Matlab\footnote{Matlab is a numerical computing environment and a
programming language developed by MathWorks. A trial version could be
requested at \url{https://www.mathworks.com}.}/Octave\footnote{GNU Octave is
an open source high-level interpreted language intended primarily for
numerical computation. Download available freely at
\url{http://www.gnu.org/software/octave/download.html}.} programming language
is available (see \ref{SUPPLMAT}).

\begin{algorithm}
\SetAlgoNoEnd
\DontPrintSemicolon
\KwIn{\\
$\bullet$ A grayscale image, $c_{0}$.\\
$\bullet$ A ground truth image, $gtc_{0}$.\\
$\bullet$ Number of resolutions to be calculated, $L$.}
\KwOut{\\
$\bullet$ Detail coefficients from starlet transform, $w_{j}$.\\
$\bullet$ An optimal image that presents fission tracks contained in the
original image, $optft$. \\
$\bullet$ Matthews Correlation Coefficient (Eq. \ref{MCC}) between $gtc_{0}$ and
the algorithm result for each level $imgft_i$, $MCC_{i}$,
with $i = 3,...,L$.}
\textbf{mirroring(}$c_{0}$\textbf{);}
\For{$j \gets 1$ \textbf{to} $L$} {
	$h \gets$ \textbf{hgen(}j\textbf{);} \\
	$c_{j} \gets$ \textbf{convolution(}$c_{j-1},h$\textbf{);} \\
	$w_{j} \gets c_{j-1}-c_{j}$\textbf{;} \\
	\textbf{unmirroring(}$c_{j}$\textbf{);} \\
	\textbf{increment(}$j$\textbf{);} \\
}
\textbf{initialize} $sum$ \textbf{to} 0\textbf{;}

\For{$i \gets 3$ \textbf{to} $L$} {
	\For{$j \gets 3$ \textbf{to} $i$} {
		$sum \gets sum + w_{j};$ \\
		$imgft_{i} \gets sum;$
	}
	$MCC_{i}(imgft_{i}, gtc_{0})$; \\
}
$optft \gets \max(MCC_{i})$; \\
\Return{$w_{j},optft,MCC_{i}$}\;
\caption{Pseudocode for automatic determination of fission tracks in an
image, based on starlet algorithm application (adapted from \cite{SIQUEIRA2014,STARCK2011}).}
\label{ALGORITHM1}
\end{algorithm}

Starlet transform of $c_{0}$ is given by $W = \{w_{1}, \cdots, w_{L}, c_{L}\}$.
\textbf{hgen()} (Algorithm \ref{ALGORITHM2}), referenced on Algorithm
\ref{ALGORITHM1}, is applied when $j$ is incremented. For $j>1$, $h$ has $2^{j-1}$
zeros between its elements, characterizing the \textit{à trous} transform.

\begin{algorithm}
\SetAlgoNoEnd
\DontPrintSemicolon
\KwIn{\\
$\bullet$ $h_{1D}$ filter, given by Eq. \ref{H1D}.\\
$\bullet$ Current resolution level, $j$.}
\KwOut{\\
$\bullet$ Filter $h_{2D}$, h.}
\If{$j = 0$}{
	$h \gets h_{1D}$;\;
}
\Else{
	$M \gets$ \textbf{size(}$h_{1D},2$\textbf{)};
	\textbf{initialize} $k$ \textbf{to} $0$;\;
	\For{$i \gets 1$ \textbf{step} $2^{i-1}$ \textbf{to} $M+2^{i-1}*(M-1)$}{
		\textbf{increment(}$k$\textbf{)};\;
		$h(i) \gets h_{1D}(k)$;\;
	}
}
\textbf{initialize} $aux$ \textbf{to} 0;\;
$aux \gets sum(sum(h'*h))$;\;
$h \gets (h'*h)/aux$;\;
\Return{$h$}\;
\caption{\textbf{hgen:} $h$ filter generation and zero-inserting\cite{SIQUEIRA2014}.}
\label{ALGORITHM2}
\end{algorithm}

\section{Experimental results}
\label{RESULTS}

In order to present the proposed method results, six dataset images with
different size and fission track distribution are shown (Fig.
\ref{FIGSAMPLES}). The darker regions of these images correspond to
fission tracks in epidote surface.

\begin{figure*}[h]
\centering
\includegraphics[scale=0.4]{Fig1.pdf}
\caption{Photomicrographs of epidote phenocrystals. Fission tracks are
shown in sample surface as dark segments. Nominal magnification factor:
1000X (dry).}
\label{FIGSAMPLES}
\end{figure*}

The proposed method was applied in the test images with $L = 3$ to $L = 9$.
The optimal segmentation level was obtained from MCC, for each image.
Also, precision, recall and accuracy were obtained for each level (Fig.
\ref{FIGPRECRECACC}), in order to evaluate the method performance. For a
satisfactory segmentation degree, an optimal ratio between precision and
recall becomes necessary.

\begin{figure*}[h]
\centering
\includegraphics[scale=0.45]{Fig2.pdf}
\caption{Precision, recall and acuracy values obtained from Fig. \ref{FIGSAMPLES}.}
\label{FIGPRECRECACC}
\end{figure*}

One could see from Fig. \ref{FIGPRECRECACC} that accuracy and precision
increases until level $L = 7$. On the other hand, recall decreases as
$L$ increases.

Fig. \ref{FIGSAMPLES}(a) will be used to introduce the proposed method.
According to the first step, starlet transform detail levels are obtained
from the input image (Fig. \ref{FIGSTARLET}).

\begin{figure*}[htb]
\centering
\includegraphics[scale=0.4]{Fig3.pdf}
\caption{Starlet detail decomposition levels of Fig. \ref{FIGSAMPLES}(a).
$D = 1$ and $D = 2$ were disregarded in method application due to noise
amount. Although higher detail levels tend to aggregate, reducing
segmentation accuracy, these levels present better information about
fission tracks.}
\label{FIGSTARLET}
\end{figure*}

After starlet application, Algorithm \ref{ALGORITHM1} is applied seven
times, from $L = 3$ to $L = 9$. For example, method application for
$L = 6$ consists of:
\begin{itemize}
	\item disregard $D = 1$ (Fig. \ref{FIGSTARLET}(a)) and D = 2 (Fig. \ref{FIGSTARLET}(b)).
	\item sum $D = 3$ (Fig. \ref{FIGSTARLET}(c)) to $D = 6$ (Fig. \ref{FIGSTARLET}(f)): $\sum D_{i}, i = 3,\cdots, 6$.
\end{itemize}
Results of the proposed method are shown as binary images, where fission
tracks are represented by the white color and background by the black color
(Fig. \ref{FIGOUTPUT}).

\begin{figure*}[htb]
\centering
\includegraphics[scale=0.4]{Fig4.pdf}
\caption{Segmentation output of the proposed method. Different levels,
from $L = 3$ to $L = 9$ were considered.}
\label{FIGOUTPUT}
\end{figure*}

The next step is to determine the optimal segmentation level using MCC
(Table \ref{TABLE1}). Values presented here are given in
percentages in order to ease results comprehension. The optimal
segmentation level for images of Fig. \ref{FIGSAMPLES} is $L = 7$, according
to MCC; then results obtained with $L = 7$ using the proposed method
were compared to ground truth (GT).

GT images obtained from Fig. \ref{FIGSAMPLES} are used to evaluate the
method performance. These images were obtained manually by a specialist
using GIMP\footnote{Available freely at http://www.gimp.org/downloads/.},
an open source graphics software. TP pixels are shown as green, FN
pixels as blue and FP pixels as red, to facilitate visualization of
differences between the images (Fig. \ref{FIGCOMPARISON}).

\begin{table*}[htb]
\caption{MCC obtained from method application for levels $L = 3$ to
$L = 9$. $L = 7$ was the optimal segmentation level for Fig. \ref{FIGSAMPLES}
images, according to MCC values.}
{\small
\begin{center}
\begin{tabular}{lccccccc}
\toprule
\bf{MCC (\%)} & \bf{L = 3} & \bf{L = 4} & \bf{L = 5} & \bf{L = 6} & \bf{L = 7} & \bf{L = 8} & \bf{L = 9}\\
\midrule
\bf{Fig. 1 (a)} & 19.32996 & 25.30251 & 35.37489 & 49.55525 & 57.85511 & 56.06013 & 52.62559\\
\bf{Fig. 1 (b)} & 24.22276 & 32.30995 & 45.43835 & 61.52304 & 66.14956 & 55.69191 & 43.78881\\
\bf{Fig. 1 (c)} & 29.07748 & 40.25502 & 56.84173 & 74.15157 & 79.22872 & 71.01649 & 49.15956\\
\bf{Fig. 1 (d)} & 16.09310 & 20.57131 & 27.80563 & 36.43833 & 41.70653 & 39.22664 & 30.58268\\
\bf{Fig. 1 (e)} & 24.60903 & 32.38622 & 45.38681 & 63.17591 & 67.75544 & 60.42723 & 41.72874\\
\bf{Fig. 1 (f)} & 22.82558 & 28.93141 & 40.72572 & 62.05960 & 69.84571 & 64.87040 & 60.59735\\
\bottomrule
\end{tabular}
\end{center}
}
\label{TABLE1}
\end{table*}

\begin{figure*}[htb]
\centering
\includegraphics[scale=0.5]{Fig5.pdf}
\caption{First column: Ground truths for Fig. \ref{FIGSAMPLES} images.
Second column: proposed method output for the optimal decomposition level.
Third column: a comparison between GT and method output for optimal level
(Green: TP pixels; blue: FN pixels; red: FP pixels).}
\label{FIGCOMPARISON}
\end{figure*}

Most fission tracks presented in GT images were located using level
$L = 7$. While precision and recall values vary, accuracy is higher than
$89\%$ for Fig. \ref{FIGSAMPLES} images (Table \ref{TABLE2}). Comparing
these results with near levels, $L = 6$ and $L = 8$, one can see that
accuracy method for Fig. \ref{FIGSAMPLES} using $L = 7$ is higher than
with $L = 6$, but lower than $L = 8$ for Fig. \ref{FIGSAMPLES}(a),(d),(e)
and (f). However, the difference is small for these cases, less than
$1.5\%$. Furthermore, level 8 recall values are smaller than recall for
$L = 7$ (difference between 14 and $22.5\%$).

\begin{table*}[htb]
\caption{Precision, recall and accuracy values for starlet application
levels 6, 7 and 8.}
{\small
\begin{center}
\begin{tabular}{lccccccccc}
\toprule
& \multicolumn{3}{c}{\bf{Precision (\%)}} & \multicolumn{3}{c}{\bf{Recall (\%)}} & \multicolumn{3}{c}{\bf{Accuracy(\%)}}\\
\midrule
& \bf{L = 6} & \bf{L = 7} & \bf{L = 8} & \bf{L = 6} & \bf{L = 7} & \bf{L = 8} & \bf{L = 6} & \bf{L = 7} & \bf{L = 8}\\
\midrule
\bf{Fig. 1 (a)} & 30.047 & 40.674 & 45.556 & 91.918 & 88.860 & 74.529 & 90.899 & 94.232 & 95.306 \\
\bf{Fig. 1 (b)} & 46.044 & 55.792 & 56.656 & 91.219 & 84.953 & 61.034 & 92.163 & 94.420 & 94.197 \\
\bf{Fig. 1 (c)} & 74.977 & 94.820 & 95.545 & 80.685 & 70.104 & 56.958 & 93.731 & 95.428 & 93.805 \\
\bf{Fig. 1 (d)} & 16.186 & 21.400 & 22.303 & 97.641 & 92.448 & 79.438 & 85.188 & 89.878 & 91.330 \\
\bf{Fig. 1 (e)} & 46.118 & 56.291 & 60.873 & 95.075 & 87.678 & 65.457 & 92.369 & 94.711 & 94.960 \\
\bf{Fig. 1 (f)} & 45.815 & 57.202 & 62.156 & 92.491 & 90.864 & 72.802 & 92.581 & 95.105 & 95.461 \\
\bottomrule
\end{tabular}
\end{center}
}
\label{TABLE2}
\end{table*}

Accurate results were obtained when images contained better visual state;
for example, images without grain fractures (as seen in Fig.
\ref{FIGSAMPLES}(d)). Also, textures in the background may lead to incorrect
segmentation, thus lowering method accuracy. The red agglomerated regions
in Fig. \ref{FIGOUTPUT} exhibit this issue.

\section{Conclusion}
\label{CONCLUSION}

In this study we present a fission track automatic segmentation method
for photomicrographs, based on starlet wavelets. This method uses starlet
decomposition detail levels to determine edges of objects in an input
image. Levels corresponding to noise are discarded and remain levels
are summed. The method presented here can help the user to determine fission
tracks in photomicrographs of mineral samples.

Automation is achieved using precision, recall and accuracy, together with
Matthews correlation coefficient. MCC has proved to be a satisfactory measure
to the method automation, representing a good balance between precision,
recall and accuracy.

An application of this method is given here, in epidote crystal images
obtained by optical microscopy. In this application, the proposed method
presents a high accuracy degree, even for challenging images.

Algorithms used in this study are available for download. In future studies,
this methodology will be used in images obtained by different materials,
in order to estimate related features. Also, from this algorithms, an
open source software aimed to analyze fission tracks in microscopical
images will be built.

\section*{Acknowledgements}
\label{ACKNOWLEDGEMENTS}

The authors would like to acknowledge the Brazilian foundations of research
assistance CNPq, CAPES and FAPESP. This research is supported by FAPESP
(Procs 2010/20496-2 and 2011/09438-3).