The Nonlocal 𝑝-Laplacian Evolution for Image Interpolation

Abstract

This paper presents an image interpolation model with nonlocal 𝑝-Laplacian regularization. The nonlocal 𝑝-Laplacian regularization overcomes the drawback of the partial differential equation (PDE) proposed by Belahmidi and Guichard (2004) that image density diffuses in the directions pointed by local gradient. The grey values of images diffuse along image feature direction not gradient direction under the control of the proposed model, that is, minimal smoothing in the directions across the image features and maximal smoothing in the directions along the image features. The total regularizer combines the advantages of nonlocal 𝑝-Laplacian regularization and total variation (TV) regularization (preserving discontinuities and 1D image structures). The derived model efficiently reconstructs the real image, leading to a natural interpolation, with reduced blurring and staircase artifacts. We present experimental results that prove the potential and efficacy of the method.

1. Introduction

Digital image interpolation is an important technology in digital photography, TV, multimedia, advertising, and printing industries, which is applied to obtain higher-resolution image with better perceptual quality. The key task in image interpolation is to remove zigzagging, blurry, and other artifacts producing visually pleasing resulting image. Many literatures are devoted to address these problems [1–10]. These methods, usually known as edge directed, level set based, or isophote oriented, vary considerably.

The methods based on edge direction were proposed to obtain smooth edges of the resulting images [1–3]. However, these methods suffer from degradations on edge-free regions because they rely on local directions estimation, creating false edges in uniform regions. Wang and Ward have developed an interesting technique based on the detection of ridges (straight edges) in images [4], which allows them to interpolate directionally only pixels situated on straight edges and avoid the apparition of false edges. In [5], the zigzagging artifacts are reduced by restricting curvature of the interpolated isophotes (equi-intensity contours). Minimum curvature is required on isophotes of the interpolated images. Cha and Kim proposed a method based on the TV energy in order to remove the so-called checkerboard effect and to form reliable edges [6]. Morse and Schwartzwald [7] presented a scheme that uses existing interpolation techniques as an initial approximation and then iteratively reconstructs the isophotes using constrained smoothing. However, they are complex compared to traditional methods and thus computationally expensive. Malgouyres and Guichard [8] proposed to choose as solution of the interpolation the image that minimizes the TV. This method leads to resulting images without blurring effects, as it allows discontinuities and preserves 1D fine structures. Aly and Dubois [9] proposed a model-based TV regularization image up sampling methods. Image acquisition process is modeled after a lowpass filtering followed by sampling. However, TV minimization is based on the assumption that the desirable image is almost piecewise constant, which yields a result with over smoothed homogeneous regions. Belahmidi and Guichard [10] have improved the TV-based interpolation by developing a nonlinear anisotropic PDE, hereafter referred to as BG interpolation method. In order to enhance edge preservation, this PDE performs a diffusion with strength and orientation adapted to image structures. This method balances linear zooming on homogeneous regions and anisotropic diffusion near edges, trying to combine the advantages of these two processes. The anisotropic diffusion scheme, including edge-directed or isophote-oriented method, uses the gradient to extract the image feature (edge) direction, that is, the gradient direction is considered to be the direction across the image feature. Nevertheless, the information contained in the gradient is local, not good to determine the edge directions. Nonlocal information should be considered in determining edge directions.

Recently, a nonlocal evolution equation and variations of it have been widely used to model diffusion processes in many areas [11, 12]. Let us briefly introduce some references of nonlocal problem considered along this work. A nonlocal evolution equation corresponding to the Laplacian equation is presented as follows𝑢_𝑡(𝑥,𝑡)=𝐽∗𝑢−𝑢(𝑥,𝑡)=_𝑅𝑁𝐽𝑥−𝑦𝑢𝑦,𝑡−𝑢(𝑥,𝑡)𝑑𝑦. (1.1) It is called a nonlocal diffusion equation since the diffusion of the density at a point 𝑥 and time 𝑡 does not only depend on 𝑢(𝑥;𝑡), but on all the values of 𝑢 in a neighborhood of 𝑥 through the convolution term 𝐽∗𝑢. This equation shares many properties with the classical heat equation 𝑢_𝑡=Δ𝑢. More precisely, as stated in [13], if 𝑢(𝑡,𝑥) is thought of as the density at the point 𝑥 at time 𝑡, and 𝐽(𝑥,𝑦) is thought of as the probability distribution of jumping from location 𝑦 to location 𝑥, then the convolution (𝐽∗𝑢)(𝑥,𝑡)=∫_𝑅𝑁𝐽(𝑦−𝑥)𝑢(𝑦,𝑡)𝑑𝑦 is the rate at which individuals are arriving to position 𝑥 from all other places, and −𝑢(𝑥,𝑡)=−∫_𝑅𝑁𝐽(𝑦−𝑥)𝑢(𝑦,𝑡)𝑑𝑦 is the rate at which they are leaving location 𝑥 to travel to all other sites. This nonlocal evolution can be thought of as nonlocal isotropic diffusion.

In this paper, we propose a new method for image interpolation based on nonlocal 𝑝-Laplacian evolution. The nonlocal 𝑝-Laplacian and TV act as a regularizer to restrict edges of resulting image. The evolution is similar to an anisotropic energy dissipation process. The diffusion performs accurately along the direction of edges curves and its orthogonal direction. The magnitude of |𝑢(𝑦,𝑡)−𝑢(𝑥,𝑡)|^𝑝−2 determines diffusion directions. It suppresses diffusion across the image feature direction and enhances diffusion along the image feature direction.

The rest of the paper is organized as follows. In Section 2, we review the method based on image up sampling with TV regularization in [9] and BG image interpolation [10]. The proposed nonlocal 𝑝-Laplacian evolution for image interpolation is presented in Section 3. In Section 4, we demonstrate the experimental results to verify the effectiveness of our method, and the last section is for conclusion.

2. Background

We can think that a digital lower-resolution image 𝑢₀ (input image) defined on some lattice is obtained by transforming a high-resolution image 𝑢 (output image) defined on some better precision lattice, that is,𝑢₀=𝐻𝑢, (2.1) where 𝐻 is a sparse matrix that combines both filtering and down sampling process. The goal of image interpolation is to solve the inverse problem (2.1), an ill-posed inverse problem. This ill-posed inverse problem is generally approached in a regularization-based framework, which would be formulated as an energy functional [9],𝐸(𝑢)=𝐽_𝑑𝑢,𝑢₀+𝜆𝐽_𝑟(𝑢), (2.2) where 𝜆 is a regularization parameter that controls the tradeoff between 𝐽_𝑑 and 𝐽_𝑟. The data fidelity function 𝐽_𝑑 generally is formulated in the classical least-squares sense as 𝐽_𝑑(𝑢,𝑢₀)=(1/2)|𝐻𝑢−𝑢₀|². The TV regularizer 𝐽_𝑟 is taken as 𝐽_𝑟(𝑢)=∫_Ω|∇𝑢(𝑥)|𝑑𝑥. The formula (2.2) is rewritten as follows𝐸(𝑢)=𝜆_Ω||1∇𝑢𝑑𝑥+2𝐻𝑢−𝑢₀². (2.3) Using Euler's equation, the minimizer is the steady-state solution of the nonlinear PDE given by𝑢_𝑡=𝜆|∇𝑢|div∇𝑢|∇𝑢|+𝐻^𝑇𝐻𝑢−𝑢₀. (2.4)

On the other hand, Belahmidi and Guichard solve the ill-posed inverse problem based on the classical heat diffusion model [10]. Let 𝜂 denote the direction of local gradient, and 𝜉 the direction perpendicular to the gradient, namely,1𝜂=||𝑢∇𝑢^𝑥,𝑢_𝑦1,𝜉=|∇𝑢|−𝑢_𝑦,𝑢_𝑥. (2.5)

The second-order directional derivatives of the image 𝑢 along the directions of 𝜂 and 𝜉 are easily computed as follows:𝜕²𝑢𝜕𝜉²𝑢=𝜉^𝑥𝑥𝑢_𝑥𝑦𝑢_𝑥𝑦𝑢_𝑦𝑦𝜉^𝑇=𝑢_𝑥𝑥𝑢²_𝑦−2𝑢_𝑥𝑦𝑢_𝑥𝑢_𝑦+𝑢_𝑦𝑦𝑢²_𝑥|∇𝑢|²,𝜕²𝑢𝜕𝜂²𝑢=𝜂^𝑥𝑥𝑢_𝑥𝑦𝑢_𝑥𝑦𝑢_𝑦𝑦𝜂^𝑇=𝑢_𝑥𝑥𝑢²_𝑥+2𝑢_𝑥𝑦𝑢_𝑥𝑢_𝑦+𝑢_𝑦𝑦𝑢²_𝑦|∇𝑢|². (2.6) The interpolation scheme based on heat diffusion model is formulated as follows:𝑢_𝑡=||𝜕∇𝑢²𝑢𝜕𝜉²𝜕+𝑔(|∇𝑢|)²𝑢𝜕𝜂²−𝐻^𝑇𝐻𝑢−𝑢₀. (2.7) In this equation, the function 𝑔(𝑠) is typically defined as𝑔(𝑠)=11+(𝑠/𝐾)², (2.8) with 𝐾>0 is a constant to be tuned for a particular application. The role of the diffusion coefficient 𝑔(|∇𝑢|) is to control the smoothing adaptively.

When 𝑔≡0, (2.7) reduces to (2.4) with 𝜆=1. All these models can be viewed as interpolation schemes based on nonlinear diffusion model. The two regularizers |∇𝑢|div(∇𝑢/|∇𝑢|) and |∇𝑢|(𝜕²𝑢/𝜕𝜉²)+𝑔(|∇𝑢|)(𝜕²𝑢/𝜕𝜂²), respectively, in (2.4) and (2.7) result in different interpolation effects. In fact, |∇𝑢|div(∇𝑢/|∇𝑢|)=|∇𝑢|(𝜕²𝑢/𝜕𝜉²) is the second-order directional derivative in the direction that is orthogonal to the gradient |∇𝑢|, and 𝜕²𝑢/𝜕𝜂² is the second-order directional derivative in the direction of the gradient |∇𝑢|.

From the viewpoint of geometry, the evolution processes in the artificial time 𝑡 given by these models are seen as energy dissipation processes in two orthogonal directions 𝜂 and 𝜉 [9]. The diffusion process of 𝑢(𝑥,𝑡) along 𝜉 will preserve the location and the intensity transitions of the contours, while smoothing along them maintaining their crispness. This diffusion term is used to maintain edges with smooth isophotes in [9, 10]. The diffusion of the grey values along 𝜂 walks across both sides of the local image contour. This process blurs crisp contours as in the case of linear interpolators. Two divided means are used to deal with it. The diffusion process along 𝜂 is cast aside in [9], while controlled by edge-stopping function 𝑔(|∇𝑢|) to balance the two diffusion terms in [10].

TV regularization in (2.3) does an excellent job at preserving edges while reconstructing images [14]. This phenomenon can also be explained physically, since the resulting diffusion is strictly orthogonal to the gradient of the image. But TV-based interpolation favors solutions that are piecewise constant. This sometimes causes a staircasing effect in homogeneous regions, which are long observed in the literature in denoising, for example, [15, 16]. Not only having blocky solutions, but they can also develop false edges in resulting image.

The function g in (2.7) is to be chosen with values between 0 and 1. The energy dissipation process (2.7) is adaptively controlled in the direction 𝜂 and 𝜉, that is, minimal smoothing in the directions 𝜂 (across the image features) preserving sharp edges, and maximal smoothing in the directions 𝜉 (along the image features) obtaining smooth contours. The anisotropic diffusion scheme that uses the gradient to extract the image feature direction can mistakenly give maximal smoothing to the across feature direction and severely damage the image features, especially the image lines and textures [17]. And blurry and/or oscillatory edges are introduced in interpolated image. The drawback of this model is that the gradient used to extract the image feature direction is too local. The information contained in the gradient is limited to a point and its immediate neighbors, while the edge curve that determines the edge directions is not a local event. The interpolation direction extraction should base on a larger neighborhood.

3. Nonlocal 𝑝-Laplacian Image Interpolation

In this section, we adopt nonlocal 𝑝-Laplacian evolution to overcome the local limit. Our proposed energy functional for regularized image interpolation is given by𝐸(𝑢)=𝛼_Ω||𝛽∇𝑢𝑑𝑥+2𝑝_Ω𝐽𝑥−𝑦𝑢𝑦−𝑢(𝑥)^𝑝1𝑑𝑦𝑑𝑥+2𝐻𝑢−𝑢₀². (3.1) The first part (the sum of the first term and the second) of right-hand side is regularizer, and the other is data fidelity function. The gradient flow associated to the functional 𝐸(𝑢) is𝑢_𝑡(𝑥,𝑡)=𝛼div|∇𝑢|⁻¹∇𝑢+𝛽𝑃^𝐽_𝑝(𝑢)−𝐻^𝑇𝐻𝑢−𝑢₀,𝑢(𝑥,0)=𝐻^𝑇𝑢₀, (3.2) where𝑃^𝐽_𝑝(𝑢)=_Ω𝐽𝑥,𝑦𝑢𝑦,𝑡−𝑢(𝑥,𝑡)^𝑝−2𝑢𝑦,𝑡−𝑢(𝑥,𝑡)𝑑𝑦. (3.3) The kernel 𝐽∶Ω→𝑅 is assumed to be nonnegative, bounded continuous radial function, with sup(𝐽)⊂𝐵(0,𝑑) and ∫_Ω𝐽(𝑧)𝑑𝑧=1.

The nonlocal energy dissipation is implemented mostly by 𝑃^𝐽_𝑝(𝑢) in our model. It is necessary to investigate the relation between the heat diffusion equation related to (2.7) and the 𝑝-Laplacian equation. The 𝑝-Laplacian evolution equation 𝑢_𝑡=div(|∇𝑢|^𝑝−2∇𝑢) is well studied in image processing, which can be represented as [18]𝑢_𝑡=|∇𝑢|^𝑝−2𝜕²𝑢𝜕𝜉²+𝑝−1|∇𝑢|^𝑝−2𝜕²𝑢𝜕𝜂². (3.4) When 𝑝=1, it is the TV flow keeping the edges but suffering from the staircase effect. When 𝑝=2, 𝑢_𝑡=Δ𝑢, this is isotropic diffusion because of the same diffusion coefficients. This model can smooth image, while bluring sharp edges. When 1<𝑝<2, the grey values of 𝑢(𝑥,𝑡) in (3.4) diffuse along the directions 𝜂 and 𝜉, respectively, as in the following equation:𝑢_𝑡=||𝜕∇𝑢²𝑢𝜕𝜉²𝜕+𝑔(|∇𝑢|)²𝑢𝜕𝜂², (3.5) which has different control factor compared with (3.4). That is to say, (3.4) and (3.5) all preserve advantage of adaptive smoothing.

A nonlocal improvement to 𝑝-Laplacian equation was studied in [19],𝑢_𝑡(𝑥,𝑡)=_Ω𝐽𝑥,𝑦𝑢𝑦,𝑡−𝑢(𝑥,𝑡)^𝑝−2𝑢𝑦,𝑡−𝑢(𝑥,𝑡)𝑑𝑦. (3.6) Moreover, nonlocal problems of type (3.6) have been used recently in the study of deblurring and denoising of images [20]. With Neumann boundary conditions, the solutions to problem (3.6) converge to the solution of the classical 𝑝-Laplacian if 𝑝>1 [19]. The nonlocal 𝑝-Laplacian evolution (3.6) improves the limit of the diffusion direction extraction depending on gradient (local information). The diffusion of the density at a point 𝑥 and time 𝑡 depends on all the values of 𝑢 in a larger neighborhood of 𝑥. More precisely, if 𝑢(𝑡,𝑥) is thought of as the density at the point 𝑥 at time 𝑡, and 𝐽(𝑥,𝑦) is thought of as the probability distribution of jumping from location 𝑦 to location 𝑥, then 𝑃^𝐽_𝑝(𝑢) is the rate at which individuals are arriving to position 𝑥 from all other places. In image interpolation, 𝑃^𝐽_𝑝(𝑢) is also the rate at which individuals are devoting to interpolated pixel 𝑥 from all other pixels. The evolution process in the artificial time 𝑡 given by (3.6) is seen as an anisotropic energy dissipation process. The direction of anisotropic diffusion is indicated by the magnitude of |𝑢(𝑦,𝑡)−𝑢(𝑥,𝑡)|^𝑝−2 in a larger neighborhood. It approximates to the direction of edge curve more accurate than the direction indicated by gradient. When 1<𝑝<2, the energy dissipation process is adaptively controlled by |𝑢(𝑦,𝑡)−𝑢(𝑥,𝑡)|^𝑝−2 along the direction of edge curves and the orthogonal direction to edge curves. The diffusion process along the direction of edge curves is suppressed for small |𝑢(𝑦,𝑡)−𝑢(𝑥,𝑡)|^𝑝−2, and the diffusion along the orthogonal direction is enhanced for larger |𝑢(𝑦,𝑡)−𝑢(𝑥,𝑡)|^𝑝−2. This results in minimal smoothing in the directions across the image features preserving sharp edges and maximal smoothing in the directions along the image features reducing zigzagging artifacts and oscillatory.

4. Numerical Algorithm and Experimental Results

In this section, we develop a fully discrete numerical method to approximate problem (3.2). We recall first the notations in the finite differences scheme used in our paper. Let ℎ and Δ𝑡 be the space and time steps, respectively, and let (𝑥_1𝑖;𝑥_2𝑗)=(𝑖ℎ;𝑗ℎ) be the grid points. Let 𝑢^𝑛(𝑖;𝑗) be an approximation of the function 𝑢(𝑛Δ𝑡;𝑥_1𝑖;𝑥_2𝑗), with 𝑛≥0. Equation (3.2) can be discretized as follows:𝑢^𝑛+1𝑖,𝑗−𝑢^𝑛𝑖,𝑗=Δ𝑡𝛼div∇𝑢|∇𝑢|−𝐻^𝑇𝐻𝑢−𝑢_0𝑖𝑗+𝛽_{()𝑘,𝑙∈Ω}𝐽(𝑘,𝑙)−𝑖,𝑗𝑢^𝑛(𝑘,𝑙)−𝑢^𝑛𝑖,𝑗^𝑝−2𝑢^𝑛(𝑘,𝑙)−𝑢^𝑛𝑖,𝑗. (4.1) In all numerical experiments, we choose the following kernel function:𝐽(𝑥)1∶=𝐶exp|𝑥|²−𝑑²if0𝑥<𝑑,if𝑥≥𝑑. (4.2) The constant 𝐶>0 is selected such that ∫_Ω𝐽(𝑥)𝑑𝑥=1.

We tested the proposed interpolation method on a variety of images. Some of the results are shown in Figures 1–3. Images are expanded by a factor of 2×2 in Figures 1 and 2 and by a factor of 10×10 in Figure 3. For comparison, we also show images interpolated using the BG image interpolation proposed in [10], the edge-guided image interpolation (EGI) in [21], and the edge-directed interpolation (NEDI) proposed in [3]. The choice of the parameters is based on subjective quality of the results assessed informally by our personal preference as human viewers in terms of edge sharpness, contour crispness, no ringing in smooth regions, and no ringing near edges. We use the following parameters: 𝛼=0.5, 𝛽=0.0001, and 𝑝=1.8 for the proposed interpolation method, 𝑘=0.0001 for the heat diffusion model, and time step Δ𝑡=0.15 for these experiments. There is a non visible improvement on subjective or objective quality of the results when the parameters are not badly changed. The iteration is terminated, when |𝑢^𝑛+1−𝑢^𝑛|²<10⁻⁶, normally within a few decade iterations.

(a)

(b)

(c)

(d)

Figure 1 

Barbara reduced and expanded by 2×2 (portion). (a) NEDI interpolation; (b) EGI interpolation; (c) BG method; (d) proposed method.

(a)

(b)

(c)

(d)

Figure 2 

Flower expanded by 2×2 (portion). (a) NEDI interpolation; (b) EGI interpolation; (c) BG method; (d) proposed method.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

A portion of Barbara, Mandrill, and house images expanded by 10×10. Top row: BG method; bottom row: proposed method.

In the first experiments, Barbara image with a size of 512×512 was lowpass filtered and subsampled by a factor of 2×2, then the subsampled image was interpolated to the original image size. The interpolation was performed by four different methods, and a portion of the results is shown in Figure 1. Figure 2 shows a portion of a result interpolating a given flower image with a size of 320×240 by a factor of 2×2 without subsample. From the two examples, the NEDI interpolation method tends to introduce the zigzagging artifacts (Figures 1(a) and 2(a)). The BG results produce slight blurry edges (stripes in Figure 1(c)) and ringing in smooth regions (as shown in Figure 2(c)), because the direction decided by gradient misses their real directions as stated in Section 2. Our proposed method and the EGI method produce sharp edges and smooth contours, but the EGI method and the NEDI interpolation are applied only by a factor of 2×2.

The second experiment directly interpolates a portion of images (Barbara, mandrill, and house) by a factor of 10×10 shown in Figure 3. This experiment is performed by the BG method and our proposed method since the other methods only resize image by a factor of 2×2. It is clear from the figures that results obtained with the proposed approach are better than the results by the BG method in [10]. The proposed method generates sharper and crisper stripes in Figure 3(d) compared to the result in Figure 3(a). The BG method produces blurry edges (the stripes in Figure 3(a), the beard in Figure 3(b)), zigzagging artifacts, and oscillatory (the line in Figure 3(c)), and images tend to be less natural. In the BG interpolation result, the boundaries of the text suffer from artifacts that make visualization difficult. It can be seen that our method results in an interpolated image with the fewer spurious patterns.

We use two measures, the classic PSNR and the mean structural similarity (MSSIM) index [22], to characterize the difference between the reference image and the output of a method. The MSSIM seems to approximate the perceived visual quality of an image better than PSNR or various other measures [23]. MSSIM index takes values in [0,1] and increases as the quality increases. It is calculated by the code available at http://www.cns.nyu.edu/lcv/ssim/, using the default parameters. We use several test images of size 512×512 including Lena, mandrill, and Barbara of size 256×256 including cameraman, barche, peppers, and resolution test. To show the true power of the interpolation algorithms, we first downsampled each image by a factor of 2×2 and then interpolated the result back to its original size. The PSNR is shown on Table 1 and the MSSIM on Table 2. From the two tables, the proposed method yields improved PSNR (except for the NEDI), and MSSIM results in all the experiments. This improvement may be attributed to the fact that the nonlocal 𝑝-Laplacian evolution works better than other methods.

Table 1 

Comparison of different interpolation algorithms using PSNR values.

Table 2 

Comparison of different interpolation algorithms using MSSIM values.

5. Conclusion

In this paper, a new image interpolation model based on TV and nonlocal 𝑝-Laplacian regularization is proposed. It combines the advantages of TV regularizer and nonlocal 𝑝-Laplacian regularizer, that is, allowing discontinuities and preserving 1D image structures and the diffusion of the grey values of images along image feature direction. It overcomes the drawbacks of the anisotropic diffusion proposed by Belahmidi and Guichard. The direction of anisotropic diffusion is indicated by the information of image feature in a larger neighborhood. This results in minimal smoothing in the directions across the image features preserving sharp edges and maximal smoothing in the directions along the image features reducing zigzagging artifacts and oscillatory. We have shown improvement over nonlocal 𝑝-Laplacian on a subjective scale, and in many cases with an improvement in PSNR and MSSIM. We expect to prove convergence of the evolution equation in future work.

Acknowledgments

This research was partially supported by Natural Science Foundation Project of CQ CSTC (Grant no. cstcjjA40012) and the National Natural Science Foundation of China (Grant no. 10871217).