Robust Estimation for Two-Dimensional Autoregressive Processes Based on Bounded Innovation Propagation Representations

Robust methods have been a successful approach to deal with contaminations and noises in image processing. In this paper, we introduce a new robust method for two-dimensional autoregressive models. Our method, called BMM-2D, relies on representing a two-dimensional autoregressive process with an auxiliary model to attenuate the effect of contamination (outliers). We compare the performance of our method with existing robust estimators and the least squares estimator via a comprehensive Monte Carlo simulation study which considers different levels of replacement contamination and window sizes. The results show that the new estimator is superior to the other estimators, both in accuracy and precision. An application to image filtering highlights the findings and illustrates how the estimator works in practical applications.



There are no comments yet.


page 3

page 15

page 16


Conditional Maximum Lq-Likelihood Estimation for Regression Model with Autoregressive Error Terms

In this article, we consider the parameter estimation of regression mode...

More on the restricted almost unbiased Liu-estimator in Logistic regression

To address the problem of multicollinearity in the logistic regression m...

Consistency Results for Stationary Autoregressive Processes with Constrained Coefficients

We consider stationary autoregressive processes with coefficients restri...

Robust Estimation in High Dimensional Generalized Linear Models

Generalized Linear Models are routinely used in data analysis. The class...

Robust graphical lasso based on multivariate Winsorization

We propose the use of a robust covariance estimator based on multivariat...

Objective Bayesian Analysis of a Cokriging Model for Hierarchical Multifidelity Codes

Autoregressive cokriging models have been widely used to emulate multipl...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Robust statistics methods arise in a wide range of applications; in particular, the words “robust” and “robustness” frequently appear in the context of image and signal processing. Filters based on robust statistics such as the median filter, provide basic mechanisms in image processing (Ays ; Hua ; Tzu ). In general, image processing based on robust strategies has shown remarkable ability in image restoration and segmentation (see, for example, Ben ; Ji ; Tar ). High-level image analysis, or vision, also benefits from the use of robust estimation techniques, as it can be seen in Com ; Got ; Kim ; Pra and Sin .

There is a large number of examples such as signal analysis and processing on the plane, imaging remote sensing and design of experiments in agronomy, in which data is recorded in a grid or lattice in . A class of 2D autoregressive processes has been proposed (Whi ) as a set of plausible models for the spatial correlation in such data (Tjo ). These models are natural extensions of the autoregressive processes used in time series analysis (Bas ).

Thus, most of the robust procedures suggested for time series parametric models have been generalized for spatial parametric models when the process has been contaminated with innovation or additive outliers (

Kas ). Because a single extreme value is capable of introducing significant distortions in estimators, most initiatives are focused on providing estimators that are robust to the appearance of anomalous data.

In this context, there are at least three classes of robust estimators that have been studied. They are the estimators M, GM and RA. The M estimators were successfully used to develop and implement a computational algorithm for image restoration, based on 2D autoregressive models (Kas ). Later, All1 analyzed the development and implementation of Generalized M (GM) estimators for the same class of models. The robust estimators of Residual Autocovariance (RA) were presented by Bus3 in the context of the time series, where in the procedure of recursive estimation the residues are cleaned through the application of a robust function. The extension of the RA estimators for two-dimensional autoregressive models and their corresponding computational treatment were developed by Oje4 . Monte Carlo simulation studies show that the performance of the RA estimator is higher than the M estimator and slightly better than the GM estimator when the model has been contaminated with additive type outliers. In addition, Bus5 studied the asymptotic behavior of the RA estimator for unilateral two-dimensional autoregressive processes, generalizing the results for the asymptotic behavior of one-dimensional series previously established by Bus4 . We do not yet know the asymptotic properties associated with the M and GM estimators when the model is contaminated.

One of the reasons why the spatial autoregressive model (AR-2D) has been extensively used in image analysis and processing is its remarkable ability to represent a variety of real scenarios without the need to use a large number of parameters. However, the robust estimators developed so far for the parameters of the AR-2D model have been constructed only under the assumption of innovative or additive type random noise; there are no proposals for parameter estimation when the model is contaminated from another more general pattern of noise. In this work we define a new class of robust estimators for the AR-2D contaminated models. This class of estimators is robust under replacement contamination that includes additive type contamination. To avoid the propagation of the effect of an outlier when calculating the innovative residuals in the AR-2D model we propose a new approach that consists of defining this residuals using an auxiliary model. For this reason, bounded innovation propagation AR-2D (BIP-AR-2D) models are proposed. With the help of these models we suggest an estimator for the AR-2D model that is robust when the process contains outliers.

Our proposal is the generalization of the two- dimensional case of the BMM estimators, developed by Mull for ARMA time series models. The approach of generalizing one-dimensional proposals to the case of two or more dimensions is the strategy that is naturally used in many areas of science. The work is relevant because it offers a new way of robustly estimating the parameters in the two-dimensional autoregressive model. It is a contribution that can be beneficial in the field of signal processing, spatial statistics and, in general, in areas of applied mathematics that use this model intensively, for example, to represent texture images such as robotic vision, recognition of patterns and matching problems. Although estimators of the model parameters have been developed so far, the central objective of this work is to provide a new estimator that betters the estimation proposals already known.

The rest of the paper is organized as follows. In Section 2, the basic definitions are presented. We first presented some background material on bidimensional autoregressive processes (AR-2D) and parameter estimators of the model. We also define procedures for generating replacement contamination in such models. In Section 3, the new model BIP-AR-2D for spatial processes and the new estimator of the AR-2D model parameters are presented. In Section 4, several Monte Carlo studies are carried out to evaluate the performance of the new estimator against different contamination schemes, compared to the LS, M, GM and RA estimators. Section 5 presents two applications to real images that demonstrate the capabilities of the BMM-estimator to represent, segment and restore contaminated images. Conclusions and future works appear in Section 6. The results of the Monte Carlo studies (Section 4) are shown in the Appendix.

2 Preliminaries

2.1 The spatial ARMA models

A great variety of texture can be generated through the two-dimensional autoregressive moving average models (ARMA-2D models). For example, Figure 1 shows textures generated for the particular case of ARMA-2D processes, the AR-2D process, with two and three parameters (Fig. 1 (a)-(b) and (c)-(d) respectively). Besides, in the last years, several theoretical properties have been studied, and multiple applications have been developed for these processes (Bar ; Vall3 ; Bus6 ; Qui ; Zie ; Yao ; Sad ).

(a) (b) (c) (d)

Figure 1: Autoregressive Processes. (a) , , ; (b) , , ; (c) , , ; (d) , , .

Spatial autoregressive moving average (ARMA) processes can be defined over random fields indexed over where considers the usual partial order, i.e., for in if for . We define and for such that and .

A random field is said to be a spatial ARMA of order if it is weakly stationary and satisfies the equation


where and denote the autoregressive and moving average parameters respectively with and

denotes a sequence of i.i.d. random variables with variance

Note that if the process is called spatial autoregressive AR() random field. An ARMA random field is called causal if it can be represent for the following equation:


with Similar to the time series case, there are conditions on the (AR or MA) polynomials for them to be stationary and invertible respectively (Bas ). As an example, consider a first-order autoregressive process as in (1) with , and . Then and the model is of the form


where to simplify the notation it took , and . In equivalent form, (3) can be expressed as


where and are the backward operators given by , and in (4), . In the case that has inverse, we can write equation (2) as

Bas studied the correlation structure of a process like (3). They obtained conditions to guarantee the existence of the stationary representation of the model (3) as in (2). In that case, the use of a multinomial expansion for implies the convergent representation


where with are the coefficients of this multinomial expansion. The increase in the number of parameters in the model also increases the diversity of possible textures but in contrast, the calculations become more complex. In this paper we worked with the AR 2D model with three parameters as in (3).

2.2 Types of contamination in AR-2D processes


described (Chapter 8) some probability models for time series outliers, including additive outliers (AOs), replacement outliers (ROs) and innovation outliers (IOs).

In this section we generalized the notion of replacement outliers for spatial processes. Let be a stationary process, for example an AR-2D process as (3) and let be the observed process. It is said that process behaves like a two-dimensional Replacement Outlier model (RO) if it is given by


where is a zero-one process such that and and is a replacement process that is not necessarily independent of . The fraction is positive and small.
A particular case of the RO models is the two-dimensional Additive Outlier model (AO), in which

is a stationary process independent of and is a Bernoulli process. This type of contamination is very important for satellite image processing; for example, it is present in optical images such as those from Landsat satellites. When does not follow the pattern of AO, we say that the contamination process follows a Strictly Replacement Outlier model (SRO).

2.3 Robust Parametric Estimation

Because LS estimators are very sensitive to the presence of atypical values (Martin, 1980), several alternative estimators arise to mitigate the impact of contaminated observations on estimates. Most of these proposals are natural extensions of robust estimators studied in time series.

Robust estimators have been defined for models with a small amount of parameters. Here, we summarize the well-known robust estimators for the model (3); however, a more general development can be found for the AR and MA models in Kas , All4 , Oje4 , Vall3 and Bus5 .

Note that model (3) can be rewritten in the linear model form:


is a parameter vector and

To obtain the LS estimator of , we minimize the following function:



Similarly, the class of M estimators for 2D autoregressive processes (Kas ) is defined by the minimization of the function


Because the M estimators are very sensitive when the process is contaminated with additive outliers is that other robust estimators arise to reduce the effects of additive outliers. Alternatively, All4 developed the generalized M estimators (GM) for spatial AR processes. A GM estimator of is the solution to the problem of minimizing the non-quadratic function defined by:

where is as a function like the one that was defined in (8), and are weights corresponding to the respective

In Oje4 , the authors presented the robust autocovariance (RA) estimator for AR-2D processes which was first defined for time series models by Bus3 . This estimator is determined by the following equations

where is a function that depends of the residues, are coefficients that depend on the parameters and is a independent estimate

3 A new approach

3.1 BIP-AR 2D models

A new class of bounded nonlinear AR-2D models is presented in this work: the bounded innovation propagation AR-2D model (BIP-AR 2D). This model arises from the need to estimate the best possible parameters of an autoregressive central model when a contaminated process is observed. The BIP-AR 2D model is a two-dimension generalization of the model presented for time series by Mull .

Given a stationary and invertible AR-2D model like in (3), it supports a stationary representation as in (5). We define the BIP-AR 2D auxiliary model as:


with , ’s are i.i.d. random variables with symmetric distribution and

is an odd and bounded function. Besides,

is a robust M-scale of

, which coincides with the standard deviation if

are normal, and is defined as the solution of the equation .

Note that (9) can also be written as

and multiplying both members by , we get

which is equivalent to

3.2 BMM estimator for AR-2D processes

In time series, Mull introduced the MM- estimators for ARMA models based in the definition of MM-estimate for regression where the residuals are calculated as in the BIP-ARMA model instead of just as in the pure ARMA model. The idea of MM-estimators in regression is to compute a highly robust estimator of the error scale in a first stage, and this estimated scale is used to calculate an M-estimator of the regression parameters in a second stage. However, in time series this differs somewhat because an MM-estimate is not enough to guarantee robustness.
In the same way that residues of AR-2D model exist (7), there are residues obtained from BIP-AR 2D model:


for all . With this residue the objective function that must be minimized to obtain the M-estimator of the parameters under a model BIP-AR 2D is defined:


where is a robust estimate of .
One way of robustly estimating the scale was introduced in 1964 by Hub as follows: Given a sample , with , an M-estimate of scale is defined by any value satisfying


where is a continuous and non- constant function, non-decreasing in and symmetric around zero as well. To make the M-scale estimate consistent with the standard deviation when the data are normal, it requires that

under the standard normal distribution. Taking

, we get a maximum breakdown point of 0.5. With all this, we can define the new BMM-2D estimator by following the two steps given below:

First Step: At this stage, an estimate of is obtained. For this purpose, two estimates are considered: one using an AR-2D model, another using a BIP-AR 2D model; then we choose the smallest of them.
Let a continue, non-constant, non-decreasing in , bounded and symmetric function and such that: . This guarantees that for a normal random sample, the M-scale estimator based on converges to the standard deviation and the breakdown point of is 0.5. Put

for some . Then, we defined an estimate of :

and the corresponding estimate of :


where ,
with as in (7) and is the M-estimate of scale based on and defined as in (12). Later, we described the estimate corresponding to the BIP-AR model. Define by the minimization of over . The value is an estimate of computed as if were the true parameters and the ’s were normal. Then, from (9), because in this case the M-scale coincides with the standard deviation of , we had:

where and . Let a robust estimate of and where . Then, we defined

The scale estimate corresponding to the BIP-AR-2D model is defined by


where, for simplify, we denoted and

with defined as in (3.2). Finally, our estimate was

Second Step: We considered a bounded function that satisfies the same properties as but also . This function was chosen such that the corresponding M-estimator is highly efficient under normal innovations. Given the objective functions defined in (8) and (11) with the scale obtained in the first step ():


The corresponding M-estimators of the parameters for each function are:

Then, we defined the BMM-estimator 2D as:

4 Monte Carlo Results

The aim of this section is to analyze the performance of the new BMM estimator to estimate the parameters in the model (3) compared to the LS, M, GM and RA estimators. We performed several experiments. Each experiment is based on different Monte Carlo studies. We set the parameter values of (3) as:


It is important to mention that the set parameters were chosen at random satisfying the conditions of invertibility of Bas . We performed our study under five different conditions of the model (Cases); in Case I, the model was non-contaminated, while in Cases II, III, IV and V, the model was contaminated according to (6):

  • Case I) Non-contaminated model like in (14), where is a normal distribution process with and for all .

  • Case II) AO, where the process is independent of the process and follows a normal distribution with zero mean and variance 50.

  • Case III) SRO, where the replacement process follows a t-student distribution with 2.3 f.d.

  • Case IV) SRO, where the replacement process is another autoregressive process, independent of the process, with parameters , and .

  • Case V) SRO, where the replacement process

    is a white noise process with normal distribution with zero mean and variance 50.

In each of the five variants of the model (14), the parameters and were estimated by the five procedures presented in the previous sections. In each experiment, 500 simulations of the model were generated, and the mean value, the mean square error (MSE) and the sample variance were computed. For the contaminated models we considered four levels of contamination: 5%, 10%, 15% and 20%. Besides, we performed our study considering different window sizes: , , , and . For the calculation of the BMM estimator, a robust estimator of the scale was obtained as in (13). was selected according to the same criteria that was taken for the definition of the BMM estimator for time series present in Mull , where the function is given by:

The same function was used to calculate the M estimators. In addition, for the implementation of the GM estimator the weights were set:

where is the following version of the Huber function:

Finally the RA estimators were implemented according to the details formulated in Oje . To facilitate the paper reading, only the boxplots of the simulations have been included in the body of the work; the numerical Monte Carlo results are shown in the Appendix.

4.1 Experiments

In a first experiment, we studied the performance of the BMM estimator for the non-contaminated model (Case I). All the methods estimated the parameters quite well. Table 4 shows the results obtained for the four different window sizes considered. In Figure 2, the corresponding boxplots are shown. In this case, it is convenient to use the LS method due to its simplicity and calculation speed.
The second experiment was developed in the context of Case II. We analyzed the ability of the BMM method to estimate the parameters of the model, considering a 10% of additive contamination, and for window sizes , , and , in comparison with the LS, M, GM and RA methods. Table 4 shows the estimated values for , and , using the different five procedures analyzed. Figure 3 exhibits the corresponding boxplots. For window size and , it can be seen that the BMM estimator is the best because its values are closer to the parameter than the estimates produced by the other methods mentioned. In addition, BMM estimator has the lowest variance and the lowest MSE. When the window size was or , the best performance corresponded to the GM and RA estimators; however, the BMM estimation values were similar to the RA and GM estimations. An analogue affirmation is valid to the sample variance and MSE of BMM. We also noted that for any window size, the M estimator had a very small sample variance but their estimations were wrong when compared to the ones in the other methods.
Figure 2: Case (I): LS, M, GM, RA and BMM estimation boxplots for (a) , (b) and (c) ; model (14) without contamination, varying the window sizes.


Figure 3: Case (II): LS, M, GM, RA and BMM estimation boxplots for (a) , (b) and (c) , varying the window sizes. Model (14) with additive contamination level, with a normal noise.
Figure 4: Case (II): LS, M, GM, RA and BMM estimation boxplots, for (a) , (b) and (c) in model 14, with additive contamination, varying the contamination level, with a window size. The contamination process is a normal noise, with .


Figure 5: Case (III): LS, M, GM, RA and BMM estimation boxplots, for (a) , (b) and (c) in model 14, with contamination of replacement, varying the contamination level, with a window size. The process of contamination follows a t-Student distribution with 2.3 d.f.

The third experiment also refers to Case II. We set the window size at and varied the additive contamination level, considering four levels: , , and . The BMM method was the best in most of the cases studied, followed by the RA estimator. This behavior is deduced from the comparison of the values estimated by the BMM method with the respective estimates obtained by the other procedures. The values of the dispersion measures also point out that the BMM estimator is the most accurate methodology. These results can be seen in Table 6. In addition, from Figure 4 we can note that using any of the five estimators, the parameter was estimated with less precision than and as the percentage of contamination increases. Besides, while was underestimated by all methods, for all levels of contamination, the RA estimator was the only method that overestimated and independently of the contamination level. The fourth experiment is related to Case III. The process of contamination is a replacement contamination where the replacement process follows a t-student distribution with 2.3 f.d.. The simulations were performed for a window. Table 6 and plot 5 show the results. Boxplots show that the BMM method is the best performing estimator, followed by the RA, GM, M and LS methods, in that order. We also noted that in all methods the estimates deteriorate as the level of contamination increases. Additionally, the classic LS estimator presented greater dispersion of the data. The fifth experiment was performed in the context of Case IV, where the replacement process was a autoregressive process. We set the window size at , varying the level of contamination ( and ). Table 8 shows these results. Besides, Figure 7 displays the boxplots corresponding to these tables. We can see a pattern similar to the fourth experiment; except that in this case, the variances seem very much alike. Finally, the sixth experiment was carried out according to Case V. The replacement process of the contamination was a white noise with variance 50. As in the fifth experiment, we set the window size at , varying the level of contamination. Table 8 presents the estimated values obtained. The corresponding boxplots are shown in Figure 8. The parameter values and were underestimated for all methods, excluding the RA estimator that overestimated the values of and parameters. The BMM estimator was less affected by the contamination process. The LS and M estimators are less accurate than the GM, RA and BMM estimators. Comparatively, the RA estimator presented the highest variance, whereas the GM estimator, although quite accurate, deteriorates more than the BMM estimator as the level of contamination increases.

4.2 Computational Time Evaluation

All the computational routines were developed in R and were carried out on the server JupiterAce of FaMAF-UNC. It has 12-cores 2.40GHz Intel Xeon E5-2620v3 processor, with 128 GiB 2133MHz of available DDR4 RAM. Running time as time logarithm of a single simulation to each estimator vs the window size in Case I is presented in Figure 6. Time was expressed in seconds. The graph shows that the computational cost of the RA estimator is the highest; for example, in a window size, the RA running time was 43.812 seconds, while BMM, GM, M and LS computational costs were 2.936, 0.552, 0.516 and 0.436 seconds, respectively. This results show that, although RA estimator is one of the major competitors of BMM estimator, due to its accuracy and good asymptotic properties, it exhibits its computational cost as a disadvantage. This makes RA an unattractive estimator for the processing of big size images.

Figure 6: Logarithm of the estimation time (in seconds) when the process has additive contamination of according to the window size.


Figure 7: Case (IV): LS, M, GM, RA and BMM estimation boxplots, for (a) , (b) and (c) in model 14, varying the contamination level, with a window size. The contamination process is of replacement type, by an AR process with , and parameters.


Figure 8: Case (V): LS, M, GM, RA and BMM estimation boxplots, for (a) , (b) and (c) in model 14, varying the contamination level, with a window size. The contamination process is of replacement type, with a white noise, of variance 50.

5 Application to real images

The analysis of contaminated images is of great interest in several areas of research. For example, the reconstruction of contaminated images is relevant in modeling of images (All2 ; Vall2 ), and, in general, the reduction of the noise produced by interferences taking place in the processes of obtaining the physical image and transmitting it electronically plays an important role in the literature (Bus1 ).
In Oje2 , two algorithms for image processing based on the unilateral AR-2D model with two parameters were presented. The foundations of the algorithms are random field theory and robustness for spatial autoregressive processes. The first one produces a local approximation of images, and the second one, is a segmentation algorithm. In this work, we proposed to use a variant of these algorithms using a unilateral AR-2D process with three parameters (model (3)), instead of two parameters as it was originally proposed. We called the modified algorithms as Algorithm 1 and Algorithm 2. We applied them to reconstruction and segmentation of images using the LS, GM and BMM estimators of the parameters in the model (3). Later, we inspected and compared the performance of these estimators in Algorithms 1 and 2 on contaminated images. To compare the images generated by the algorithms and, therefore, the performance of the different estimators, we calculated three indexes used in the literature; the SSIM index (Wan ), the CQ index (Oje3 ), and CQmax index (Pis ). Next, we present two numerical experiments using the image “Lenna”, which was taken from the USC-SIPI image database In Figure 9-(I), the original image is shown.
In the following we present Algorithms 1 and 2, for more details about the notation you can see the work Oje2 .

0:  Original image .
0:  Approximated image of the original image
1:  Define as
2:  Generate block
3:  Compute the estimations , , of , and corresponding to the block extended to
4:  Define on the block by
where and
5:  Define as
Algorithm 1 Local approximation of images by using AR-2D processes
0:  Original image
0:  Segmentated image
1:  Generate an approximated image of with the Algorithm 1.
2:  Compute the residual image defined as .
Algorithm 2 Segmentation

In the first experiment, Algorithm 1 was applied to image representation. We locally adjusted an AR-2D process to the original image for different window sizes, and estimated the parameters of the model with the BMM estimator. Fig. 9, (a), (b), (c) and (d) exhibits the BMM reconstructed images obtained respectively using the window sizes , , and . For all window sizes, the BMM recostructed images are visually good; although a quantitative analysis of the similarity between each BMM reconstructed image and the original image showed differences. We calculated the SSIM, CQ (1,1) and CQ index, between each reconstructed image and the original image. The three indexes revealed that the similarity decreases as the size of the window increases (Table 1); so the best fits were obtained with small window sizes. This result reflects the assumption that the two-dimensional autoregressive model is a local adjustment model. Next, we applied Algorithm 2 and generated four difference images (e), (f), (g) and (h) shown in Fig. 9. We observed that the difference image (h) highlights the edges more than the others do. This shows that when we performed the reconstruction with a window size, (Fig. 9 (d)), a lot of information got lost and this is reflected by the difference image (Fig. 9 (h)).

(a) (b) (c) (d)
(e) (f) (g) (h)
Figure 9: Image (I) Lena original image. (Right) The first row has the reconstructions done for BMM estimator adjusting window sizes , , and respectively ((a) to (d)). The second row has the respective differences ((e) to (h)) with respect to the original image (I).
Window size SSIM CQ(1,1) CQ
0.9948914 0.8582201 0.9706984
0.9827996 0.8309626 0.9544317
0.9779204 0.8151581 0.9462133
0.9762065 0.8073910 0.9423786
Table 1: SSIM, CQ and CQ index between the original image and each one of the BMM reconstructed images (a), (b), (c) and (d) in Figure 9.

In the second experiment, the original image was 10% additively contaminated (Fig. 10 (II)), and we used it as input in Algorithm 1. We obtained four reconstructed images using the LS, GM and the BMM estimators. Next, the Algorithm 2 was performed. The studies were carry out considering and window sizes. In the first two columns in Figure 10, we can observe the results obtained considering windows. Visually, there are not great differences between the different images reconstructed. When analyzing Table 2, we verify this because the measured indexes are comparable to each other. On the other hand, the third and fourth columns of Figure 10 show the results obtained by adjusting windows. It is observed that the image (l), corresponding to the difference between the image restored with BMM (l) and the one contaminated with additive noise, highlights the edges slightly more.

(I) (II)
(a) (d) (g) (j)
(b) (e) (h) (k)
(c) (f) (i) (l)
Figure 10: Image (I) Lena original image; Image (II) image with 10% additive contamination and . In the first row, adjustments made with LS; in the second row, with GM; and in the third row, with BMM. Columns 1 and 2 correspond to adjustments with windows, and columns 3 and 4 to windows. Columns 1 and 3 are the reconstruction for Algorithm 1 and columns 2 and 4 are the segmented images for Algorithm 2.
Estimate SSIM CQ(1,1) CQ
LS 0.9836079 0.8416351 0.9588826
GM 0.9390820 0.7821954 0.9103257
BMM 0.9846007 0.8328356 0.9577393
Table 2: Similarity between the original image and the reconstructions of Lena with additive contamination using window size (Figure 10-(I) vs. figures 10-(a),(b),(c)).

6 Conclusions and discussions

A new estimator called BMM was proposed to estimate the parameters in first-order two-dimensional autoregressive models with three parameters. The new estimator is a two-dimensional extension of the BMM estimator proposed by Mull , for autoregressive models of time series. We also extended the definition of replacement contamination, given for one-dimensional models (Maro ), to the case of AR-2D models; this type of contamination includes additive-type contamination. The performance of the proposed estimator for AR-2D models with replacement contamination and without contamination was analyzed. Besides, the new estimator was compared with the classical least square estimator (LS) and robust estimators M, GM and RA. The comparative analysis was performed from six experiments, each of which involved several Monte Carlo studies, considering different replacement contamination patterns and varying the level of contamination and window sizes of observation. The LS estimator produced estimates that are very sensitive to the presence of atypical values, while the other estimators had better results. Using Monte Carlo simulation, we observed that the GM, RA and BMM estimators are highly superior to the M and LS estimators. However, the new estimator showed the best behavior, in both accuracy and precision, followed by the RA estimator in accuracy and the GM in precision. An analysis of the computational cost showed that the RA estimator is the most expensive, followed by the BMM, GM, M and LS estimators, in that order. Finally, in an application to real data, we introduced a variant of the algorithm developed by Oje2 , to perform image segmentation, using an AR-2D model with three parameters, and BMM estimators. In the light of the examples shown in Section 5, we conclude that the adapted algorithm is able to highlight the borders and contours in the images.
The following proposals outline some directions for future work. In Oje , the author established the asymptotic normality and consistency of the robust RA estimator for the parameter of a two-dimensional autoregressive process. Although the estimators M, GM and BMM are reasonable to estimate parameter , their asymptotic behavior is still an open problem. In the context of image processing, in this work, the difference between a real image and an BMM approximated image was computed. The resulting image could contribute to solve problems like that border detection, classification and restoration of images. It would be interesting to explore the limitations of a segmentation method based on the difference image between a real and a fitted image. It is also important to analyze the behavior of the BMM estimator in combination with image restoration techniques. The same relevance has the study of properties of BMM estimator, in particular, and robust estimators, in general, as alternatives to the least square estimators under non causal and semi causal AR-2D models.


We thank to Ph. D. Oscar Bustos and Ph. D. Ronny Vallejos for helpful comments and suggestions. The authors were supported by Secyt-UNC grant (Res. Secyt 313/2016.), Argentina. The first author was partially supported by CIEM-CONICET, Argentina.


  • (1) Achim A, Kuruoglu EE, Zerubia J., SAR image filtering based on the heavy-tailed rayleigh model. IEEE Transactions on Image Processing 15(9):2686-2693 (2006).
  • (2)

    Allende H, Galbiati J. A non-parametric filter for digital image restoration, using cluster analysis. Pattern Recognition Letters 25(8):841-847 (2004).

  • (3) Allende H, Galbiati J, Vallejos R. Digital image restoration using autoregressive time series type models. Bulletin European Spatial Agency 434:53-59 (1998).
  • (4) Allende H, Galbiati J, Vallejos R. Robust image modeling on image processing. Pattern Recognition Letters 22(11):1219-1231 (2001).
  • (5) Aysal TC, Barner KE. Quadratic weighted median filters for edge enhancement of noisy images. IEEE Transactions on Image Processing 15(11):3294-3310 (2006).
  • (6) Baran S, Pap G, van Zuijlen MC. Asymptotic inference for a nearly unstable sequence of stationary spatial ar models. Statistics and probability letters 69(1):53-61 (2004).
  • (7) Basu S, Reinsel GC. Properties of the spatial unilateral first-order arma model. Advances in applied Probability 25(3):631-648 (1993).
  • (8) Bhandari A, Kumar A, Singh G. Improved knee transfer function and gamma correction based method for contrast and brightness enhancement of satellite image. AEU-International Journal of Electronics and Communications 69(2):579-589 (2015).
  • (9) Bustos O. Robust statistics in sar image processing. In: Image Processing Techniques, First Latino-American Seminar on Radar Remote Sensing, vol 407, p 81 (1997).
  • (10) Bustos, Oscar and Fraiman, Ricardo and Yohai, Victor J. (1984). Asymptotic behaviour of the estimates based on residual autocovariances for ARMA models. Robust and Nonlinear Time Series Analysis, Springer, pp. 26-49.
  • (11) Bustos O, Ojeda S, Vallejos R. Spatial arma models and its applications to image filtering. Brazilian Journal of Probability and Statistics pp 141-165 (2009).
  • (12) Bustos OH, Yohai VJ. Robust estimates for arma models. Journal of the American Statistical Association 81(393):155-168 (1986).
  • (13) Bustos OH, Ruiz M, Ojeda S, Vallejos R, Frery AC. Asymptotic behavior of ra-estimates in autoregressive 2d processes. Journal of Statistical Planning and Inference 139(10):3649-3664 (2009).
  • (14) Comport AI, Marchand E, Chaumette F. Statistically robust 2-d visual servoing. IEEE Transactions on Robotics 22(2):415-420 (2006).
  • (15) Gottardo R, Raftery AE, Yee Yeung K, Bumgarner RE. Bayesian robust inference for differential gene expression in microarrays with multiple samples. Biometrics 62(1):10-18 (2006).
  • (16) Hamza AB, Krim H. Image denoising: A nonlinear robust statistical approach. IEEE transactions on signal processing 49(12):3045-3054 (2001).
  • (17) Huang HC, Lee TC. Data adaptive median filters for signal and image denoising using a generalized sure criterion. IEEE Signal Processing Letters 13(9):561-564 (2006).
  • (18) Huber PJ, et al. Robust estimation of a location parameter. The Annals of Mathematical Statistics 35(1):73-101 (1964).
  • (19)

    Ji Z, Huang Y, Xia Y, Zheng Y. A robust modified gaussian mixture model with rough set for image segmentation. Neurocomputing (2017).

  • (20) Kashyap RL, Eom KB. Robust image modeling techniques with an image restoration application. IEEE Transactions on Acoustics, Speech, and Signal Processing 36(8):1313-1325 (1988).
  • (21) Kim JH, Han JH Outlier correction from uncalibrated image sequence using the triangulation method. Pattern recognition 39(3):394-404 (2006).
  • (22) Lin TC. A new adaptive center weighted median filter for suppressing impulsive noise in images. Information Sciences 177(4):1073-1087 (2007).
  • (23) Maronna R, Martin RD, Yohai V. Robust statistics. John Wiley and Sons, Chichester. ISBN (2006)
  • (24) Muler N, Pena D, Yohai VJ, et al. Robust estimation for arma models. The Annals of Statistics 37(2):816-840 (2009).
  • (25) Ojeda S. Robust RA estimators for bidimensional autoregressive models. PhD thesis, Universidad Nacional de Córdoba. Facultad de Matemática, Astronomı́a y Fı́sica (1999).
  • (26)

    Ojeda S, Vallejos R, Bustos O. A new image segmentation algorithm with applications to image inpainting. Computational Statistics and Data Analysis 54(9):2082-2093 (2010).

  • (27) Ojeda SM, Vallejos RO, Lucini MM. Performance of robust ra estimator for bidimensional autoregressive models. Journal of Statistical Computation and Simulation 72(1):47-62 (2002).
  • (28)

    Ojeda SM, Vallejos RO, Lamberti PW. Measure of similarity between images based on the codispersion coefficient. Journal of Electronic Imaging 21(2):023019 (2012).

  • (29) Ojeda SM, Britos GM, Vallejos R. An image quality index based on coefficients of spatial association with an application to image fusion. Spatial Statistics 23:1-16 (2018).
  • (30) Pistonesi S, Martinez J, Ojeda SM, Vallejos R. A novel quality image fusion assessment based on maximum codispersion. In: Iberoamerican Congress on Pattern Recognition, Springer, pp 383-390 (2015).
  • (31) Ponomarenko N, Lukin V, Zelensky A, Egiazarian K, Carli M, Battisti F. Tid2008-a database for evaluation of full-reference visual quality assessment metrics. Advances of Modern Radioelectronics 10(4):30-45 (2009).
  • (32) Ponomarenko N, Jin L, Ieremeiev O, Lukin V, Egiazarian K, Astola J, Vozel B, Chehdi K, Carli M, Battisti F, et al. Image database tid2013: Peculiarities, results and perspectives. Signal Processing: Image Communication 30:57-77 (2015).
  • (33)

    Prastawa M, Bullitt E, Ho S, Gerig G. A brain tumor segmentation framework based on outlier detection. Medical image analysis 8(3):275-283 (2004).

  • (34) Quintana C, Ojeda S, Tirao G, Valente M. Mammography image detection processing for automatic micro-calcification recognition. Chilean Journal of Statistics (ChJS) 2(2) (2011).
  • (35) Sadabadi MS, Shafiee M, Karrari M. Two-dimensional arma model order determination. ISA transactions 48(3):247-253 (2009).
  • (36) Schabenberger O, Gotway CA. Statistical methods for spatial data analysis. Text in Statistical Sciences. CHAPMAN and HALL/CRC press (2005).
  • (37)

    Singh M, Arora H, Ahuja N. A robust probabilistic estimation framework for parametric image models. In: European Conference on Computer Vision, Springer, pp 508-522 (2004).

  • (38) Tarel JP, Ieng SS, Charbonnier P. Using robust estimation algorithms for tracking explicit curves. Computer VisionECCV 2002 pp 492-507 (2002).
  • (39) Theodoridis S, Koutroumbas K. Pattern recognition. Elsevier (2003).
  • (40) Tjøstheim D. Statistical spatial series modelling. Advances in Applied Probability 10(1):130-154 (1978).
  • (41) Vallejos RO, Garcı́a-Donato G. Bayesian analysis of contaminated quarter plane moving average models. Journal of Statistical Computation and Simulation 76(2):131-147 (2006).
  • (42) Vallejos RO, Mardesic TJ. A recursive algorithm to restore images based on robust estimation of nshp autoregressive models. Journal of Computational and Graphical Statistics 13(3):674-682 (2004).
  • (43) Wang Z, Bovik AC. Image and multidimensional signal processing-a universal image quality index. IEEE Signal Processing Letters 9(3):81-84 (2002).
  • (44) Whittle P. On stationary processes in the plane. Biometrika pp 434-449 (1954).
  • (45) Yao Q, Brockwell PJ. Gaussian maximum likelihood estimation for arma models ii: spatial processes. Bernoulli 12(3):403-429 (2006)
  • (46) Zielinski J, Bouaynaya N, Schonfeld D. Two-dimensional arma modeling for breast cancer detection and classification. In: Signal Processing and Communications (SPCOM), 2010 International Conference on, IEEE, pp 1-4 (2010).