A Comparative Study of Some Markov Random Fields and Different Criteria Optimization in Image Restoration

The present chapter illustrates the use of some recent alternative methods to deal with digital image filtering and restoration. This collection of methods is inspired on the use of Markov Random Fields (MRF), which introduces prior knowledge of information that will allow, more efficiently, modeling the image acquisition process. The methods based on the MRF are analyzed and proposed into a Bayesian framework and their principal objective is to eliminate noise and some effects caused by excessive smoothness on the reconstruction process of images which are rich in contours or edges. In order to preserve object edges into the image, the use of certain convexity criteria into the MRF is proposed obtaining adequate weighting of cost functions in cases where discontinuities are remarked and, even better, for cases where such discontinuities are very smooth.

Some image analysis and processing tasks involve the filtering or image recovery x (e.g. restoration) after x passes by a degradation process giving the observed image y (see equation (1)). Such degradation covers the noise perturbations, blurring effects due to focus of the acquisition system lenses or to the bandwidth capacity, and other artifacts that may be added to the correct data. The use of Bayesian methods proposed in the seventies (Besag, 1974;1986;Geman and Geman, 1984), are nowadays essential at least in the cases of image segmentation and image restoration (Andrews and Hunt, 1977). The basic idea of these methods is to construct a Maximum a posteriori (MAP) estimate of the modes or so called estimator of true images by using MRF into a Bayesian framework. The general approach is based on a robust scheme which could be adapted to reject outliers, tackling situations where noise is present in different forms during the acquisition process (Bertaux et al., 2004;Cai et al, 2008;Chan et al., 2006;Durand and Nikolova, 2006a;b;Nikolova, 2006;Nikolova and Ng, 2005;Pan and Reeves, 2006).
The restoration or recuperation of an image to its original condition given a degraded image, passes thus by reverting the effects caused by noise or / and a distortion function. In fact, the degradation characteristic is a crucial source of information and its mathematical

A Comparative Study of Some Markov
Random Fields and Different Criteria Optimization in Image Restoration function must be known or estimated during the inversion procedure. Typically, this is a point spread function (PSF) which can be linked with the probability distribution of the noise contamination or likelihood function. In the case of MAP filters, usually the additive Gaussian noise is considered, however in some applications this noise is non-Gaussian (Bertaux et al., 2004;Cai et al, 2008) or unknown noise probability distribution function (pdf), p(e)=p(y|x), with some partial knowledge. In the latter case, we propose to use the data itself to obtain a non-parametric entropy estimate (EE) of the log-likelihood pdf, p n,h (e) (De la Rosa and Fleury, 2002;De la Rosa et al., 2003;De la Rosa, 2004). Then the log-likelihood will be optimized (e.g. minimized) together with a log-MRF to obtain the MAP image estimation. A variety of applications in signal processing and instrumentation, are based in statistical modeling analysis. One of the most used is the linear regression model (simple or multi-variable in the case of images) where y represents the response (observed data, acquired data or explained variables), to x explicative variables for i = 1,...,N and j = 1,...,M, and a system response parameterized by θ (being θ a linear operator which models the image degradation (or PSF)) which is associated to data (y, x). In some applications θ are functional parameters which will be estimated by an identification procedure if x are known, but if θ are known and x are unknown, the estimation is made about x, or the estimation can be made for both cases (e.g. blind deconvolution). The e variables are independent random processes identically distributed accordingly to p(e). A natural extension of the linear regression model is the non-linear regression model based on a parameterized function f (x, θ) (see equation (2)). This function is nonlinear with respect to the parameters, and its use is also considered because it has been shown in a large variety of signal processing and control applications that the modeling could be more realistic when using nonlinear functions. The perturbations affecting the analyzed system are also modeled as stochastic processes. Then, for this case: Another source of information which imposes a key rule in the image processing context, is the contextual or spatial information that represents the correlation between the intensity values of a neighborhood of pixels well specified. The modellig, when using MRF, takes into account such spatial interaction and it was introduced and formalized in (Besag, 1974), and extended in some pioneering works (Besag, 1986;Bouman and Sauer, 1993;Geman and Geman, 1984;Sauer and Bouman, 1992). Combining both kinds of information into a statistical framework, the restoration is led by an estimation procedure given the MAP of the true images when the distortion functionals are known. The algorithms implemented in this chapter were developed considering a degraded signal, where the resulting non-linear recursive filters show excellent characteristics to preserve all the details contained in the image and, on the other hand, they smooth the noise components. Some existent optimization methods are commented, since the optimization (e.g. maximization or minimization) task is central during the reconstruction and its performance depends on the structure of the cost functionals. Finally, a comparison between some existent Bayesian schemes and new proposed for noise filtering and image deblurring is presented. Particularly, five MAP estimation schemes are implemented where the first two schemes use two MRF; the semi-Huber (SH) potential A Comparative Study of Some Markov Random Fields and Different Criteria Optimization in Image Restoration 3 function which is proposed as an extension of some previous works (De la Rosa and Fleury, 2006;De la Rosa et al., 2007), and the generalized Gaussian MRF (GGMRF) introduced in (Bouman and Sauer, 1993) in both cases, Gaussian noise is assumed. The other three MAP schemes correspond to new MAP entropy estimators (MAPEE) for image filtering, where unknown noise pdf assumptions are made and the generalized Gaussian MRF is also used.
The rest of this chapter is organized as follows. Section 2 describes the general definition of a MRF and the proposal of the MAP estimator. The potential functions must be proposed to conduct adequately the inversion process, some functions are described in section 3 where the convexity level is the key to formulate an adequate criterion to be maximized or equivalently minimized. In sections 4 and 5 are deduced some MAP estimators and commented some optimization techniques resulting from different MRF structures for Gaussian noise assumptions, and MAPEE estimators for unknown noise pdf. Section 6 introduces the Kernel structures used for MAPEE estimators, while in section 7 results for Gaussian and non-Gaussian restoration are discussed. Finally, some conclusions and comments are given in section 8.

MAP estimation and Markov random fields
The problem of image estimation (e.g. restoration) into a Bayesian framework deals with the solution of an inverse problem, where the estimation process is carried out in a whole stochastic environment. The most popular estimators used nowadays are:

Maximum Likelihood estimator (ML):
A classical procedure to estimate x when θ is known (from equations (1) and (2)), is based in a cost function or criterion J (x) which varies in function ψ(·) of residuals or noise e(x),where: and so This is, for example, the case of the maximum likelihood (ML) estimator: The ML estimator is optimal when all information about the distribution p(e) is accessible, which is the case for Gaussian noise. When the knowledge about p(e) is imprecise or wrong, the estimator x ML is possibly suboptimal. Moreover, under certain circumstances, in image processing restoration, it results in an ill-posed problem or produces excessive noise and also causes smooth of edges. The regularization of the ML estimator gives both a more effective approach to reduce noise and smoothness (in fact, it depends on the mathematical structure of the regularization term). In a probabilistic framework, the regularization of the ML estimator leads to the Bayesian approach, where it is important In this case, the estimator is regularized by using a MRF function g(x) which models all prior information as a whole probability distribution, where X is the set of pixels x capable to maximize p(x|y),andp(y|x) is the likelihood function from y given x.
The Markov random fields (MRF) can be represented in a generic way by using the following equation: where Z is a normalization constant, C is a set of "cliques" c or local neighborhood of pixels, and V c (x) is a weighting function given over the local neighborhood. Generally, the "cliques" correspond to the possible set of neighborhoods of pixels. A theorem introduced by Hammersley-Clifford (Besag, 1974;Geman and Geman, 1984) probes the equivalence between the Gibbs distribution and the MRFs. The Markov random fields have the capacity to represent various images sources. The main partial disadvantage using MRFs is that the estimation procedure deals with local minimization schemes (generally it grows the computation time), while the proposal of global minimization schemes has been recently introduced by modifying the local structure of the MRFs in order to establish a better trade-off between quality of reconstruction and computation time. There is a variety of MRF models which depend on the cost functions also known as potential functions that can be used. Each potential function characterizes the interactions between pixels in the same local group. As an example, the following family represents convex functions: where Δ = λ[x s − x r ], λ is a constant parameter to be chosen, and p takes constant values such as p ≥ 1 accordingly to the theorem 2 in (Bouman and Sauer, 1993).

MRFs and convexity of potential functions
In recent works (Champagnat and Idier, 2004;Ciuciu and Idier, 2002;Ciuciu et al., 2001;Giovannelli et al., 2002;Idier, 2001;Rivera and Marroquin, 2003;2004;Rivera, 2005) some new potential functions were introduced. Such proposed functions are semi-quadratic functionals or half-quadratic and they characterize certain convexity into the regularization term (Geman and Reinolds, 1992;Geman and Yang, 1995)(eg. extension of penalization) which permits to build efficient and robust estimators in the sense of data preservation which is linked to the original or source image. Also, the necessary computation time decreases with respect to other proposed schemes as shown in (Chan et al., 2006;Durand and Nikolova, 2006a;b;Nikolova, 2006;Nikolova and Ng, 2005;Nikolova and Chan, 2007) and  (Labat and Idier, 2006;Labat, 2006) in the case of non-convex and convex potential functions.
On the other hand, in previous works from (Gibbs, 2000) it has been proposed a way to obtain the posterior distribution of the images; in this case, it is necessary to use sophisticated stochastic simulation techniques based on the Markov Chain Monte Carlo (MCMC) methods (Neal, 1993;Robert and Casella, 2004). If it is possible to obtain the posterior distribution of any image, thus, it is also possible to sample from such posterior distribution and obtain the MAP estimator, or other estimators such as the median estimator which is sensible to be as well as the MAP estimator under certain circumstances imposed by the noise structure, the MAP and the median estimators search the principal mode of the posterior distribution.
In the present chapter two already reported potential functions are revised and their performance is compared. The generalized Guassian MRF introduced in (Bouman and Sauer, 1993;Sauer and Bouman, 1992) is compared with respect to semi-Huber functional used in problems of one dimensional estimation in (De la Rosa and Fleury, 2006).

Generalized Gaussian MRF (GGMRF)
If one considers to generalize the Gaussian MRF (when p = q = 2 one has a GMRF, see equation (15)) as proposed by Bouman, where the generalized potential functions correspond to where ct is a constant term, and theoretically a s > 0a n db sr > 0, s is the site or pixel of interest and S is the set of sites into the whole MRF, and r corresponds to the local neighbors. In practice and for Gaussian noise assumptions, it is recommended to take a s = 0, assuring the unicity of x MAP , thus the GGMRF can be written as and from equation (2), log p(y|x) is strictly convex and so x MAP is continuous in y,a n dinp. On the other hand, the MAPEE approach proposed here takes into account that the noise pdf is unknown or non-Gaussian. For this particular case, it is recommended to take a s = 1, since the likelihood term is not given in terms of quadratic q = 2 functional, and in order to relax the possible non-convexity problem it has been used

147
A Comparative Study of Some Markov Random Fields and Different Criteria Optimization in Image Restoration www.intechopen.com 6 Image Acquisition also ensuring that log p(y|x) is convex and so x MAPEE is continuous in y,andinp. The choice of p is capital, since it constrains the convergence speed of the local or global estimator, and the quality of the restored image.

Semi-Huber potential function
In order to assure completely the robustness into the edge preserving image filtering, diminishing at the same time the convergence speed, the Huber-like norm or semi-Huber (SH) potential function is proposed as a half-quadratic (HQ) function. Such functional has been used in one dimensional robust estimation as described in (De la Rosa and Fleury, 2006) for the case of non-linear regression. This function is adjusted in two dimensions according to the following equation: where and where Δ 0 > 0 is a constant value, b sr is a constant that depends on the distance between the r and s pixels, and ϕ 1 (x)=e 2 where e =(x s − x r ). The potential function ρ 2 (x) fulfills the following conditions

MAP estimators for Gaussian noise and optimization
In this section some estimators are deduced. The single problem of filtering noise to restore an observed signal y leads to establish the estimators.

Image filtering and optimization
The observation equation could be where I is the identity matrix. For this particular problem when using the GGMRF and under hypothesis of Gaussian noise with variance σ 2 n ,theMAPestimatorisgivenby , the minimization in equation (15) or optimization problem leads to consider various methods (Allain et al., 2006;Chan et al., 2006;Nikolova, 2006;Nikolova and Ng, 2005;Nikolova and Chan, 2007): • global iterative techniques such as: the descendent gradient (Nikolova and Chan, 2007), conjugate gradient (Rivera and Marroquin, 2003) (for recent propositions one can consult the work (Labat and Idier, 2006;Labat, 2006)), Gauss-Seidel, Over-relaxed methods, etc. • local minimization techniques: minimization at each pixel x s (which generally needs more time, but from our point of view are more precise), where also some of the above methods can be used.
In this work the local techniques were used, the expectation maximization (EM) was not implemented, nor complete half-quadratic scheme as proposed in (Geman and Reinolds, 1992) and in (Geman and Yang, 1995), since all hyper-parameters included into the potential functions were chosen heuristically or according to values proposed in references. Only, the step of minimization with respect to x was implemented to probe convergence of estimators. For example, the local MAP estimator for the GGMRF is given by where, according to the value of parameters p and q, the performance of such estimator varies. For example, if p = q = 2, one has the Gaussian condition of the potential function where the obtained estimator is similar to the least-squares one (L 2 norm) since the likelihood function is quadratic, with an additional quadratic term of penalization which degrades the estimated image On the other hand, in the case of p = q = 1, the criterion is absolute (L 1 norm), and the estimator converges to the median estimator, which in practice, is difficult to implement in a precise way x s = median (y s , x r 1 ,...,x r I ) .
This criterion is not differentiable at zero and this fact causes instability in the minimization procedure. For intermediate values of p and q these estimators become sub-optimal, and the iterated methods can be used to minimize the obtained criterions. Such iterative methods are the sub-gradient, or the Levenberg-Marquardt method of MATLAB 7, the last method was used as optimization technique in this work. For cases where q = p,f o re x a m p l eq = 2 and 1 < p < 2, some studies and different propositions of a priori functions have been proposed in (Durand and Nikolova, 2006a;b;Nikolova, 2006;Nikolova and Ng, 2005;Nikolova and Chan, 2007), particulary in (Durand and Nikolova, 2006a;b;Nikolova, 2006;Nikolova and Ng, 2005) where non-convex regularized least-squares schemes are deduced and its convergence is analyzed. 2) moreover, if p = q, the criterion is not homogeneous, but: x(αy, λ)=α x(y, α 1−q/p λ) is convex, assuring the convergence and existence of the estimator which is continuous with respect to p.
Moreover, when the noise is Gaussian, the value of q = 2 is a good choice, but if the noise is non-Gaussian, and if the structure of noise is known, thus the likelihood term changes giving a particular estimator such as that proposed in (Bertaux et al., 2004), or in (Cai et al, 2008), giving some properties of robustness to the minimization procedure. If the structure of noise is supposed unknown, still one could reconsider the GGMRF with values for q ∈ (1, 2], or consider also another type of MRFs as those described in some works of Idier (Champagnat and Idier, 2004;Ciuciu and Idier, 2002;Ciuciu et al., 2001;Giovannelli et al., 2002) and Nikolova (Nikolova, 2006;Nikolova and Ng, 2005;Nikolova and Chan, 2007).
A second MAP estimator can be obtained when using the SH potential function, the global estimator can be described by the equation As in the previous MAP estimator, it has been proposed to implement the local estimator which leads to a similar expression for the first local MAP estimator (see equation (16)), that is The use of a prior distribution function based on the logarithm, with any degree of convexity and quasi-homogeneous, permits to consider a variety of possible choices of potential functions. Maybe, the most important challenges that must be well solved are: the adequate selection of hyper-parameters from potential functions, where different versions of the EM algorithms try to tackle this problem (Champagnat and Idier, 2004;Giovannelli et al., 2002;Idier, 2001;Nikolova and Ng, 2005), another is the minimization procedure which in any sense will regulate the convergence speed as proposed in (Allain et al., 2006;Geman and Reinolds, 1992;Geman and Yang, 1995;Labat and Idier, 2006;Labat, 2006;Nikolova, 2006;Rivera and Marroquin, 2003), etc.

Image deconvolution
On the other hand, for the problem of image deblurring to restore an observed signal y,t h e observation equation used is given by Now, for the two previous MAP estimators the likelihood term changes, such that, 150 Advanced Image Acquisition, Processing Techniques and Applications www.intechopen.com where the matrix Θ is known and given by the following truncated Gaussian blurring function, as used also in (Nikolova et al., 2010), with σ b = 1.5. And k = 1, 2 according to the GGMRF, and the SH potential functions. Here, the results were improved combining ideas introduced in a similar Bayesian way by Levin (Levin et al., 2007a;b) adding a Sparse prior (SP) for filtering and then reconstructing the image.

MAP entropy estimators (EE) for unknown noise pdf
Our proposition for a new MAP scheme is to use a generalized Gaussian MRF introduced in section 3.1, together with three kernel estimators used in (De la Rosa and Fleury, 2002;De la Rosa et al., 2003;De la Rosa, 2004) to obtain cost functionals or criterions based on the entropy of the approximated likelihood function (first term of equation (6)) p n,h (e).T h u s , − log p(y|x) is built on the basis of the entropy of an estimate (EE) version p n,h (e) of the distribution p(e). A first proposition is due to Pronzato and Thierry (Pronzato and Thierry, 2000a;b;2001). The approximation is obtained using the classical kernel estimators which uses the empirical distribution of the random vector e 1 (x),...,e n (x). The next expression denote such estimators: This expression assumes the hypothesis that p(e) is symmetric, two times differentiable and positive, indeed, it is assumed that K(z) is a kernel weighted function which satisfies some imposed conditions treated in (Masry, 1983) and subsequently taken back by Devroye (Devroye, 1992;1989;Devroye and Krzyzȧk, 1999;Devroye, 1997), Berlinet (Berlinet and Devroye, 1994), and Loader (Loader, 1999) in some of their research work. The bandwidth h = h n is given in function of the sample size n. This parameter could be considered as a sequence of positive numbers that must satisfy: h n → 0andnh n → ∞ when n → ∞. The strong uniform consistency of p n,h (e) and its convergence toward p(e), depend on a convenable procedure of bandwidth selection. A simple and faster procedure which has been retained in this work is the technique proposed and developed by Terrell (Terrell, 1990;Terrell and Scott, 1985). In the two dimensional kernel cases, the previous idea has been extended in this work accordingly to the equation: If the convergence and consistence of p n,h (e) are assumed, such that p n,h (e) → p(e), then the entropy criterion over p n,h (e) can be approximated to − log p(y|x). The fact that the entropy of any probability density function is invariant by translation, leads to consider one practical artifact to build a suitable criterion. An extended criterion p n,h (e E ) is based on the residuals or noise extended vector which is given by: e E = {e 1,1 (x),...,e n,n (x), −e 1,1 (x),...,e n,n (x)} where Finally, if we assume that the EE is a version of the log-likelihood function into the MAP estimator, then a general version of the MAP-entropy estimator (MAPEE) which assumes unknown noise pdf can be constructed from the fact that − log p(y|x) can be approximated by the entropy of an estimate version p n,h (e) of the distribution p(e),thatisH A p n,h (e E ) ,thus: Now, substituting a particular log g(x) such as the GGMRF into the equation (27) one could obtain at least three MAPEE estimators given by for m = 1, 2, 3 according to the three following kernels: Normal, cosine and Hilbert.

Af u n c t i o no ft h ef o r mK(z) is assumed as a fixed kernel K h (z)=1/(h d )K(z/h),w h e r e
h is called the kernel bandwidth and h > 0. The fundamental problem in kernel density estimation lies in both the selection of an appropriate value for h and the selection of the kernel structure. The choice of K(z) could depend on the smoothness of p(e). Three different kernels or nonparametric schemes are reviewed in this section to approximate p n,h (e E ).T h e gaussian kernel, which has proved to give good performance when h is selected by using the over-smoothed principle introduced by Terrell (Terrell, 1990;Terrell and Scott, 1985) when the errors vector is generally of finite length n. Second kernel is the cosine based weights functions, where h is viewed in a different way. Finally, the third kernel is obtained from a recent class of Hilbert kernels (Devroye and Krzyzȧk, 1999). It avoids the bandwidth h selection and its performance depends on other parameters, whose selection is very easy (parameters d and k are defined in section 6.3).

Normal or gaussian kernel
Among the different classical kernels (Berlinet and Devroye, 1994), the gaussian kernel has been chosen since it is one of the most popular and easy to implement estimator. The following expression resumes this estimator by a sum of exponentials: In such a case, and considering that a fixed kernel structure has been chosen, Terrell (Terrell, 1990) proposes to use an over-smoothed bandwidth h that corresponds to: this bandwidth value guarantees the minimization of the Mean Integrated Squared Error (MISE), σ is the standard deviation of the sample z,a n d K(z) 2 dz = 1 (2 . Under mild conditions, the kernel density estimates based on the over-smoothing principle are consistent and for sufficiently large sample sizes, they will display all information present in the underlying errors density p(e),wechosen × n array of samples z k,l near z.

Cosine based weights functions
A sequence of special weight cosine functions c q (z)=1/(A q ) cos 2q (z/2) requires only O(q 2 n) arithmetic operations, where q is a slowly growing function that increases without bound together with n (i.e. q(n) < n or q(n) > n). This estimator is similar to the kernel classical estimators and evaluates a series expansion in a more efficient way. The role of h is replaced by q (c q (z) ≡ K h (z)) and the efficiency is attained (according to the time of calculus), since the selection of q is generally more simple than h. Thus (25) is equivalent to p n,q (z)= 1 where the A q value could be approximated by Sufficient consistency and convergence conditions are given in (Egecioglu and Srinivasan, 2000) for the one dimensional case.

The Hilbert kernel
Finally, other kernel density estimates called the Hilbert kernel estimate is used. The equation where the smoothing factor h is canceled obtaining: The Hilbert estimates are viewed as a universally consistent density estimate whose expected performance (L 1 , L ∞ , pointwise) is monotone in n (at least in theory) for all densities. The consistency of this class of estimators is proved in Devroye and Krzyzȧk (1999) infinite peaks produced during estimation, in one dimensional case and using the value of k = 2 the kernel estimate is given by: where Den i,j = z − z i 2d + z − z j 2d and V d is the volume of the unit ball in R b .T h i s last expression is also called Cauchy density estimate, due to its similarity to the multivariate Cauchy density, · denotes the L 2 metric on R d . Finally, it is assumed that p n (z) → p(z) at least in probability for almost all z. For a suitable choice of A n and alternatively of h n , q,o rd and k, these estimators could be "blind asymptotically efficient". The asymptotic properties and the strong consistency of the truncated entropy estimators were analyzed in (Pronzato and Thierry, 2000a;b;2001). More over, in some other papers, the power of these nonparametric tools have been largely used for different signal processing problems , (Wolsztynski et al., 2004a) - (Wolsztynski et al., 2005c).

Denoising and deblurring experiments
Results presented in this section were obtained experimenting extensively with four images: synthetic, Lena, Cameraman, and Boat (shown in Figure 1), to probe the performance of the presented estimators.

Image filtering for Gaussian noise
Some estimation results are presented when images are contaminated by Gaussian noise, and there are no other type of distortions (all θ i,j = 1). The first experiment was made considering that σ n = 5, 10, 15. In the work (De la Rosa et al., 2003) some results were previously presented based on the analysis of a synthetic image and the standard image of "Lena", different levels o fn o i s ew e r ea d d e dt ot h ei m a g e s : e ∼N(0, Iσ 2 n ),t h ev a l u e so fσ n are given such that the obtained degradation is perceptible and difficult to eliminate. The obtained results were compared using different values for p and λ preserving q = 2 (MAP1), and different values for Δ 0 and with λ = 1 (MAP2). Generally, with the two MAP estimators the filtering task gives good visual results preserving contours and smoothing noise (see Figures 2 and 3), but the computation time is different between them, the fastest estimator in general is the MAP2 (see image 2(d)), while the slowest is the MAP1 with p = 1.2. Table 1 shows the performance of the MAP estimators for the problem of filtering Gaussian noise for different images, where an objective evaluation is obtained accordingly to the peak signal to noise ratio (PSNR). Also the computation times in MATLAB are shown in Table 1. Such comparative evaluation shows that our proposed approach MAP2, gives better or similar performance with respect to MAP1. Some interesting applications of robust estimation are particulary focused in phase recovery from fringe patterns as presented in the work Villa et al. (2005), phase unwrapping, and some other problems in optical instrumentation. In this sense, some filtering results were thus obtained using the presented MAP estimators with similar results to those reported here. On the other hand, the use of half-quadratic potential functions permits flexibility for the computation time (Durand and Nikolova, 2006a;b;Nikolova and Ng, 2005). But still is a challenge to tune correctly the hyper-parameters to obtain a better performance in the sense 154 Advanced Image Acquisition, Processing Techniques and Applications www.intechopen.com  of restoration quality, the most simple potential function to tune is the semi-Huber (MAP2). Also, making the correct hypothesis over the noise could help to improve the performance of the estimator. This could be directly reflected by proposing a more adapted likelihood function, as proposed in section 5 for cases of non-Gaussian noise or unknown pdf. More over, also the problem could be solved by using variational and partial differential equations as illustrated in the famous work of Perona and Malik (Perona and Malik, 1990), and some recent related works.

Image deconvolution
Now, for the problem of image deconvolution some estimation results are presented when images are contaminated by Gaussian noise, and Gaussian distortion (with σ b = 1.5) blurring the image. This second experiment was carried out considering that σ n = 3, 5, 7. The results were also compared using different values for Δ 0 and with λ = 1 (MAP2), different values for p and λ preserving q = 2 (MAP1). Figure 4 shows a comparison of results obtained for the Cameraman image accordingly to the two MAP estimators. One can notice that preserving values of hyper-parameters near those for the filtering case, the estimators smooth the noise but does not made a good recuperation of edges into the image. One must change the hyper-parameters values searching a trade of between the granularity of the noise and the sharpness of the image. Figure 5 shows the results obtained on the Cameraman image when using a combination of proposed estimators together with a Sparse prior (SP) deconvolution technique introduced in (Levin et al., 2007a;b), the improvement in restoration is visible. In the first stage only the SP deconvolution is obtained (see image 5(b)), where some ringing and granularity artifacts occur due to noise and they cannot be correctly smoothed. Then, applying first the filtering of noise and then the SP deconvolution, the images 5(c) and 5(d) are obtained. In Table 2 it is shown the performance of the MAP estimators together with the SP deconvolution for three images giving similar results. An objective evaluation is made using the PSNR measure and also computation times in MATLAB are shown.

Image filtering for non-Gaussian noise
Now, for the problem of filtering non-Gaussian noise, some estimation results were obtained when images are contaminated by Gamma, Beta, Uniform and impulsive noise. Here, it is assumed that there are not other type of distortions (all θ i,j = 1). The observation equation then could be y = x + e,w h e r ee ∼G(α, β), e ∼B(α, β),...
A third experiment was made considering Gamma noise where α = 1.5, 3.5, 7 and β = 1, 2, 3, and also two factors of amplification of noise were used σ a = 5, 10 (σ a G(α, β)). The values of α and β are given in such a way that the obtained degradation is perceptible and difficult to eliminate.    Figure 6 shows the filtering of the synthetical image using the three MAPEE estimators introduced in section 5 for parameters λ = 10, p = 1.2 of the GGMRF. Generally, with the three estimators the filtering task gives good performing results (see also Figure 7 where a three dimensional-3D view is shown), but the computation times are different between them, the fastest estimator was MAPEE3 (image 6(d)). Some objective measurements were obtained for all probe images, such as the PSNR depicted in Table 3 which agree with the visual results. Also all times of computation are shown in Table 3. On the other hand, in the filtering case of the Cameraman and the Boat images the results coincide with those obtained for synthetic and Lena images. Visually one can see in Figures 8 and 9 the robustness performance of MAPEE estimators. Some other experiments were conducted for Beta and impulsive noise. For the case of Beta noise generation we choose α = 1.5, 2.5, 5, 7.5 and β = 1, and also two factors of amplification of noise were used σ a = 10, 20 (σ a B(α, β)). The filtering results are promising for non-Gaussian noise, the obtained results for the Beta noise filtering also confirm the robustness of the proposed approach. Moreover, as one can see for impulsive noise with a 15% of outliers in Figures 10 and 11, the robustness is assured when using the three MAPEE estimators, where for instance, in the synthetic image case the PSNR was improved from 11.45 dB to 17.1 dB.

Conclusions and comments
Some advantages on the use of GGMRF as prior information into the Bayesian estimation (MAP1) are: the continuity of the estimator is assured as a function of the data values when 1 < p ≤ 2 even for Gaussian and non-Gaussian noise assumptions. The edge preserving is also assured, over all when p → 1 and obviously it depends on the choice between the interval 1 < p < 2. The use of other semi-quadratic or half-quadratic potential functions such as SH (MAP2) also presents some advantages: the convexity of these functions is relaxed with respect to the GGMRF. This fact gives as a result the decrement of the computation time, which is a good advantage over the GGMRF, since the tuning of hyper-parameters in the GGMRF is more complicated, one has more degrees of freedom. However, this problem can be solved as argued by Idier (Champagnat and Idier, 2004) and Rivera (Rivera and Marroquin, 2003) by implementing more sophisticated algorithms. The robustness in both cases is assured. Moreover, the advantages presented by those estimators could overcome the disadvantages, by experimenting largely with them, and establishing the aim to implement more sophisticated algorithms with the compromise to reduce computation time and better quality in restoration as recently exposed by Ruimin Pan (Pan and Reeves, 2006).
On the other hand, the performance and improvement of the MAPEE estimators for non-Gaussian noise filtering, could be classified in terms of simplicity and in terms of filtering quality. The obtained results for the three proposed estimators are favorable in general in the sense of robustness, and the fastest convergence is obtained for the MAPEE1 and MAPEE3 estimators (in the case of the MAPEE2 estimator q = 40, and for MAPEE3, k = 2, d = 1, and V d = 0.7071. The value of q was changed alternatively in the range of 30 ≤ q ≤ 90, and some performance improvement has been noticed). A general scheme for MAPEE estimators has been introduced and it works in real frameworks, where the nonlinearity conditions of some systems could be present. For future works it exists the interest to implement procedures of Bayesian estimation into high level programming that will be characterized into algorithms to be used in DSP cards, and tasks such as image reconstruction and segmentation. Also, one can probe to use MAPEE to processing real signals issued from real instrumentation or signal processing problems. The final objective of this work has been to contribute with a series of software tools for image analysis, focused for instance, in optical instrumentation tasks such as those treated in the works (Villa et al., 2005)- (Villa et al., 2010b) and (Rivera and Marroquin, 2004;Rivera, 2005) obtaining competitive results in filtering and reconstruction.

Acknowledgements
We would like to express our deep gratitude to the Secretaría de Educación Pública (SEP) of Mexico, this work was partially supported through the PIFI Mexican Program (PIFI 2010-2011).

References
Allain, M.; Idier, J. and Goussard, Y. (2006). On global and local convergence of half-quadratic algorithms, IEEE Trans. Image Processing, Vol. 15, No. 5, pp. 1130-1142. Andrews, H. C. and Hunt, B. R. (1977. Digital image restoration, Prentice-Hall, Inc., New Jersey. "Advanced Image Acquisition, Processing Techniques and Applications" is the first book of a series that provides image processing principles and practical software implementation on a broad range of applications. The book integrates material from leading researchers on Applied Digital Image Acquisition and Processing. An important feature of the book is its emphasis on software tools and scientific computing in order to enhance results and arrive at problem solution.

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: José I. De la Rosa, Jesús Villa, Ma. Auxiliadora Araiza, Efrén González and Enrique De la Rosa (2012). A Comparative Study of Some Markov Random Fields and Different Criteria Optimization in Image Restoration, Advanced Image Acquisition, Processing Techniques and Applications I, Dr. Dimitrios Ventzas (Ed.), ISBN: 978-953-51-0342-4, InTech, Available from: http://www.intechopen.com/books/advanced-image-acquisitionprocessing-techniques-and-applications-i/a-comparative-study-of-some-markov-random-fields-and-differentcriteria-optimization-for-image-resto