Image magnification based on a blockwise adaptive Markov random field model

Image and Vision Computing 26 (2008) 1277–1284

Contents lists available at ScienceDirect

Image and Vision Computing journal homepage: www.elsevier.com/locate/imavis

Image magnification based on a blockwise adaptive Markov random field model Xiaoling Zhang a,b,c, Kin-Man Lam a,*, Lansun Shen b a b c

Centre for Signal Processing, Department of Electronic and Information Engineering, The Hong Kong Polytechnic University, Hong Kong, China Signal and Information Processing Lab, Beijing University of Technology, Beijing, China Department of Communication Engineering, Xiamen University, Xiamen, Fujian, China

a r t i c l e

i n f o

Article history: Received 31 May 2006 Received in revised form 7 January 2008 Accepted 4 March 2008 Available online Keywords: Image magnification MAP estimation Markov random field Iterated function systems

a b s t r a c t In this paper, an effective image magnification algorithm based on an adaptive Markov random field (MRF) model with a Bayesian framework is proposed. A low-resolution (LR) image is first magnified to form a high-resolution (HR) image using a fractal-based method, namely the multiple partitioned iterated function system (MPIFS). The quality of this magnified HR image is then improved by means of a blockwise adaptive MRF model using the Bayesian ‘maximum a posteriori’ (MAP) approach. We propose an efficient parameter estimation method for the MRF model such that the staircase artifact will be reduced in the HR image. Experimental results show that, when compared to the conventional MRF model, which uses a fixed set of parameters for a whole image, our algorithm can provide a magnified image with the well-preserved edges and texture, and can achieve a better PSNR and visual quality. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction In most electronic imaging applications such as scientific, medical and space systems, consumer electronics, and entertainment applications, images with a high resolution are desired and often required. A high-resolution (HR) image can not only give viewers a pleasing picture, but also provide additional details of the picture that are important for recognition and analysis in many applications. Given a low-resolution (LR) image, a magnification technique can provide a HR image with the visual contents of the original image preserved as much as possible, without changing the resolution of the image sensor. To compute an estimate of the HR image of a given LR observation, the commonly used image magnification approaches are based on nearest-neighbor, bilinear, and bicubic interpolations [1,2]. The nearest-neighbor interpolation copies from the corresponding neighboring pixels to an empty location in a HR image, which tends to produce blocky images. The bilinear and bicubic techniques usually smooth an image in discontinuous regions, which then produce a HR image with a blurred appearance. Some non-linear interpolation techniques [3,4] have also been proposed to maintain edge sharpness in recent years. Other more sophisticated solutions to the magnification problems include the Bayesian maximum a posteriori (MAP) approach [5], wavelet-based approach [6,7], fractal-based approach [8–13], and PDEs-based approach [14,15]. * Corresponding author. Tel.: +852 27666207. E-mail addresses: [email protected] (X. Zhang), [email protected] (K.-M. Lam). 0262-8856/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.imavis.2008.03.003

The conventional fractal-based method causes serious blocky artifacts due to the independent and lossy coding of the range blocks. The MAP approach in [5] can preserve the edges in an image using the Huber-MRF prior model, which is a constrained optimization problem with a unique minimum. The Bayesian MAP estimation method is also very promising in image super-resolution technology [7,16–18]. Therefore, our proposed algorithm will magnify a LR image in the Bayesian MRF framework. However, the conventional MRF model used is usually characterized by a single set of fixed parameters. This cannot describe the whole image effectively, as natural images are usually inhomogeneous. Consequently, a conventional MRF model will often produce blocky or staircase artifacts at the edge regions. These artifacts will become much more serious when an image is magnified 4  4 (16) times or more. Image magnification is an ill-posed problem, which can be solved using a Bayesian MAP framework with a probabilistic image model. This MAP approach can identify an optimal solution through the use of a prior image model. In other words, it can provide a flexible and convenient way to model prior knowledge. In this paper, we employ this approach to search for an optimal realization of the HR image given a LR image when the magnification factor is 16 or more, and we use the gradient-descent algorithm to compute the HR estimate. In addition, our MAP approach is different from the previous ones in the following two ways. First, although the MRF can provide a general model for images, the conventional MRF uses a single set of fixed parameters to model a whole image. Most of the parameter estimation methods for the MRF model are used for homogeneous images only, and the estimation of the parameters is usually computationally intensive. In

1278

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

order to apply the MRF model to characterize prior knowledge for image magnification, we propose a new, efficient parameter estimation method based on blockwise adaptive estimation. In our method, an initial estimated HR image is first divided into small blocks, and the parameters of the MRF model for each block are defined as inversely proportional to the energy in its corresponding direction. The particular MRF prior model is used to further regularize the estimated HR image in the next iteration. With this approach, we can avoid artifacts at the edges while retaining and enhancing sharpness in the HR images. Second, the initial HR image is calculated based on the multiple fractal codes, i.e. the multiple partitioned iterated function system (MPIFS) proposed in [13]. Using different range block partitioning schemes, a number of magnified images are generated, which are then averaged to obtain the initial HR image. The averaging process in the MPIFS method can remove the noise and preserve the edges well. In summary, our proposed method first employs the MPIFS to generate an initial HR image, whose quality is then improved by a regularization procedure with a block-based adaptive MRF model. Experiments show that our algorithm can produce HR images with significantly better visual quality and better PSNR quality than the conventional MRFMAP approach can. This paper is organized as follows. In the next section, we will present the LR image acquisition model and the Bayesian MAP analysis framework. In Section 3, we introduce the MRF prior model and describe our proposed parameter estimation method used for image magnification. Section 4 will give a brief introduction to the MPIFS magnification algorithm used for generating an initial HR image. Section 5 provides the implementation details of our algorithm in the MAP framework. Section 6 gives the experimental results, and Section 7 contains the discussion and conclusion. 2. LR image acquisition model and bayesian map approach In this section, we will describe the LR image formation model from its HR version, and show how image magnification is an optimization problem that can be solved using the Bayesian MAP approach. 2.1. The LR image formation

2.2. The maximum a posteriori (MAP) estimation The MAP estimate is a statistical approach which has the appealing attribute of yielding the most likely image given the observed data. This approach can work very well with a variety of problems. Given a LR image Y, the MAP estimate of the HR image X is obtained by maximizing the following conditional probability density function P(XjY): b ¼ arg max PðXjYÞ: X X

ð3Þ

Using the Bayesian rule and taking the logarithm of the a posteriori probability, (3) can be written as follows: b ¼ arg minf log PðYjXÞ  log PðXÞ þ log PðYÞg; X X

ð4Þ

where P(Y) is a constant because Y is the given LR observation, P(X) is the image prior model, and P(Y|X) is the probability of reconstructing the LR image Y given the HR image X. For noise-free cases, P(Y|X) is also a constant. Therefore, the MAP estimate of the HR image is given as follows: b ¼ arg minf log PðXÞg; X X

ð5Þ

where X is the set of all images that satisfy the constraints. In the MAP optimization (5), the a priori knowledge represented by P(X) provides a regularized HR image estimate. Thus, the specific choice of a prior distribution P(X) is critical. The Markov random field (MRF) has the capability of representing different images and the context-dependent nature of each image, so it is adopted to model the HR images. As a matter of fact, the MRF has been widely used to model images in the Bayesian framework for characterizing a stochastic pattern or prior knowledge. 3. The MRF prior model and parameter estimation

Suppose that X represents the HR image of size qN1  qN2, and Y represents the corresponding LR image of size N1  N2, where q is the down-sampling factor in both the horizontal and vertical directions. If we consider the magnification of the LR observation Y to the HR image X, q will be called the zooming factor or the magnification factor. In this paper, the corresponding block of q  q pixels obtained by magnifying a pixel in the LR image is called a zooming block. In a LR image acquisition procedure, the point spread function (PSF) of the LR sensor is usually modeled by spatial averaging [5], so the LR image Y is related to its HR version X, as described by (1). Here, we assume that no noise is introduced in the acquisition process. ! qðiþ1Þ1 X qðjþ1Þ1 X 1 Y i;j ¼ 2 X k;l ; i ¼ 0; . . . ; N 1  1; q ð1Þ k¼qi l¼qj j ¼ 0; . . . ; N 2  1: Eq. (1) defines the constraint of each pixel value Yi,j of the LR image, which should be equal to the average of the corresponding q  q pixels, i.e. the zooming block, in the HR image. The formulation (1) can be written in matrix form as follows: Y ¼ DX;

of possible HR images X for a given LR image Y, the image magnification is therefore an ill-posed inverse problem in the sense of Hadamard [5,7]. A well-posed problem can be formulated using a stochastic regularization technique by means of the Bayesian MAP estimation; this results in a constrained optimization problem with a unique minimum using a variety of established techniques.

ð2Þ

where D is the down-sampling matrix of dimension N1N2  qN1qN2, which generates aliased LR images. Since there is an infinite number

It is well known that the neighboring pixels in an image tend to have a strong correlation to each other, as the intensity values of the pixels usually vary gradually, except in the boundary regions. The MRF is a popular approach used in image processing to reflect this property of smoothness. The joint distribution of an MRF is characterized by a Gibbs function with a set of clique potential parameters, where the parameters specify the weight of smoothness in an image [19]. Typically, these parameters allow the MRF prior model to be adjusted for the best performance. In this section, we will first present the MRF model, and then propose an efficient blockwise parameter-estimation method. 3.1. The Markov random field (MRF) model In the 1920s, mostly inspired by the Ising model [20], a new type of stochastic process appeared in probability theory, namely the Markov random field (MRF). The use of MRF in image processing became popular after the work of Geman and Geman [21] on image restoration in 1984. The MRF theory provides a convenient and consistent way of modeling context-dependent entities such as image pixels and other spatially correlated features. The practical use of MRF models was largely ascribed to the equivalence between MRF and Gibbs distributions established by Hammersley and Clifford and further developed by Besag [22].

1279

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

Let S = {s = (i,j)|1 6 i 6 M1,1 6 j 6 M2} denote an M1  M2 pixel site lattice with a neighborhood system N = {Ns:s 2 S}, where Ns  S is a set that contains the neighboring sites of s. Let c 2 S be a clique, i.e. every site in the clique is a neighbor of all the other sites in the clique, and let C denote the set of cliques. Let X = {Xs:s 2 S} be a random field defined on S. The random field or the image X can be modeled as an MRF if PðX s jX Sjs Þ ¼ PðX s jX Ns Þ and PðXÞ > 0;

ð6Þ

where Sjs denotes the whole image except the site s. Ref. [19] has also shown that the joint probability distribution of X is equivalent to the Gibbs distribution. Gibbs distributions are used to explicitly describe the distributions of MRF. A Gibbs distribution is any distribution which can be expressed in the following form: ! X 1 1 PðXÞ ¼ expðUðXÞÞ ¼ exp  V c ðXÞ ; ð7Þ Z Z c2C where Z is a normalizing constant known as the partition function, U(X) is an energy function, and Vc(X) is the potential function contributed from the sites in the clique c. In fact, the Hammersley–Clifford theorem [19,21,22] states that if X has a strictly positive density function, then X is a MRF if and only if its distribution is in the form of a Gibbs distribution. Using the MRF model (7), the MAP estimate of a HR image in (5) can be written as follows: ( ) X b X ¼ arg minfUðXÞg ¼ arg min V c ðXÞ : ð8Þ X

X

c2C

Each potential function Vc(X) characterizes the interactions among a group of local pixels. With a different choice of potential function, a variety of distinct models exist within the class of MRF. In this paper, we consider pairs of symmetrical cliques within the 3  3 neighborhood of each pixel, and a quadratic cost function is used. We employ the particular second-order neighborhood model in [18], with the Gibbs energy function given by X

V c ðX; hÞ ¼

c2C

q:N 1 1 q:N 2 1 X X k¼0

fb1 ½ðxk;l  xk;lþ1 Þ2 þ ðxk;l  xk;l1 Þ2 

l¼0

þ b2 ½ðxk;l  xk1;l Þ2 þ ðxk;l  xkþ1;l Þ2  þ b3 ½ðxk;l  xk1;lþ1 Þ2 þ ðxk;l  xkþ1;l1 Þ2  þ b4 ½ðxk;l  xk1;l1 Þ2 þ ðxk;l  xkþ1;lþ1 Þ2 g;

ð9Þ

where qN1  qN2 represents the size of the HR image, and the parameter set h = [b1, b2, b3, b4] is called the clique parameters. Fig. 1 shows the second-order neighborhood of the pixel xk,l. It can be seen that the four terms in (9) measure the edge intensities in the vertical, horizontal, diagonal and anti-diagonal directions,

⋅ ⋅ ⋅

⋅⋅⋅

xk−1,l−1

xk−1,l

xk−1,l+1

xk,l−1

xk,l

xk,l+1

xk+1,l−1

y

x

xk+1,l

xk+1,l+1

⋅ ⋅ ⋅ Fig. 1. The neighborhood of pixel xk,l.

respectively, in an image. An MRF prior model is specified in terms of a set of the parameters. 3.2. Adaptive MRF parameter estimation The MRF model can be used to improve the quality of a reconstructed image, but errors may occur if the estimated parameters of the model do not accurately characterize the image. Therefore, it is desirable to estimate the corresponding parameters that allow the prior model to be adjusted in order to achieve an elegant solution for a given data set. During past years, many methods have been proposed to estimate the MRF parameters, such as the least squares (LS) [19,22], maximum likelihood (ML) [19,23], and Markov Chain Monte Carlo (MCMC) [24,25] methods. However, most of the previous parameter estimation methods are not appropriate for image magnification in the MAP framework for three reasons. First, these MRF parameter-estimation methods are computationally expensive. Second, they are usually used to estimate the parameters of the MRF for statistically homogeneous images. As natural-scene images are usually inhomogeneous, a single set of fixed parameters is not sufficient to characterize a whole image effectively. Finally, the HR image is not available for estimation of the parameters, and the parameters cannot be estimated from the LR image as the field property is not preserved across the scale or the resolution pyramid [18,19]. In this paper, we propose an effective and computationally efficient scheme for estimating the MRF model parameters used for image magnification. As a natural image is usually not homogeneous, and its edges are approximately straight in a small region, we therefore estimate the parameter set in a local manner using a blockwise adaptive parameter-estimation method. In our algorithm, we first partition an estimated HR image into sub-blocks of size p  p, then we estimate the respective parameter values of each sub-block. As these sub-blocks are small in size, we may assume that each sub-block is homogeneous and has a fixed set of parameters. The Gibbs energy function (9) contains four energy terms, each of which is associated with a weighting factor bi. These four terms are used to measure the vertical, horizontal, diagonal and anti-diagonal edges, respectively, in a sub-block. Therefore, the parameter set h can be used to control the degree of smoothness for the different directions in the image block. For example, if the vertical edge intensity in the block is very large, the corresponding values of b1 should be small so that the vertical edge will not be penalized in the minimization process. This will allow the resulting estimated image block to contain the vertical edge. Similarly, when the block contains strong diagonal or anti-diagonal edge intensity, the corresponding parameters b3 and b4 should be small. In this setting, the staircase artifact can be reduced when diagonal patterns appear in the image. In other words, the parameter of a certain direction should be set as inversely proportional to its corresponding energy value. The parameter estimation is performed as described below. (1) Compute the energy function in each of the four edge directions for each sub-block as follows:

⋅⋅⋅

u1 ¼

mþp1 X X nþp1

½ðxk;l  xk;lþ1 Þ2 þ ðxk;l  xk;l1 Þ2 ;

k¼m

u2 ¼

l¼n

mþp1 X X nþp1

½ðxk;l  xk1;l Þ2 þ ðxk;l  xkþ1;l Þ2 ;

k¼m

u3 ¼ 0:5 

l¼n mþp1 X X nþp1 k¼m

u4 ¼ 0:5 

l¼n

mþp1 X X nþp1 k¼m

ð10Þ ½ðxk;l  xk1;lþ1 Þ2 þ ðxk;l  xkþ1;l1 Þ2 ; and

l¼n

½ðxk;l  xk1;l1 Þ2 þ ðxk;l  xkþ1;lþ1 Þ2 ;

1280

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

where m and n represent the coordinates of the top-left location of a sub-block in an image, p  p is the size of the subblocks, and u1, u2, u3, and u4 are the energy values for the horizontal, vertical, diagonal, and anti-diagonal directions, respectively. As the distance between a pair of pixels in the diagonal direction is longer than that in the vertical or the horizontal direction, an appropriate weighting factor is included in the four energy terms in (10). We assign an appropriate weighting factor for each of the energy terms according to the square of the Euclidean distance between the related pixels. Thus, the weight is set at 1 for the horizontal and vertical directions, and at 0.5 for the diagonal and anti-diagonal directions. (2) As analyzed above, the parameter for each sub-block is set as inversely proportional to its energy value ui. For the sake of simplicity, we define a parameter u0i as follows: u0i

ð11Þ

(3) The estimated parameters u0i are normalized by dividing them by their sum for each block, so that all the parameters in the whole image are consistent. 0 j¼1 uj

g k ¼ W k ðg 0 Þ:

:

ð12Þ

Using this method, we can calculate the corresponding parameters for each small block. The above MRF model with estimated parameters ensures that the edges in the four particular directions can be preserved. The most important characteristics of this estimation method are that: (1) it can reflect the smoothness and sharpness in the four directions in local regions of the estimated HR image, and (2) we can compute the model parameters on the fly with a low computational cost. In our MAP framework, an initial magnified image is needed, which can be produced using various methods. To achieve a highquality HR image, our algorithm adopts a fractal-based method, namely MPIFS [13], which can provide better quality than that produced using the conventional interpolation methods. Then, the estimated parameters for the MRF are used to further regularize the initial HR image. MPIFS will be briefly described in the next section.

ð13Þ

Thus, K magnified images, g1,g2, . . . , gK, can be generated. Finally, the magnified image g is acquired by averaging them as follows: g¼

1 ¼ : ui

u0 bi ¼ P4 i

ceed at any image size, regardless of the size of the image that was used for encoding, a magnified image can be acquired when the PIFS is applied iteratively to an arbitrary larger size image. In our previous proposed method [13], different range-block partitioning schemes are applied to the image f. The partitioning schemes are denoted as fRnðkÞ g16n6NðkÞ ; 1 6 k 6 K, where K is the total r number of partitioning schemes used, and N ðkÞ r is the corresponding total number of range blocks for the kth partitioning scheme. The PIFS Wk of the image f is computed based on the kth partitioning scheme, which can then be applied to an arbitrary image g0 of size qN1  qN2 to generate a corresponding magnified image gk, i.e.

K 1X g : K k¼1 k

ð14Þ

By using this multiple PIFS (MPIFS) method, the blocky artifact in the magnified image is reduced significantly. Although averaging is performed in (14), this will not blur the edges in g. This is because the averaging is performed over the pixels at the same pixel position in the respective K magnified images, rather than over the neighborhood of the pixels concerned. Fig. 2 shows the magnified images generated by three different methods: bilinear interpolation, conventional fractal magnification, and MPIFS magnification. To compare the visual quality of the respective methods, only a part of the magnified images is shown. It can be seen that the MPIFS method can remove the noise and preserve the edge very well when compared to the other two methods. Fractal image coding involves a lot of computations; most of the encoding time is spent on finding the best-matched domain block from a large domain pool to represent an input range block. Many fast algorithms, such as the canonical classification [11] and the kd tree search [26], have been proposed. In our algorithm, we first apply the canonical classification to the domain pool as in [11], and then use the k-d tree to search for the best-matched domain block in each class. For an image of size 256  256, the computa-

4. Image magnification using multiple pifs (MPIFS) Fractal image coding, which uses the iterated function system (IFS) and collage theorem, was first proposed by Barnsley [10]. One of the unique features of fractal image coding is that its decoding is resolution independent. Therefore, if an image is exactly selfsimilar, it can be magnified to any size without any distortion. Unfortunately, natural images are usually not exactly self-similar. Hence, the image generated using the fractal codes of an original image will not be exactly identical, but an approximation of the original image only. A fully automated, fractal-based, image coding technique is based on the partitioned iterative function system (PIFS). First, the image f is partitioned into non-overlapping range blocks fRn g16n6Nr of size nr  nr, and Nr is the total number of range blocks in the image. Let D ¼ fDn g16n6Nd be the collection of all possible domain blocks of size nd  nd within the image support, and Nd be the total number of domain blocks. For each range block Rn, the corresponding best matched domain block Dn with transformation wn will be found. wn is called the fractal code for the range block Rn, and the set W ¼ [ wn is the PIFS of the given image. According to the contraction n mapping theorem, the image can be decoded by iteratively applying the PIFS W to an arbitrary starting image. As the decoding can pro-

Fig. 2. The image ‘‘Lena” of size 128  128, magnified with a zoom factor of 4.

1281

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

tion of the IFS takes only one or two seconds, with a slight degradation in quality. However, with the averaging process in our algorithm, the resulting noise in the magnified images can be reduced significantly. 5. Image magnification algorithm In our algorithm, image magnification in a Bayesian framework is formulated as an optimization problem whose solution is obtained by minimizing a cost function. During the iteration, the magnified image should always satisfy the constraint from the LR image. In the LR image acquisition model (1), as one pixel in the LR image of size N1  N2 corresponds to a zooming block (of size q  q pixels) in the HR image of size qN1  qN2, the average value of the zooming block must be in constraint and equal to the value of the corresponding pixel value in the LR image. In other words, when a pixel in the LR image is magnified to q  q pixels, the pixels in the magnified HR image should also be subject to this constraint. However, with the use of this LR image constraint, as in [5], blocky artifacts at the boundaries of the zooming blocks will appear in the estimated HR images. In our Bayesian MAP estimator, we also consider and incorporate another regularization constraint in order to achieve a solution with high visual quality. Therefore, we add a smoothness term S(X) into the cost function (15) to reduce the blocky artifacts. Substituting (9) into (8), the final cost function to be minimized is given as follows: X EðXÞ ¼ V c ðX; hÞ þ kSðXÞ c2C

¼

qN 1 =p qN 2 =p ðmþ1Þp1 X X X ðnþ1Þp1 X n

ðm;nÞ b1 ½ðxk;l

m¼0

n¼0

k¼mp

2

2

 xk;lþ1 Þ þ ðxk;l  xk;l1 Þ 

l¼np

þ b2 ½ðxk;l  xk1;l Þ2 þ ðxk;l  xkþ1;l Þ2  þ b3 ½ðxk;l  xk1;lþ1 Þ2  o ðm;nÞ þ xk;l  xkþ1;l1 Þ2  þ b4 ½ðxk;l  xk1;l1 Þ2 þ ðxk;l  xkþ1;lþ1 Þ2  ( "ðnþ1Þq1 #) ðmþ1Þq1 N1 X N2 X X X 2 2 ; ðxmq;l  xmq1;l Þ þ ðxk;nq  xk;nq1 Þ þk ðm;nÞ

m¼0 n¼0

ðm;nÞ

l¼nq

k¼mq

ð15Þ where k is a regularization parameter, p  p is the size of a sub-block, q  q is the size of a zooming block, and (m, n) is the index to the subblocks in the first term and to the zooming blocks in the second term. Note that the size of the zooming blocks may be different to that of the sub-blocks, but in our algorithms they are the same. ðm;nÞ ðm;nÞ ðm;nÞ ðm;nÞ hðm;nÞ ¼ ½b1 ; b2 ; b3 ; b4  denotes the MRF parameters of the sub-block (m, n), which is different for different sub-blocks. The above cost function consists of two terms: the first is based on the MRF prior model, while the second term S(X) measures the difference between the pixels at the borders of adjacent blocks. Regularizing an image may be seen as the minimization of the function E(X). As the cost function (15) is convex, the gradient-descent algorithm can be employed for minimization, i.e. X ðnþ1Þ ¼ X ðnÞ  aðnÞ rEðX ðnÞ Þ; (n)

ð16Þ (n)

where X is the estimated image in the nth iteration, a is the step size which can be fixed or updated adaptively during the iteration, and $E(X(n)) is the gradient. In our algorithm, we first use the MPIFS scheme to form the initial HR image, which can provide more edge information. However, some texture information for the initial HR image is lost during the averaging process, so we regularize the initial HR image by using the MRF prior with an adaptive scheme for parameter estimation, as well as with the LR image constraint and the border smoothness constraint. The optimization procedure using the gradient-descent algorithm is summarized in the following steps:

(1) Compute the initial HR image X(0) using the MPIFS algorithm, as described in Section 4; (2) Initialize all the parameters of the MRF model to 1; (3) Estimate the model parameters h(n) based on the estimated HR image X(n) for each sub-block, as defined in (10)–(12); (4) Compute the gradient $E(X(n)) of the cost function (15) with the parameters h(n); (5) Project the gradient into the LR image constraint space using the projection matrix P defined in [5]: ðnÞ

d

¼ PðrEðX ðnÞ ÞÞ:

ð17Þ

(6) Compute the step size using the Hessian matrix H of the cost function E(X(n)), i.e. ðnÞ

aðnÞ ¼

ðnÞ

½d T ½d  ðnÞ T

ðnÞ

½d  H½d 

:

ð18Þ

(7) Update the image X(n+1) as (16); and (8) If kX ðnþ1Þ  X ðnÞ k kX ðnÞ k

6 e;

ð19Þ

b ¼ X ðnþ1Þ , otherwise, return to Step 3 for the next iteration. then X After a number of iterations, the regularized estimated HR image will be obtained. 6. Experiments In this section, we demonstrate the efficacy of our proposed algorithm based on the MPIFS and the adaptive MRF model. To measure the quality of a magnified image, we down-sample a number of standard images of size 512  512 to 128  128 to form the LR images, using a 4  4 moving-average prefilter, as shown in Fig. 3. In other words, we will use the LR images and magnify them by a zooming factor of 4 to produce the corresponding HR images. b are then compared to the original HR These magnified images X images X. The performance of a magnification algorithm is measured by the PSNR, which is defined as follows: PSNR ¼ 10log10

2552 ; MSE PM PN 1

ð20Þ

b ij Þ2 , and M and N are the image where MSE ¼ MN i¼1 j¼1 ðX ij  X width and height, respectively. In the MPIFS magnification process, the sizes of the range blocks and the domain blocks are set to nr = 4 and nd = 8, respectively. Therefore, 16 partition schemes are available, and the respective starting points of the sub-images for MPIFS are (i, j), where 0 6 i, j < nr. These 16 LR sub-images are magnified and then averaged to form the initial estimated HR image. As fractal coding involved in this part is computationally intensive, the k-d tree is applied to speed up the search with a slight loss in quality. In the MAP optimization process, the regularization parameter k is set to 0.1, which is an optimal value determined by experiment. The HR image is partitioned into blocks of size 4  4. With this small block size, we can estimate the parameter sets in a local manner with the assumptions that each block is homogeneous and its edges are straight. If a larger block size is used, the above assumptions may no longer be valid. However, for a smaller block size, the local energy is sensitive to noise, and so is the local MRF parameter’s estimation. Our experiments have shown that the use of a block size of 4  4 in our algorithm can achieve the best overall performance. The MRF model parameters are then estimated for each block. Since the MPIFS can provide an accurate estimate of the HR images, we found that the condition in (19) will be satisfied with no more than 10 iterations by experiment. Therefore, in our algorithm, the maximum number of

1282

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

Fig. 3. The testing images.

iterations is set to 10. Table 1 shows the quality of the magnified images in terms of PSNR with different methods, where FMRF denotes the fixed-parameter MRF model, while AMRF denotes the adaptive parameter MRF model proposed in this paper. NN, BL and MPIFS mean that the initial HR image is produced by the nearestneighbor interpolation, by the bilinear interpolation, and by the MPIFS scheme, respectively. Therefore, MPIFS-AMRF is the method that uses an initial HR image generated by MPIFS, and then employs the AMRF model to improve the image quality. The performances of the nearest-neighbor (NN) interpolation and the MPIFS methods are also shown. Figs. 4 and 5 illustrate the respective magnified images for the images ‘‘Lena” and ‘‘Boat” based on the different methods, where we show only a portion of the magnified images. Table 1 shows that MPIFS-AMRF outperforms all other methods in terms of PSNR. When compared to MPIFS-FMRF, our adaptive MRF model can have a 0.34 dB gain on average. The performance of our algorithm is particularly good for those images with sharp edges in diagonal directions, such as the images Lena (0.5 dB), Boat (0.3 dB), Peppers (0.56 dB) and Woman (0.41 dB). From the magnified images (d), (e), and (f) in Figs. 4 and 5, the images generated by MPIFS are visually sharper, but some texture information is smoothed during the averaging process. Serious staircase artifacts appear at the diagonal boundaries if an image is regularized by the fixed-parameter MRF model. Using our proposed method, the staircase effect is reduced significantly, and a much better visual quality can be achieved as sharp edges and texture information can also be preserved very well. This can be observed from the rim of the hat in Fig. 4(f), and the masts in Fig. 5(f).

It should be noted that instead of the MPIFS, the NN interpolation method (or other simple interpolation methods) can also be used to generate the initial HR image. Our adaptive MRF model can then be applied to improve its quality. In this case, the computational complexity can be reduced. However, the quality of the initial image generated using NN is lower than that of MPIFS, so it will need a longer time, i.e. more iterations, to converge to the optimal solution. Thus, the maximum number of iterations is set to 20 in order to achieve a good result. The regularization parameter k is set to 0.5, as the blocky artifact is more serious when the nearest-neighbor interpolation is used (Figs. 4(a) and 5(a)). Table 1 also compares the PSNR of the different images produced by NN with the FMRF and the AMRF method, and BL with AMRF method. The performance of the MAP method proposed by Schultz et al. [5], with the initial image generated by NN interpolation, is also shown in Table 1. In this method [5], the energy function is slightly different to that of (9), and it uses a fixed set of parameters. The discontinuities are preserved by using the Huber function. However, the method requires knowledge of the edge magnitude for the selection of the threshold value for the Huber function. In addition, a single threshold value cannot accurately describe edges in a real image, so this reduces its efficacy. The results presented in Table 1 are based on a Huber threshold value of 10, which was selected by experimentation. We can also observe that its performance is better than that of the NN with the fixed parameter MRF model (9), i.e. NN-FMRF. However, the NN-AMRF method outperforms [5], with a PSNR improvement of 0.14 dB on average. From Table 1, we can also see that our proposed method performs

Table 1 The performance of the magnified images in terms of PSNR (dB) based on different methods with the zooming factor = 4 Images

NN

Schultz [5]

NN-FMRF

NN-AMRF

BL-AMRF

MPIFS

MPIFS-FMRF

MPIFS- AMRF

Lena_128 Boat_128 Peppers_128 Crowd_128 Goldhill_128 Woman_128 Man_128 Airplane_128 Average

26.92 25.07 25.99 24.65 26.59 32.02 25.58 24.69 26.44

29.46 26.77 28.43 27.27 27.95 36.62 27.32 26.74 28.82

29.37 26.71 28.28 27.16 27.92 36.48 27.23 26.68 28.73

29.79 26.94 28.55 27.47 27.94 36.72 27.52 26.76 28.96

29.76 26.86 28.55 27.41 27.91 36.67 27.43 26.89 28.94

27.76 25.23 27.54 25.16 26.53 34.12 25.89 25.08 27.16

29.40 26.73 28.35 27.19 27.92 36.54 27.25 26.73 28.76

29.90 27.03 28.91 27.48 27.98 36.95 27.57 26.96 29.10

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

1283

Fig. 4. The magnified images of size 512  512, generated from the LR image ‘‘Lena” of size 128  128.

Fig. 5. The magnified images of size 512  512, generated from the LR image ‘‘Boat” of size 128  128.

especially well for those images with sharp edges. Most importantly, our adaptive MRF approach can always produce HR images with a superior visual quality, as illustrated in (c) and (f) in Figs. 4 and 5. The visual quality improves significantly at the edge regions. Nevertheless, the hair-like structure of the hat in Fig. 2(d) is over-smoothened. Our MRF model cannot help to enhance this region. This is because during the MRF parameter estimation, the parameters are estimated based on the energy of four different

directions. Our algorithm can therefore achieve a better visual quality for the edges in all the four directions. However, within a textural area such as the hair-like structure of the hat, the parameter estimation is no longer effective, so the visual quality of this kind of regions is less satisfactory. In general, the computational complexity of our method, based on MPIFS and AMRF, is higher than that of other methods. The experiments were conducted on a Pentium IV 2.4 GHZ computer

1284

X. Zhang et al. / Image and Vision Computing 26 (2008) 1277–1284

Experiments have shown that the proposed method outperforms the method based on the MRF model with a fixed set of parameters in terms of both PSNR and visual quality. The PSNR of the method MPIS-AMRF is 0.06–0.56 dB higher than when the MPIFS-FMRF is used, while the NN-AMRF is 0.02–0.42 dB higher than when the NN-FMRF is used. Acknowledgement The work described in this paper was supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. PolyU 5199/06E). References

Fig. 6. Comparison of visual quality with the use of (a) a single set of fixed parameters and (b) the adaptive MRF model.

with 512 MB memory. On average, the runtime for generating an initial HR image using MPIFS from a LR image of size 128  128 is 9.4 s. The average runtimes required for the AMRF and FMRF models are 11.2 s and 5.8 s, respectively. The runtime required for AMRF is longer than that for FMRF because the MRF parameters of each sub-block are computed at each iteration. Although the NN method, which is computationally simple, can also be used to generate the initial HR image, it will require more iterations in the AMRF process and the quality of the resulting image is lower than that attained by using MPIFS-AMRF. Based on our experiments, the average runtimes for MPIFS-AMRF and NN-AMRF are 20.6 s and 22.4 s, respectively. When the zooming factor is 8, our method can also produce a satisfactory visual quality. Fig. 6 compares two magnified images with a zooming factor of 8 using the fixed and our adaptive MRF model together with MPIFS. 7. Discussion and conclusion In this paper, we have proposed an image magnification method based on MPIFS and an adaptive MRF prior model in the MAP framework. The initial HR image is generated by the MPIFS method, which can maintain the edge information effectively. Then the MRF parameters are estimated from the initial HR image, and the specified MRF model is used to further regularize the estimated HR image, which can simultaneously preserve the edge with the adaptive estimated parameters. At each iteration, the model parameters will be updated depending on the estimated HR image of the last iteration. In addition, using a simple interpolation method to generate the initial image, our proposed adaptive MRF model can still provide a magnified image of high visual quality.

[1] W.K. Pratt, Digital Image Processing, John Wiley & sons Inc., New York, 2001. [2] M. Unser, Splines: a perfect fit for signal and image processing, IEEE Signal Processing Magazine 16 (6) (1999) 22–38. [3] L. Zhang, X. Wu, An edge guided image interpolation algorithm via directional filtering and data fusion, IEEE Transactions on Image Processing 15 (8) (2006) 2226–2238. [4] X. Li, M.T. Orchard, New edge-directed interpolation, IEEE Transactions on Image Processing 10 (10) (2001) 1521–1527. [5] R.R. Schultz, R.L. Stevenson, A Bayesian approach to image expansion for improved definition, IEEE Transactions on Image Processing 3 (3) (1994) 233– 242. [6] W.K. Carey, D.B. Chuang, S.S. Hemami, Regularity preserving image Interpolation, IEEE Transactions on Image Processing 8 (9) (1999) 1293–1297. [7] S. Chaudhuri, Super-Resolution Imaging, Kluwer Academic Publishers, Boston, 2001. [8] S.K. Mitra, C.A. Murthy, M.K. Kundu, A technique for image magnification using partitioned iterative function system, Pattern Recognition 33 (7) (2000) 1119– 1133. [9] K.H. Chung, Y.H. Fung, Y.H. Chan, Image enlargement using fractal, in: Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, 2003, pp. 273–276. [10] M. Barsley, Fractal Everywhere, Academic Press, New York, 1988. [11] Y. Fisher, Fractal Image Compression: Theory and Application, Springer-Verlag, Berlin, Germany, 1995. [12] E. Reusens, Overlapped adaptive partitioning for image coding based on the theory of iterated functions systems, in: Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, 1994, pp. 357–364. [13] C.M. Lai, K.M. Lam, W.C. Siu, An efficient fractal-based algorithm for image magnification, proceedings, in: Proceedings of the International Symposium on Intelligent Multimedia, Video and Speech Processing, October 2004, pp. 571–574. [14] A. Tsai, A. Yezzi Jr., Curve evolution implementation of the Mumford–Shah functional for image segmentation, denoising, interpolation, and magnification, IEEE Transactions on Image Processing 10 (8) (2001) 1169– 1186. [15] D. Tschumperle, R. Deriche, Vector-valued image regularization with PDEs: a common framework for different applications, IEEE Transactions on Pattern Analysis and Machine Intelligence 27 (4) (2005) 506–517. [16] S. Borman, R.L. Stevenson, Super-resolution from image sequences – a review, in: Proceedings of the 1998 Midwest Symposium on Circuits and Systems, vol. 5, April 1998. [17] S.C. Park, M.K. Park, M.G. Kang, Super-resolution image reconstruction: a technical overview, IEEE Signal Processing Magazine 20 (3) (2003) 21–36. [18] M.V. Joshi, S. Chaudhuri, R. Panuganti, A learning-based method for image super-resolution from zoomed observation, IEEE Transactions on System Man and Cybernatics B 35 (3) (2005) 527–537. [19] S.Z. Li, MRF Modeling in Computer Vision, Springer-Verlag, 1995. [20] Z. Kato, Modélisations markoviennes multirésolutions en vision par ordinateur application à la segmentation d’images SPOT, Ph.D. Thesis, 1994. [21] S. Geman, D. Geman, Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images, IEEE Transactions on Pattern Analysis Machine Intelligence 6 (1994) 721–741. [22] J.E. Besag, Spatial interaction and the statistical analysis of lattice systems, Journal of the Royal Statistical Society B 36 (2) (1974) 192–236. [23] S.S. Saquib, C.A. Bouman, K. Sauer, ML parameter estimation for Markov random fields with applications to Bayesian tomography, IEEE Transactions on Image Processing 7 (7) (1998) 1029–1044. [24] L. Wang, J. Liu, S. Li, MRF parameter estimation by MCMC method, Pattern Recognition 33 (11) (2000) 1919–1925. [25] Y. Yu, Q. Cheng, MRF parameter estimation by an accelerated method, Pattern Recognition Letters 24 (9–10) (2003) 1251–1259. [26] D. Saupe, Accelerating fractal image compression by multi-dimensional nearest neighbor search, in: Proceedings of the Data Compression Conference, March 1995.