Image Segmentation and Time Series Clustering Based on Spatial and Temporal ARMA Processes

1. Introduction

During the past decades, image segmentation and edge detection have been two important and challenging topics. The main idea is to produce a partition of an image such that each category or region is homogeneous with respect to some measures. The processed image can be useful for posterior image processing treatments.

Spatial autoregressive moving average (ARMA) processes have been extensively used in several applications in image/signal processing. In particular, these models have been used for image segmentation, edge detection and image filtering. Image restoration algorithms based on robust estimation of a two-dimensional process have been developed (Kashyap & Eom 1988). Also the two-dimensional autoregressive model has been used to perform unsupervised texture segmentation (Cariou & Chehdi, 2008). Generalizations of the previous algorithms using the generalized M estimators to deal with the effect caused by additive contamination was also addressed (Allende et al., 2001). Later on, robust autocovariance (RA) estimators for two dimensional autoregresive (AR-2D) processes were introduced (Ojeda, 2002). Several theoretical contributions have been suggested in the literature, including the asymptotic properties of a nearly unstable sequence of stationary spatial autoregressive processes (Baran et al., 2004). Other contributions and applications of spatial ARMA processes have been considered in many publications (Basu & Reinsel, 1993, Bustos 2009a, Francos & Friendlaner1998, Guyon 1982, Ho 2011, Illig & Truong-Van 2006, Martin1996, Vallejos & Mardesic 2004).

A new approach to perform image segmentation based on the estimation of AR-2D processes has been recently suggested (Ojeda 2010). First an image is locally modeled using a spatial autoregressive model for the image intensity. Then the residual autoregressive image is computed. This resulting image possesses interesting texture features. The borders and edges are highlighted, suggesting that the algorithm can be used for border detection. Experimental results with real images clarify how the algorithm works in practice. A robust version of the algorithm was also proposed, to be used when the original image is contaminated with additive outliers. Applications in the context of image inpainting were also offered.

Another concern that has been pointed out in the context of spatial statistics is the development of coefficients to compare two spatial processes. Coefficients that take into account the spatial association between two processes have been proposed in the literature. (Tjostheim, 1978) suggested a nonparametric coefficient to assess the spatial association between two spatial variables. Later on, (Clifford et al. 1989) proposed an hypothesis testing procedure to study the spatial dependence between two spatial sequences. Rukhin & Vallejos (2008) studied asymptotic properties of the codispersion coefficient first introduced by Matheron(1965). The performance and impact of this coefficient to quantify the spatial association between two images is currently under study Ojeda et al. (2012). An adaptation of this coefficient to time series analysis was studied in Vallejos (2008).

In the context of clustering time series Chouakria & Nagabhushan (2007) proposed a distance measure that is a function of the codispersion coefficient. This measure includes the correlation behavior and the proximity of two time series. They proposed to combine these distances in a multiplicative way, introducing a tuning constant controlling the weight of each quantity in the final product. This makes the measure flexible to model sequences with different behaviors, comparing them in terms of both correlation and dissimilarity between the values of the series.

The structure of this chapter consist in two parts. In the first part we review some theoretical aspects of the spatial ARMA processes. Then the algorithm suggested by Ojeda(2010), its limitations and advantages are briefly described. In order to propose a more efficient algorithm new variants of this algorithm are suggested specially to address the problem of determining the most convenient (in terms of the quality of the segmentation) prediction window of unilateral AR-2D processes. The computation of the distance between the filtered images and the original one will be done by using the codispersion coefficient and other image quality measures (Wang and Bovik 2002). Examples with real images will highlight the features of the modified algorithm. In the second part, the codispersion coefficient previously used to measure the closeness between images is utilized in a distance measure to perform cluster analysis of time series. The distance measure introduced in Chouakria & Nagabhushan (2007) is generalized in the sense that considers an arbitrary lag $h$ that allows us to capture a higher serial correlation of two temporal or spatial sequences. Examples and numerical studies are presented to explore our proposal in several different scenarios. We explore the performance of hierarchical methods to classify correlated sequences when the proposed proximity measure is used, employing the Monte Carlo simulation. An application is discussed for time series measuring the Normalized Difference Vegetation Index (NDVI) in four locations of Argentina. The clusters formed using hierarchical classification techniques with the proposed distance measure preserve the geographical location where the series were obtained, providing information that is unavailable when using hierarchical methods with conventional distance measures.

2. Image Segmentation Through Estimation of Spatial ARMA Processes

2.1. The Spatial ARMA Processes

Spatial ARMA processes have been studied in the context of random fields indexed over ${\mathbb{Z}}^d$ , $d\ge 2$ , where ${\mathbb{Z}}^d$ is endowed with the usual partial order. That is, for $s=(s_1,s_2,\dots ,s_d)$ , $u=(u_1,u_2,\dots ,u_d)$ in ${\mathbb{Z}}^d$ , $s\le u$ if for $i=1,2,.\dots ,d$ , $s_i\le u_i$ For $a,b\in {\mathbb{Z}}^d,$ such that $a\le b$ and $a\ne b,$ we define $S[a,b]=\{x\in {\mathbb{Z}}^d|a\le x\le b\}$ and $S⟨a,b]=S[a,b]\backslash \{a\}.$

A random field ${(X_s)}_{s\in {\mathbb{Z}}^d}$ is said to be a spatial ARMA $$\left(p,q\right)$$ with parameters $p,q\in {\mathbb{Z}}^d$ if it is weakly stationary and satisfies the equation

$$X_s-{\sum }_{j\in S⟨0,p]}{\varphi }_j{X}_{s-j}={\epsilon }_t+{\sum }_{k\in S⟨0,q]}{\theta }_j{\epsilon }_{s-j},$$ E1

where ${({\varphi }_j)}_{j\in S⟨0,p]}$ and ${({\epsilon }_j)}_{k\in S⟨0,q]}$ , respectively, denote the autoregressive and moving average parameters with ${\varphi }_0={\theta }_0=1,$ and ${({\epsilon }_s)}_{s\in {\mathbb{Z}}^d}$ denotes a sequence of independent and identically distributed centered random variables with variance ${\sigma }^2$ . Notice that if $p=0,$ the sum over $S⟨0,p]$ is supposed to be zero, and the process is called a spatial moving average MA $$\left(q\right)$$ random field. Similarly, if $q=0$ , the process is called a spatial autoregressive AR $\left(p\right)$ random field. The ARMA random field is labeled as causal if it has the following unilateral representation.

$$X_s={\displaystyle \sum _{j\in S[0,\infty ]}{\psi }_j}{\epsilon }_{s-j}$$

with ${\displaystyle \sum _j|}{\varphi }_j|<\infty$ . Similar to the time series case, there are conditions on the (AR or MA) polynomials that ensure stationarity and invertibility, respectively. Let $\Phi (z)=1-{\sum }_{j\in S⟨0,p]}{\varphi }_j{z}^j$ and $$\Theta (z)=1-{\sum }_{j\in S⟨0,q]}{\theta }_j{z}^j,$$ where $z=(z_1,z_2,\dots ,z_d)$ and $z^j=z_1^{j_1}{z}_2^{j_2}\dots z_d^{j_d}.$ A sufficient condition for the random field to be causal is that the AR polynomial $\Phi (z)$ has no zeros in the closure of the open disc $D^d$ in ${ℂ}^d$ . For example, if $d=2,$ the process is causal if $\Phi (z_1,z_2)$ is not zero for any $z_1$ and $z_2$ that simultaneously satisfy $\left|z_1\right|<1$ and $\left|z_2\right|<1$ (Jain et al., 1999).

Applications of spatial ARMA processes have been developed, including analysis of yield trials in the context of incomplete block designs (Cullis & Glesson 1991, Grondona et al. 1996) and the study of spatial unilateral first-order ARMA model (Basu & Reinsel, 1993). Other theoretical extensions of time series and spatial ARMA models can be found in (Baran et al., 2004, Bustos et al., 2009b, Gaetan & Guyon 2010, Choi 2000, Genton & Koul 2008, Guo 1998, Vallejos and Garccía-Donato 2006).

2.2. An Image Segmentation Algorithm

In this section, we describe an image segmentation algorithm that is based on a previous fitting of spatial autoregressive models to an image. This fitted image is constructed by dividing the original image into squared sub-images (e.g., $8\times 8$ ) and then fitting a spatial autoregressive model to each sub-image (i.e., block). Then, we generate a sub-image from each local fitted model, preserving intensities on the boundary to smooth the edges between blocks. The final fitted image is yielded by putting together all generated sub-images.

Let $$Z=Z_{m,n}$$ , $0\le m\le M-1$ , $0\le n\le N-1$ , be the original image, and let $$\text{X}={\text{X}}_{m,n}$$ , $0\le m\le M-1,0\le n\le N-1,$ where for all $0\le m\le M-1$ , $0\le n\le N-1$ , $X_{m,n}=Z_{m,n}-\overline{Z},$ and $\overline{Z}$ is the mean of $$\text{Z}$$ . Let $4\le k\le \min (M,N)$ and consider the rearrange images

$$\begin{array}{l}\text{Z}={\text{Z}}_{m,n},\\ \text{X}={\text{X}}_{{}_{m,n}},\end{array}$$

where $0\le m\le M^{\prime }-1$ , $0\le n\le N^{\prime }-1$ , $M^{\prime }=\left[\frac{M-1}{k-1}\right](k-1)+1$ , $N^{\prime }=\left[\frac{N-1}{k-1}\right](k-1)+1$ . For all $i_b=1,\cdots ,\left[\frac{M-1}{k-1}\right]$ and for all $j_b=1,\cdots ,\left[\frac{N-1}{k-1}\right]$ the $(k-1)\times (k-1)$ block $\left(i_b,j_b\right)$ of the image $X$ is defined as

$$B_X\left(i_b,j_b\right)=X_{r,s},$$

where $(k-1)(i_b-1)+1\le r\le (k-1)i_b$ and $(k-1)(j_b-1)+1\le s\le (k-1)j_b$ . Then, the approximated image $\widehat{X}$ of $X$ is provided by Algorithm 1.

Algorithm 1.

For each block $B_X\left(i_b,j_b\right)$

1. Compute estimators ${\widehat{\varphi }}_1\left(i_b,j_b\right)$ , ${\widehat{\varphi }}_2\left(i_b,j_b\right)$ of ${\varphi }_1$ and ${\varphi }_2$ corresponding to the block $B_X\left(i_b,j_b\right)$ extended to:

$$B_X^{\prime }\left(i_b,j_b\right)=X_{r,s},$$

where $(k-1)(i_b-1)\le r\le (k-1)i_b$ ,$(k-1)(j_b-1)\le s\le (k-1)j_b.$

2. Let $\widehat{X}$ be defined on the block $B_X\left(i_b,j_b\right)$ by

$${\widehat{X}}_{r,s}={\widehat{\varphi }}_1\left(i_b,j_b\right)X_{r-1,s}+{\widehat{\varphi }}_2\left(i_b,j_b\right)X_{r,s-1}$$

where $(k-1)(i_b-1)+1\le r\le (k-1)i_b$ and $(k-1)(j_b-1)+1\le s\le (k-1)j_b.$

Then the approximated image $\widehat{Z}$ of the original image $Z$ is:

$${\widehat{Z}}_{m,n}={\widehat{X}}_{m,n}+\overline{Z},\text{0}\le m\le M^{\prime }-1,0\le n\le N^{\prime }-1.$$

The image segmentation algorithm we describe below is supported by a widely known notion in regression analysis. If a fitted image very well represents the patterns on the original image, then the residual image (i.e., the fitted image minus the observed image) will not contain useful information about the original patterns because the model already explains the features that are present in the original image. On the contrary, if the model does not well represent the patterns that are present in the original image, then the residual image will contain useful information that has not been explained by the model. Thus, to implement an algorithm based on these notions, we must characterize which patterns are present in the residual image when the fitted image is not a good representation of the original one, and we must develop a technique to produce a fitting that is satisfactory in terms of segmentation but not a very good estimation in that the residual image still contains valuable information. (Ojeda et al. 2010) investigated these concerns and, based on several numerical experiments with images, determined that the residual image associated with a good local fitting is in fact poor in terms of structure (i.e., it is very similar to a white noise). However, when the fitted image is poor in terms of estimation, the residual image is useful for highlighting the boundaries and edges of the original image. Moreover, a bad fitting is related to the size of the block (or window) used in Algorithm 1. The best performance is attained for the maximum block size, which would be the size of the original image. The image segmentation algorithm introduced by (Ojeda et al. 2010) can be summarized as follows.

Algorithm 2.

1. Use Algorithm 1 to generate an approximated image $\widehat{Z}$ of $Z.$

2. Compute the residual autoregressive image given by $Z-\widehat{Z}$

Example 1. We present examples with real images to illustrate the performance of Algorithms 1 and 2. These images were taken from the database http://sipi.usc.edu/database. Figure 1(a) shows an original image of size $512\times 512$ (aerial), and Figure 1(b) shows the image generated by Algorithm 1 when a moving window of size $512\times 512$ is used to define an AR-2D process on the plane. It is not possible to visualize the differences between the original and fitted images. However, the residual image (Fig 1(c)) shows patterns that the model is not able to capture. Basically, the AR-2D model does not capture the changes in the texture produced by lines, borders and object boundaries. These features are contained in the residual image produced by Algorithm 2 such that the good performance of Algorithm 2 is associated with a moderate fitting of the AR-2D model. Another image (peepers) was processed by Algorithm 2 to show the effect of the size of the moving window. Figure 2(b) shows the segmentation produced by Algorithm 2 using a moving window of size $128\times 128.$ Another segmentation with a moving of size $512\times 512$ is shown in Figure 2. In both cases, the segmentations highlight the borders and boundaries present in the original image.

2.3. Improving the Segmentation Algorithm

In all experiments carried out in (Ojeda et al., 2010) and (Quintana et al., 2011), Algorithm 1 was implemented using the same prediction window for the AR-2D process, which contains only two elements belonging to a strongly causal region on the plane. Here, we consider other prediction windows to observe the effect on the performance of Algorithm 2. A description

Figure 1.

Images generated by Algorithms 1 and 2.

Figure 2.

(b)-(c) Images generated by Algorithm 2 with prediction windows of $128\times 128$ and $512\times 512$ respectively.

of the most commonly used prediction windows in statistical image processing is in Bustos et al., (2009a). A brief description of the strongly causal prediction windows is given below.

For all $(m,n)\in {\mathbb{Z}}^2,$ a strongly causal region at $(m,n)$ is defined as

$$S(m,n)=\{(k,l)\in {\mathbb{Z}}^2:k\le m,l\le n)\}-\{(m,n)\}$$ E2

For a given $M\in \mathbb{N},$ a strongly causal prediction window is

$$W=\{(k,l)\in s(m,n):m-M\le k\le m,n-M\le l\le n\}.$$ E3

In particular, if $M=1$ , then a strongly causal prediction window containing three elements is

$$W_1=\{(k,l)\in s(m,n):m-1\le k\le m,n-1\le l\le n\}$$ E4

The set $W_1$ is shown in Figure 3 (b). Similarly, strongly causal prediction windows can be defined by considering not only the top left quadrant on the plane ${\mathbb{Z}}^2.$ The definition of such sets generates the prediction windows $W_2,W_3,$ and $W_4$ , as shown in Figure3(b). Algorithms 1 and 2 were implemented using the prediction windows $W_1,W_2,W_3,$ and $W_4,$ with two elements each (Figure 3(a)).

Figure 3.

Strongly causal prediction windows.

Visually, the best segmentation for the aerial image is yielded by the prediction window $W_2.$ The lines and edges are better highlighted in this segmentation (Figure 4(b)) than in the other segmentations. The dark regions are also stressed, which provides a more intense and brighter partition of the original features.

To gain insight on image quality measures, the fitted images produced by Algorithm 1 associated with the images shown in Figure 4(a) -(d) were compared aerially with the original image using three coefficients described in (Ojeda et al., 2012). These coefficients are briefly described below.

Consider two weakly stationary processes, ${(X_s)}_{s\in D}and{(Y_s)}_{s\in D}$ , $D\subset {\mathbb{Z}}^d.$ For a given $h\in D,$ the codispersion coefficient is defined as

$$\rho (h)=\frac{\gamma (h)}{\sqrt{V_X(h)V_Y(h)}},$$ E5

where $s,s+h\in D,\gamma (h)=\mathbb{E}[X(s+h)-X(s)][Y(s+h)-Y(s)],V_X(h)=\mathbb{E}{[X(s+h)-X(s)]}^2$ (and similarly for $${\text{V}}_Y\left(h\right)$$ ).

For $d=2,$ the sample codispersion coefficient is defined by

$$\widehat{\rho }(h)=\frac{{\sum }_{s\text{,}s+h\in D^{\text{'}}}{a}_s{b}_s}{\sqrt{{\widehat{V}}_X(h){\widehat{V}}_Y(h)}}$$ E6

with $s=(s_1,s_2),h=(h_1,h_2),D^{\prime }\subseteq D,\#(D^{\prime })<\infty$ , $$a_s=\text{X}\left(s_1+h_1,{\text{ s}}_2+h_2\right)-X\left(s_1,{\text{ s}}_2\right),$$ $$b_s=\text{Y}\left(s_1+h_1,{\text{ s}}_2+h_2\right)-Y\left(s_1,{\text{ s}}_2\right),$$ ${\widehat{V}}_X(h)={\displaystyle {\sum }_{s,s+h\in D^{\prime }}{a}_s^2},$ and $${\widehat{V}}_Y(h)={\displaystyle {\sum }_{s,s+h\in D^{\prime }}{b}_s^2}.$$

Figure 4.

a)-(d) Images generated by Algorithm 2 with prediction windows $W_1-W_4$ respectively.

The index $Q$ (Wang and Bovik, 2002) is

$$Q=\frac{4S_{XY}\overline{X}\overline{Y}}{(S_X^2+S_Y^2)[{\overline{X}}^2+{\overline{Y}}^2]}=\frac{S_{XY}}{S_X{S}_Y}·\frac{2\overline{X}\overline{Y}}{{\overline{X}}^2+{\overline{Y}}^2}·\frac{2S_X{S}_Y}{S_X^2+S_Y^2}=C·M·V,$$ E7

where $\overline{X}$ is the mean of ${(X_s)}_{s\in D}$ , $S_X$ is the standard deviation of ${(X_s)}_{s\in D},$ and $S_{XY}$ is the covariance between ${(X_s)}_{s\in D}$ and ${(Y_s)}_{s\in D}$ (and similarly for $\overline{Y}$ and $S_Y$ ). The quantity $C=S_{XY}/S_X{S}_Y$ models the linear correlation between ${(X_s)}_{s\in D}$ and ${(Y_s)}_{s\in D}$ , $M=2\overline{X}\overline{Y}/({\overline{X}}^2+{\overline{Y}}^2)$ measures the similarity between the sample means (luminance) of ${(X_s)}_{s\in D}$ and ${(Y_s)}_{s\in D}$ , and $V=2S_X{S}_Y/(S_X^2+S_Y^2)$ measures the similarity related to the contrast between the images. Coefficient $Q$ is defined as a function of the correlation coefficient; hence, it is able to capture only the linear association between ${(X_s)}_{s\in D}$ and ${(Y_s)}_{s\in D}$ It is unable to account for other types of relationships between these sequences, including the spatial association in a specific direction $h.$ Ojeda et al. (2012) suggested by the $CQ$ index, which is defined as:

$$CQ(h)=\widehat{\rho }(h)·M·V,$$ E8

where $M$ and $V$ are defined as in (7).

The correlation coefficient and the coefficients defined in (6), (7) and (8) were computed to compare the fitted images, which were generated with a prediction window with two elements and associated with the images shown in Figure 4(a) -(f), and the original images. The results are shown in Table 1. In all cases, the highest values of the image quality measures are attained for the image fitted using the prediction window $W_2.$ This means that the residual image shown in Figure 4 (b) is the best segmentation yielded by Algorithm 2. The same

Table 1.

Image quality measures between the fitted and original (aerial) images related to the residual images shown in Figure 4.

experiment was carried out for the image shown in Figure 2(a). Table 2 summarizes the values of the image quality coefficients for the fitted images generated by Algorithm 2 with prediction windows $W_1,W_2,W_3,$ and $W_4.$ In this case, the best performance is for the fitted image

Table 2.

Image quality measures between the fitted and original (peppers) images.

generated with prediction window $W_4.$ In general, the performance of Algorithm 2 depends on the choice of the prediction window. One way to choose the prediction window that yields the best segmentation is to maximize the association between the fitted and original images. Indeed, if we denote the original image by $Z$ and the fitted image generated by Algorithm 1 with the prediction window $W_i$ by ${\widehat{Z}}_{W_i},$ then the prediction window that produces the best segmentation can be obtained by finding the maximum value of one of the quality measures (6), (7) or (8) between $Z$ and ${\widehat{Z}}_{W_i}$ . This criterion is summarized in the following algorithm.

Algorithm 3.

1. Use Algorithm 1 to generate the approximated images ${\widehat{Z}}_{W_i}$ of $Z,$ $i=1,2,3,4.$

2. Compute an image quality index between $Z$ and ${\widehat{Z}}_{W_i}$ for all $i=1,2,3,4.$ Suppose that the maximum value for the image quality index is attained for ${\widehat{Z}}_{W_j},$ $1\le j\le 4.$ Then, the best fitted image is ${\widehat{Z}}_{W_j}.$

3. Compute the residual autoregressive image $Z-{\widehat{Z}}_{W_j}.$

3. Clustering Time series

3.1. Measuring Closeness and Association Between Time Series

Let $x=(x_1,x_2,\dots ,x_p)$ and $y=(y_1,y_2,\dots ,y_q)$ be two time series. There are several conventional distance measures between time series. For example, if $p=q=n$ , then the Euclidean distance between $x$ and $y$ is defined as $d_E(x,y)={\left({\displaystyle \sum _{i=1}^n{(x_i-y_i)}^2}\right)}^{1/2}.$ As is evident, $d_E$ ignores information about the dependence between $x$ and $y$ . The Minkowski distance is a generalization of the Euclidean distance, which is defined as

$$d_M(x,y)={\left({\displaystyle \sum _{i=1}^n{(x_i-y_i)}^q}\right)}^{1/q},$$ E9

where $q$ is a positive integer. The Fréchet distance was introduced to measure the proximity between continuous curves. Let $m$ be a natural number such that $m\le \min (p,q).$ Let $M$ be the set of all mappings $r$ between $x$ and $y$ such that $r$ is a sequence of $m$ pairs preserving the order

$$r=((x_{a_1},y_{b_1}),(x_{a_2},y_{b_2}),\dots ,(x_{a_m},y_{b_m})),$$

where $a_i\in \left\{1,2,\mathrm{...},p\right\}$ ,$b_j\in \left\{1,2,\mathrm{...}q\right\}$ with $a_1=1,b_1=1,a_m=p,b_m=q$ and for $i\in \{1,2,\dots ,m-1\}$ , $a_{i+1}=(a_i\text{or}{a}_i+1)$ and $b_{i+1}=(b_i\text{or}{b}_i+1).$ Note that $\left|r\right|={\max }_{i=1,2,\dots ,m}|x_{a_i}-y_{b_i}|$ is the mapping length representing the maximum span between two coupled observations. Thus, the Fréchet distance between the series $x$ and $y$ is given by

$$d_F(x,y)=\underset{{}_{r\in M}}{\min }\left|r\right|=\underset{{}_{r\in M}}{\min }\left(\underset{{}_{{}_{i=1,2,\dots ,m}}}{\max }|x_{a_i}-y_{b_i}|\right).$$ E10

Dynamic time warping (DTW) is a variant of the Fréchet distance that considers mapping length as the sum of the spans of all coupled observations. That is,

$$\left|r\right|={\displaystyle \sum _{i=1,2,\dots ,m}|}{x}_{a_i}-y_{b_i}|.$$

Dynamic time warping is then defined as

$$d_{DTW}(x,y)=\underset{{}_{r\in M}}{\min }\left|r\right|=\underset{{}_{r\in M}}{\min }{\displaystyle \sum _{i=1,2,\dots ,m}|}{x}_{a_i}-y_{b_i}|.$$ E11

The distances defined above are based on the proximity of the values $|x_{a_i}-y_{b_i}|.$ However, these distances disregard both the temporal dependence between the sequences $x$ and $y$ and the correlation structure of each sequence.

Several distance measures that are functions of the correlation between two sequences ( $Cor(x,y)$ ) have been suggested. For example, (Golay et al.,1998) proposed

$$d_{cc}(x,y)={\left(\frac{1-Cor(x,y)}{1+Cor(x,y)}\right)}^{\beta }and{d}_{cc}^2(x,y)=2(1-Cor(x,y)),$$

where $\beta$ is a parameter related to the fuzzy $c$ -means classification algorithm (Macqueen, 1967). However, none of these measures takes into account the serial association between the sequences because the correlation coefficient is a crude measure of association. This approach requires the study of coefficients that are capable of capturing the spatial or serial correlation between two sequences.

3.2. The Codispersion Coefficient for Time Series

Consider two weakly stationary processes, $X=\{X_s:s\in D\subset \mathbb{Z}\}$ and $Y=\{Y_s:s\in D\subset \mathbb{Z}\},$ and let $x$ and $y$ be realizations of these processes as in Section 3.1. For $d=1,$ the estimator (6) becomes

$$\widehat{\rho }(h)=\frac{{\sum }_{t\in N(h)}(x_{t+h}-x_t)(y_{t+h}-y_t)}{\sqrt{{\displaystyle \sum _{t\in N(h)}{(x_{t+h}-x_t)}^2}{\displaystyle \sum _{t\in N(h)}{(y_{t+h}-y_t)}^2}}}$$ E12

where $N(h)=\{t:t+h\in D\},$ $N=|N(h)|$ is the cardinality of $N(h),$ and sequences $x$ and $y$ are realizations of processes $x$ and $y$ , respectively. The coefficient $\widehat{\rho }(h)$ is called the comovement coefficient when $h=1.$ Although $\widehat{\rho }(h)$ is not the correlation coefficient, the codispersion coefficient shares a number of its standard properties. For example, $\widehat{\rho }(h)$ is translation invariant, positive homogeneous, symmetric in its arguments, positive definite for a sequence and lagged versions of itself, and interpretable as the cosine of the angle between the vectors formed by the first difference of the sampled series. As in the case of classic correlation, a codispersion coefficient of $+1$ indicates that the compared functions or processes are rescaled or retranslated versions of one another. Similarly, a profile matched with its reflection across the time axis yields a codispersion of $-1$ . The value $\widehat{\rho }(h)=0$ expresses that there is no monotonicity between $x$ and $y$ and that their growth rates are stochastically linearly independent. More details can be found in (Rukhin & Vallejos, 2008 , Vallejos, 2008).

3.3. Dissimilarity Index for Time Series

This coefficient involves a distance measure and a correlation-type measure that addresses both the correlation behavior and the proximity of two time series. The dissimilarity index depends on similarity behaviors, which should be specified in advance. The suggested dissimilarity index $D(x,y,h)$ for the realizations $x$ and $y$ and $d=1$ is given by

$$D(x,y,h)=f(\widehat{\rho }(h))·d(x,y),$$ E13

where $f$ is an adaptive tuning function, and $d(x,y)$ is one of the conventional distances described in Section 3.1 that summarizes the closeness of sequences $x$ and $y$ . There are many possible ways to choose a function $f$ . Here, we follow the guidelines given in (Chouakria & Nagabhushan, 2007), according to which $f$ is considered an exponential adaptive tuning function given by

$$f_k(t)=\frac{2}{1+exp(kt)},$$ E14

where $k\ge 0$ modulates the contributions of the proximity with respect to values and behavior. For example, when $|\widehat{\rho }(1)|$ is large and $k=2,$ the proximity with respect to behavior contributes 76.2% to $D$ . The flexibility of $D$ allows us to choose $k$ such that for highly dependent sequences, the correlation structure can have a large weight in (13).

Note that (13) is a generalization of the dissimilarity index introduced in Chouakria & Nagabhushan, (2007). The dissimilarity index (13) can capture high-order serial correlations between the sequences because the distance lag $h$ is arbitrary, while Chouakria's index only captures the first-order correlation.

The dependence of (13) on $h$ is crucial, and in some specific cases, $h$ can be chosen using an optimal criterion. For example, for two AR(1) processes with parameters ${\varphi }_1$ and ${\varphi }_2,$ respectively, and a correlation structure between the errors (Rukhin & Vallejos,2008 ), it is possible to find ${\varphi }_1,{\varphi }_2$ such that $Var(\widehat{\rho }(1))Var(\widehat{\rho }(2)).$ In other words, for those processes in which the asymptotic variance of the codispersion coefficient is known, we suggest setting the value of $\widehat{h}$ to produce the minimum variance. That is,

$$\widehat{h}=\underset{h}{\mathrm{argmin}}[Var(\widehat{\rho }(h))].$$ E15

When the variance of the codispersion coefficient is difficult to compute, resampling methods can be use to estimate the variance of the sample codispersion coefficient (Politis & Romano, 1994, Vallejos, 2008).

In the next section, we present two simulation examples to illustrate the capabilities of the hierarchical methods using the distance measure (13) under the tuning function given by (14). All else being constant, the clusters produced using traditional distances are usually different from those yielded using the distance measure (13).

3.4. Simulations

In this example, we simulate observations from six first-order autoregressive models to illustrate the clustering produced by hierarchical methods when the sequences exhibit serial correlation. To generate the series, we consider the following models.

$$X_t^i={\varphi }_i{X}_{t-1}^i+{\epsilon }_t^i,i=1,2,\mathrm{...6},$$

where $\forall i=1,2,\mathrm{...},6,X^i={\left\{X_t^i\right\}}_{t\in \mathbb{Z}}$ define the $i-$ model, and the sequence ${\epsilon }^i={\left\{{\epsilon }_t^i\right\}}_{t\in \mathbb{Z}}$ is zero-mean white noise. Note that the sequences ${\epsilon }^1$ and ${\epsilon }^2$ have the covariance structure

$$Cov({\epsilon }_t^1,{\epsilon }_s^2)=\{\begin{array}{l}\rho \sigma \tau ,s=t,\\ 0,\text{otherwise,}\end{array}$$

with ${\sigma }^2={\tau }^2=1,$ and $\rho =0.9$ . The same covariance structure is assumed for ${\epsilon }_t^1$ and ${\epsilon }_{{}^t}^3,$ with ${\sigma }^2={\tau }^2=1$ and $\rho =\mathrm{0.7.}$ $${\epsilon }^i,i=4,5,6,$$ are assumed to be white noise processes with variance 1 and are uncorrelated with all other error sequences. If $i,j\le 3$ , the correlation structure between processes $X^i$ and $X^j$ , $i\ne j$ , is not null due to the correlation structure between ${\epsilon }^i$ and ${\epsilon }^j$ . Instead, if $i,j\ge 4$ , $i\ne j$ , the correlation structure between processes $X^i$ and $X^j$ vanishes.

Two hundred observations were generated from each model for ${\varphi }_1=-0.5,$ ${\varphi }_2=-0.3,$ ${\varphi }_3=-0.7,$ ${\varphi }_4=0.1,$ ${\varphi }_5=-0.9,$and ${\varphi }_6=-\mathrm{0.2.}$ The goal was to perform time series clustering with the Euclidean distance and (13). Hierarchical methods with complete linkage using both measures were implemented to evaluate whether the methods are capable of capturing the correlation structure between the sequences described above. We used the distance measure (13) under the tuning function (14) for $h=1,$ and $k=3$ . That is, the correlation structure contributes 90.5% to $D,$ whereas the proximity with respect to values contributes 9.5%.

Figure 5.

a) Time series clustering using the Euclidean distance, (b) Time series clustering using $D$ .

In Figure 5, we see that the dendrogram obtained using hierarchical methods with the Euclidean distance does not recognize the correlation structure between $X^1$ and $X^3.$ In this case, sequences $X^1,X^2,X^4,$ and $X^5$ are pulled together before sequence $X^3.$ However, hierarchical methods using (13) yield the expected results, combining sequences $X^1,X^2,\text{and}{X}^3$ before the rest of the series.

To obtain better insight into the classification process using the proposed distance measure (13), we carried out a second simulation study that involves clustering measures based on other distances (but using the same setup). Observations from models 1-6 were generated using Gaussian white noise sequences for the errors, thereby preserving the same correlation structure used in the first study. The goal was to explore the ability of the distance measure (13) to group strongly correlated series first. A total of 1000 runs were considered for this

Table 3.

Percentage of correct clustering of the correlated series 1,2 and 3.

experiment, and 200 observations were generated in each run. We used measure (13) under the tuning function (14) for $h=1,2$ and $k=1,2,3,4.$ We counted the number of times that the hierarchical algorithm with complete linkage was able to pull together series 1, 2 and 3 before connecting them with other sequences. The traditional distances described in Section 2 were also implemented. After the 1000 simulation runs were finished, the percentage of times that the algorithm was able to recognize the correlated series was recorded. The results of the experiment are summarized in Table 3 and 4.

Table 4.

Percentage of correct clustering of the correlated series (1,2, and 3). $D_{DTW}(h,k)$ is distance measure (13) with the DTW distance. $D_M(h,k)$ and $D_F(h,k)$ denote distance measure (10) with Minkowski and Frechet distances respectively.

Note from Table 3 that the traditional distance measures failed to group the correlated sequences, with the exception of the Minkowski distance, which correctly grouped the correlated series 99% of the time. The hierarchical algorithm that uses the distance measure (13) has a higher percentage of well-clustered correlated sequences than the same algorithm using the traditional distance measures described in Section 2 (see Table 4). The percentage of correct clusters increased in all cases with the distance measure (13), suggesting that hierarchical algorithms can be improved by including coefficients of association that consider high-order cross-correlation.

3.5. The NDVI Data Set

In this section, we consider time series from four different locations in Argentina. The data set consists of 15 monthly NDVI series measured during a period of 19 years (i.e., January 1982-December 2000). The observed values correspond to a transformation to the interval $[0,255]$ of the original NDVI series, which commonly resides in the interval $[-1,1].$ The data were collected by a NOAA sensor at 1 km resolution and provided by the Comisión Nacional de Actividades Espaciales (CONAE) in Córdoba, Argentina. Fifteen time series were collected from the following: the Amazon region in the northeast of Argentina (1, 2, 3), the Patagonia region in the south of Argentina (4, 5, 6, 7), the Pampean region in the center of Argentina (8, 9, 10, 11) and the Pre-Andean zone of Argentina (12, 13, 14, 15). The time series are shown in Figure 6.

Figure 6.

Fifteen NDVI series collected from four different regions in South America.

We can observe a variety of different patterns in Figure 6. In particular, the data collected during the period 1994-1995 show irregular behavior. Additionally, the original data lack some information (less than one percent) for all series over the period 1999-2000. An imputation technique based on moving averages, which takes into account past and future values of the series, was used to replace missing values. The series were grouped by geographical region and then plotted (Figure 7). Similar patterns are observed for the series across each group.

An exploratory data analysis was carried out for each of the 15 series. There exists significant autocorrelation of order of at least one in all series. Seasonal components are present in most of partial autocorrelations. Because there is no large departure from the weakly stationary assumptions (i.e., constant means and variances), all series can be modeled using the Box-Jenkins approach. Specifically, seasonal ARIMA models can be fitted to each single series with a small number of parameters (i.e., $p\le 5,$ $q\le 5$ and $d\le 2$ ).

Figure 7.

The fifteen NDVI series grouped by area.

3.6. Clustering

Using the NVDI data set described in Section 3.5, the distance measure $D$ in (13) was computed for all possible pairs. Then, dendrogram plots were constructed using a hierarchical procedure (i.e., complete linkage) to compare the mergers and clusters produced using $D$ with those produced using the conventional distances described in Section 2.2 and (13). In Figure 8, we observe the agglomeration produced by using the Euclidean distance (top left); the other five plots show the results produced by distance (13) for different values of $k$ and $h.$ The agglomeration algorithm using the Euclidean distance merges series 5 together with series 12-15 and thus does not preserve location when grouping series. However, with $k=1$ and $h=1$ in (13), the location dependence of the 15 series is captured. Higher values of $k$ and $h$ do not modify the original clusters formed using $D$ . In Figure 9, we see the clusters and merges yielded by using the Minkowski distance in (13). Note that series 4 is classified together with series 8, 9, 10 and 11 in the top left plot, but with $k=1$ and $h=1,$ the algorithm handles the series differently (shown in the top right plot) by merging together those series that are in the same location. The effects of higher values of $k$ and $h.$ are again negligible. The Fréchet distance produced unsatisfactory results. In this case, the hierarchical algorithm does not take into account the geographical location when using both the conventional distances measures and (13). For example, series 1, 2 and 3 were merged in different stages. Nevertheless, when $k$ and $h$ are increased, the algorithm using (13) still clusters the series by geographical location. Indeed, if our goal is to produce four clusters as before, the hierarchical algorithm with $h=2$ and $k=3$ produces a geographically consistent agglomeration (dendrograms not shown here). The same analysis was performed using the hierarchical algorithm with the dynamic time warping distance measure. In this case, this distance measure produced an agglomeration by geographical location and thus did not need to be modified to capture serial correlation. The result yielded with (13) produced the same outcome for all values of $k$ and $h.$

Figure 8.

Clusters produced by using a hierarchical method with the Minkowski distance and (13).

Figure 9.

Clusters produced by using a hierarchical method with the Fréchet distance and (13)

4. Concluding Remarks and Future Work

This chapter described two problems. The first problem involved image segmentation, while the second problem involved clustering time series. For the first problem, a new algorithm was proposed that enhances the segmentation yielded by a previous algorithm (Ojeda et al., 2010). Identifying the best prediction window improves segmentation based on the estimation of AR-2D processes and generalizes the previous algorithm to different prediction windows associated with unilateral processes on the plane. An analysis of the association between the original and fitted images relies on the selection of a suitable image quality measure. Using three image quality coefficients that are commonly used in image segmentation, we carried out experiments that support our algorithm. Specifically, a set of images belonging to the image database (http://sipi.usc.edu/database/) were processed and provided satisfactory results (not shown here) in terms of image segmentation.

This chapter also proposed an extension of the dissimilarity measure first introduced in (Chouakria & Nagabhushan,2007). The simulation experiments performed and the data analysis carried out for relevant ecological series show that the distance lag $h$ plays an important role in capturing the higher-order correlation of each series. Cluster analysis performed using the proposed distance measure produced different merges and dendrograms. Furthermore, the percentage of times that the hierarchical algorithms correctly classified the highly correlated sequences increased in all cases in which the distance measure (13) was used. For the NDVI series discussed in Section 3.5, the distance measure $D$ improved the performance of the Euclidean, Fréchet and Minkowski distances in the presence of high-order autocorrelation in the series. The dynamic time warping distance measure showed the best performance in capturing the serial correlation between the NDVI series, and thus, it was not necessary to introduce modified distance measures such as (13) to ensure agglomeration by geographical location.

Now, further research for the topics presented in this chapter is outlined.

Following the notation used in the Algorithm 3, consider the following residual image.

$$R_{W_i}=Z-{\widehat{Z}}_{W_i}.$$

One interesting open problem involves the characterization of the types of images and distributions associated with the segmentation produced by Algorithms 2 and 3. In addition, the definition and study of linear combinations of residual images produced by distinct prediction windows is also of interests. For example,

$I={\displaystyle \sum _{j=1}^4{a}_j}{R}_{W_i},$

where $a_j$ is a weight associated with the residual image $R_{W_i}.$

Regarding the clustering technique problem, the distribution of $D$ can be studied from a parametric point of view. This is an open problem that we expect to address in future research.

Acknowledgments

The first author was partially supported by Fondecyt grant 1120048, UTFSM under grant 12.12.05, and Proyecto Basal CMM, Universidad de Chile. The second author was supported in part by CIEM-FAMAF, UNC, Argentina.