Open Access is an initiative that aims to make scientific research freely available to all. To date our community has made over 100 million downloads. It’s based on principles of collaboration, unobstructed discovery, and, most importantly, scientific progression. As PhD students, we found it difficult to access the research we needed, so we decided to create a new Open Access publisher that levels the playing field for scientists across the world. How? By making research easy to access, and puts the academic needs of the researchers before the business interests of publishers.
We are a community of more than 103,000 authors and editors from 3,291 institutions spanning 160 countries, including Nobel Prize winners and some of the world’s most-cited researchers. Publishing on IntechOpen allows authors to earn citations and find new collaborators, meaning more people see your work not only from your own field of study, but from other related fields too.
To purchase hard copies of this book, please contact the representative in India:
CBS Publishers & Distributors Pvt. Ltd.
www.cbspd.com
|
customercare@cbspd.com
Department of Mathematics, Universidad Técnica Federico Santa María,, Chile
Silvia Ojeda
FAMAF, Universidad Nacional de Córdoba,, Argentina
*Address all correspondence to: ronny.vallejos@usm.cl
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 hthat 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ℤd, d≥2, where ℤdis endowed with the usual partial order. That is, fors=(s1,s2,…,sd), u=(u1,u2,…,ud)inℤd,s≤u if fori=1,2,.…,d,si≤ui For a,b∈ℤd,such that a≤band a≠b, we define S[a,b]={x∈ℤd|a≤x≤b}and S〈a,b]=S[a,b]\{a}.
A random field (Xs)s∈ℤd is said to be a spatial ARMA(p,q)with parameters p,q∈ℤdif it is weakly stationary and satisfies the equation
Xs−∑j∈S〈0,p]ϕjXs−j=εt+∑k∈S〈0,q]θjεs−j,E1
where (ϕj)j∈S〈0,p]and(εj)k∈S〈0,q], respectively, denote the autoregressive and moving average parameters with ϕ0=θ0=1, and (εs)s∈ℤddenotes a sequence of independent and identically distributed centered random variables with varianceσ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(q) random field. Similarly, ifq=0, the process is called a spatial autoregressive AR(p)random field. The ARMA random field is labeled as causal if it has the following unilateral representation.
Xs=∑j∈S[0,∞]ψjεs−j
with∑j|ϕj|<∞. Similar to the time series case, there are conditions on the (AR or MA) polynomials that ensure stationarity and invertibility, respectively. Let Φ(z)=1−∑j∈S〈0,p]ϕjzjand Θ(z)=1−∑j∈S〈0,q]θjzj,where z=(z1,z2,…,zd)and zj=z1j1z2j2…zdjd. A sufficient condition for the random field to be causal is that the AR polynomial Φ(z) has no zeros in the closure of the open discDdinℂd. For example, if d=2, the process is causal if Φ(z1,z2)is not zero for any z1and z2that simultaneously satisfy |z1|<1 and |z2|<1 (Jain et al., 1999).
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×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.
LetZ=Zm,n,0≤m≤M−1 ,0≤n≤N−1 , be the original image, and letX=Xm,n,0≤m≤M−1,0≤n≤N−1, where for all0≤m≤M−1,0≤n≤N−1 ,Xm,n=Zm,n−Z¯, and Z¯is the mean ofZ. Let 4≤k≤min(M,N)and consider the rearrange images
Z=Zm,n,X=Xm,n,
where0≤m≤M′−1,0≤n≤N′−1 , M′=[M−1k−1](k−1)+1,N′=[N−1k−1](k−1)+1. For all ib=1,⋯,[M−1k−1]and for all jb=1,⋯,[N−1k−1]the (k−1)×(k−1)block (ib,jb)of the image Xis defined as
BX(ib,jb)=Xr,s,
where (k−1)(ib−1)+1≤r≤(k−1)iband(k−1)(jb−1)+1≤s≤(k−1)jb. Then, the approximated image X^of Xis provided by Algorithm 1.
Algorithm 1.
For each block BX(ib,jb)
1. Compute estimatorsϕ^1(ib,jb), ϕ^2(ib,jb)of ϕ1 and ϕ2corresponding to the block BX(ib,jb)extended to:
BX′(ib,jb)=Xr,s,
where(k−1)(ib−1)≤r≤(k−1)ib,(k−1)(jb−1)≤s≤(k−1)jb.
2. Let X^ be defined on the block BX(ib,jb)by
X^r,s=ϕ^1(ib,jb)Xr−1,s+ϕ^2(ib,jb)Xr,s−1
where (k−1)(ib−1)+1≤r≤(k−1)iband (k−1)(jb−1)+1≤s≤(k−1)jb.
Then the approximated image Z^of the original image Zis:
Z^m,n=X^m,n+Z¯,0≤m≤M′−1,0≤n≤N′−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 Z^of Z.
2. Compute the residual autoregressive image given by Z−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 size512×512 (aerial), and Figure 1(b) shows the image generated by Algorithm 1 when a moving window of size 512×512is 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×128.Another segmentation with a moving of size 512×512is 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×128and 512×512respectively.
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)∈ℤ2, a strongly causal region at (m,n) is defined as
S(m,n)={(k,l)∈ℤ2:k≤m,l≤n)}−{(m,n)}E2
For a given M∈ℕ, a strongly causal prediction window is
W={(k,l)∈s(m,n):m−M≤k≤m,n−M≤l≤n}.E3
In particular, ifM=1, then a strongly causal prediction window containing three elements is
W1={(k,l)∈s(m,n):m−1≤k≤m,n−1≤l≤n}E4
The set W1 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 ℤ2.The definition of such sets generates the prediction windows W2,W3,andW4, as shown in Figure3(b). Algorithms 1 and 2 were implemented using the prediction windows W1,W2,W3, and W4,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 W2. 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, (Xs)s∈Dand(Ys)s∈D, D⊂ℤd.For a given h∈D, the codispersion coefficient is defined as
ρ(h)=γ(h)VX(h)VY(h),E5
where s,s+h∈D,γ(h)=E[X(s+h)−X(s)][Y(s+h)−Y(s)],VX(h)=E[X(s+h)−X(s)]2(and similarly forVY(h)).
For d=2,the sample codispersion coefficient is defined by
ρ^(h)=∑s,s+h∈D'asbsV^X(h)V^Y(h)E6
withs=(s1,s2),h=(h1,h2),D′⊆D,#(D′)<∞, as=X(s1+h1, s2+h2)−X(s1, s2),bs=Y(s1+h1, s2+h2)−Y(s1, s2),V^X(h)=∑s,s+h∈D′as2, and V^Y(h)=∑s,s+h∈D′bs2.
Figure 4.
a)-(d) Images generated by Algorithm 2 with prediction windows W1−W4respectively.
where X¯ is the mean of(Xs)s∈D, SXis the standard deviation of (Xs)s∈D, and SXY is the covariance between (Xs)s∈Dand (Ys)s∈D(and similarly for Y¯andSY). The quantity C=SXY/SXSYmodels the linear correlation between (Xs)s∈D and(Ys)s∈D, M=2X¯Y¯/(X¯2+Y¯2)measures the similarity between the sample means (luminance) of (Xs)s∈Dand(Ys)s∈D, and V=2SXSY/(SX2+SY2)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 (Xs)s∈Dand(Ys)s∈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 CQindex, which is defined as:
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 W2.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 W1,W2,W3,and W4. 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 W4.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 Zand the fitted image generated by Algorithm 1 with the prediction window Wi by Z^Wi,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 ZandZ^Wi. This criterion is summarized in the following algorithm.
Algorithm 3.
1. Use Algorithm 1 to generate the approximated images Z^Wiof Z,i=1,2,3,4.
2. Compute an image quality index between Zand Z^Wifor all i=1,2,3,4.Suppose that the maximum value for the image quality index is attained for Z^Wj,1≤j≤4.Then, the best fitted image is Z^Wj.
3. Compute the residual autoregressive image Z−Z^Wj.
3.1. Measuring Closeness and Association Between Time Series
Let x=(x1,x2,…,xp)and y=(y1,y2,…,yq) be two time series. There are several conventional distance measures between time series. For example, ifp=q=n, then the Euclidean distance between xand yis defined as dE(x,y)=(∑i=1n(xi−yi)2)1/2.As is evident, dEignores information about the dependence between xand y. The Minkowski distance is a generalization of the Euclidean distance, which is defined as
dM(x,y)=(∑i=1n(xi−yi)q)1/q,E9
where qis 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≤min(p,q). Let Mbe the set of all mappings rbetween xand ysuch that ris a sequence of mpairs preserving the order
r=((xa1,yb1),(xa2,yb2),…,(xam,ybm)),
where ai∈{1,2,...,p},bj∈{1,2,...q} with a1=1,b1=1,am=p,bm=q and fori∈{1,2,…,m−1},ai+1=(aiorai+1)and bi+1=(biorbi+1).Note that |r|=maxi=1,2,…,m|xai−ybi| is the mapping length representing the maximum span between two coupled observations. Thus, the Fréchet distance between the series xandyis given by
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,
|r|=∑i=1,2,…,m|xai−ybi|.
Dynamic time warping is then defined as
dDTW(x,y)=minr∈M|r|=minr∈M∑i=1,2,…,m|xai−ybi|.E11
The distances defined above are based on the proximity of the values |xai−ybi|. However, these distances disregard both the temporal dependence between the sequences xand 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
where β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={Xs:s∈D⊂ℤ}and Y={Ys:s∈D⊂ℤ},and let xand ybe realizations of these processes as in Section 3.1. For d=1,the estimator (6) becomes
where N(h)={t:t+h∈D},N=|N(h)|is the cardinality of N(h),and sequences xand yare realizations of processes xandy, respectively. The coefficient ρ^(h) is called the comovement coefficient when h=1. Although ρ^(h)is not the correlation coefficient, the codispersion coefficient shares a number of its standard properties. For example, ρ^(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 +1indicates 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 ρ^(h)=0expresses that there is no monotonicity between xand yand 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 xand yand d=1is given by
D(x,y,h)=f(ρ^(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 xandy. There are many possible ways to choose a function f. Here, we follow the guidelines given in (Chouakria & Nagabhushan, 2007), according to which fis considered an exponential adaptive tuning function given by
fk(t)=21+exp(kt),E14
where k≥0modulates the contributions of the proximity with respect to values and behavior. For example, when |ρ^(1)|is large and k=2,the proximity with respect to behavior contributes 76.2% toD. The flexibility of Dallows us to choose ksuch 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 his arbitrary, while Chouakria's index only captures the first-order correlation.
The dependence of (13) on his crucial, and in some specific cases, hcan be chosen using an optimal criterion. For example, for two AR(1) processes with parameters ϕ1andϕ2, respectively, and a correlation structure between the errors (Rukhin & Vallejos,2008
), it is possible to find ϕ1,ϕ2such that Var(ρ^(1))Var(ρ^(2)).In other words, for those processes in which the asymptotic variance of the codispersion coefficient is known, we suggest setting the value of h^to produce the minimum variance. That is,
h^=argminh[Var(ρ^(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.
Xti=ϕiXt−1i+εti,i=1,2,...6,
where ∀i=1,2,...,6,Xi={Xti}t∈ℤdefine the i−model, and the sequence εi={εti}t∈ℤis zero-mean white noise. Note that the sequences ε1and ε2have the covariance structure
Cov(εt1,εs2)={ρστ,s=t,0,otherwise,
with σ2=τ2=1,andρ=0.9. The same covariance structure is assumed for εt1and εt3,with σ2=τ2=1and ρ=0.7.εi,i=4,5,6,are assumed to be white noise processes with variance 1 and are uncorrelated with all other error sequences. Ifi,j≤3, the correlation structure between processes XiandXj, i≠j, is not null due to the correlation structure between εiandεj. Instead, ifi,j≥4, i≠j, the correlation structure between processes Xiand Xjvanishes.
Two hundred observations were generated from each model for ϕ1=−0.5,ϕ2=−0.3,ϕ3=−0.7,ϕ4=0.1,ϕ5=−0.9,and ϕ6=−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,andk=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 usingD.
In Figure 5, we see that the dendrogram obtained using hierarchical methods with the Euclidean distance does not recognize the correlation structure between X1and X3.In this case, sequences X1,X2,X4,and X5are pulled together before sequence X3.However, hierarchical methods using (13) yield the expected results, combining sequences X1,X2,andX3before 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,2and 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). DDTW(h,k)is distance measure (13) with the DTW distance. DM(h,k)and DF(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≤5,q≤5andd≤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 Din (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 Dwith 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 kand 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=1and h=1in (13), the location dependence of the 15 series is captured. Higher values of kand h do not modify the original clusters formed usingD. 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=1and 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 kand 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 kand hare 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=2and k=3produces 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 kand 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)
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 hplays 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 Dimproved 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.
RWi=Z−Z^Wi.
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=∑j=14ajRWi,
where ajis a weight associated with the residual image RWi.
Regarding the clustering technique problem, the distribution of Dcan be studied from a parametric point of view. This is an open problem that we expect to address in future research.
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.
References
1.AllendeH.GalbiatiJ.VallejosR.2001Robust image modeling on image processing.Pattern Recognition Letters22111219 1231
2.BaranS.PapG. ZuijlenM. C. A.2004Asymptotic inference for a nearly unstable sequence of stationary spatial AR models.Statistics & Probability Letters69153 61
3.BasuS.ReinselG.1993Properties of the spatial unilateral first-order ARMA model.Advances in Applied Probbability253631 648
4.BustosO.OjedaS.VallejosR.2009aSpatial ARMA models and its applications to image filtering.Brazilian Journal of Probability and Statistics232141 165
5.BustosO.OjedaS.RuizM.VallejosR.FreryA.2009bAsymptotic Behavior of RA-estimates in Autoregressive 2D Gaussian Processes.Journal of Statistical Planning and Inference1391036493664
6.CariouC.ChehdiK.2008Unsupervised texture segmentation/classification usind 2-D autoregressive modeling and the stochastic expectation-maximization algorithm.Pattern Recognition Letters297905917
7.ChoiB.2000On the asymptotic distribution of mean, autocovariance, autocorrelation, crosscovariance and impulse response estimators of a stationary multidimensional random field.Communications in Statistics-Theory and Methods29817031724
8.ChouakriaA. D.NagabhushanP. N.2007Adaptive dissimilarity for measuring time series proximity.Advances in Data Analysis and Classification115 21
9.CliffordP.RichardsonS.HémonD.1989Assessing the significance of the correlation between two spatial processes.Biometrics451123 134
10.CullisB. R.GlessonA. C.1991Spatial analysis of field experiments-an extension to two dimensions.Biometrics47414491460
11.FrancosJ.FriedlanderB.1998 Parameter estimation of two-dimensional moving average random fields.IEEE Transaction Signal Processing46821572165
12.GaetanC.GuyonX.2010Spatial Statistics and ModellingSpringerNew York
14.GolayX.KolliasS.StollG.MeierD.ValavanisA.1998A new correlation-based fuzzy logic clustering algorithm of FMRI.Magnetic Resonance in Medicine402249260
15.GrondonaM. R.CrossaJ.FoxP. N.PfeifferW. H.1996Analysis of variety yield trials using two-dimensional separable ARIMA processes.Biometrics522763770
16.GuoJ.BillardL.1998Some inference results for causal autoregressive processes on a plane.Journal of Time Series Analysis196681691
17.GuyonX.1982Parameter estimation for a stationary process on a d-dimensional lattice.Biometrika69195105
18.HoP. G. P.2011Image segmentation by autoregressive time series model. In Image Segmentation edited by Pei-Gee Ho. InTech.
19.IlligA.Truong VanB.2006Asymptotic results for spatial ARMA Models.Communications in Statistics-Theory and Methods354671 688
20.JainA. K.MurtyM. N.FlynnP. J.1999Data Clustering: A Review.ACM Comput. Surveys,313264323
21.KashyapR.EomK.1988Robust images techniques with an image restoration application. IEEE Trans. Acoust. Speech Signal Process3681313 1325
22.MacQueen. J. B.1967Some Methods for classification and Analysis of Multivariate Observations.Proceedings of 5 -th Berkeley Symposium on Mathematical Statistics and ProbabilityBerkeleyUniversity of California Press1281297
23.MartinR. J.1996Some results on unilateral ARMA laticce processes.Journal of Statistical Planning and Inference503395411
24.MatheronG.1965Les Variables Régionalisées et leur EstimationMassonParis.
25.PolitisD. N.RomanoJ. P.1994The stationary bootstrap.Journal of the American Statistical Association894281303 1313
26.OjedaS. M.VallejosR. O.LuciniM.2002Performance of RA Estimator for Bidimensional Autoregressive Models. Journal of Statistical Simulation and Computation72147 62
27.OjedaS. M.VallejosR.BustosO.2010A New Image Segmentation Algorithm with Applications to Image Inpainting.Computational Statistics & Data Analysis5492082 2093
28.OjedaS. M.VallejosR.LambertiW. P.2012 A Measure of Similarity Between Images, in press Journal of Electronic Imaging.
29.QuintanaC.OjedaS.TiraoG.ValenteM.2011 Mammography image detection processing for automatic micro-calcification recognitionChilean Journal of Statistics226979
30.RukhinA.VallejosR.2008Codispersion coefficient for spatial and temporal series.Statistics and Probability Letters78111290 1300
31.VallejosR.MardesicT.2004A recursive algorithm to restore images based on robust estimation of NSHP autoregressive models.Journal of Computational and Graphical Statistics133674 682
32.VallejosR.Garcia-DonatoG.2006Bayesian analysis of contaminated quarter plane moving average models.Journal of Statistical Computation and Simulation762131147
33.VallejosR.2008Assessing the association between two spatial or temporal sequences.Journal of Applied Statistics35121323 1343
34.TjostheimD.1978A measure of association for spatial variables.Biometrika651109 114
35.WangZ.BovikA.2002A universal image quality index.IEEE Signal Processing Letters9381 84
Written By
Ronny Vallejos and Silvia Ojeda
Submitted: 20 February 2012Published: 24 October 2012