The Nonlocal p-Laplacian Evolution for Image Interpolation

This paper presents an image interpolation model with nonlocal p-Laplacian regularization. The nonlocal p-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 p-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.


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 isophoteoriented 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 It is called a nonlocal diffusion equation since the diffusion of the density at a point x and time t does not only depend on u x; t , but on all the values of u in a neighborhood of x through the convolution term J * u.This equation shares many properties with the classical heat equation u t Δu.More precisely, as stated in 13 , if u t, x is thought of as the density at the point x at time t, and J x, y is thought of as the probability distribution of jumping from location y to location x, then the convolution J * u x, t R N J y − x u y, t dy is the rate at which individuals are arriving to position x from all other places, and −u x, t − R N J y − x u y, t dy is the rate at which they are leaving location x 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 p-Laplacian evolution.The nonlocal p-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 |u y, t − u x, t | p−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 p-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.

Background
We can think that a digital lower-resolution image u 0 input image defined on some lattice is obtained by transforming a high-resolution image u output image defined on some better precision lattice, that is, where H 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 , where λ is a regularization parameter that controls the tradeoff between J d and J r .The data fidelity function J d generally is formulated in the classical least-squares sense as The TV regularizer J r is taken as J r u Ω |∇u x |dx.The formula 2.2 is rewritten as follows Using Euler's equation, the minimizer is the steady-state solution of the nonlinear PDE given by 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,

2.5
The second-order directional derivatives of the image u along the directions of η and ξ are easily computed as follows:

2.6
The interpolation scheme based on heat diffusion model is formulated as follows: In this equation, the function g s is typically defined as with K > 0 is a constant to be tuned for a particular application.The role of the diffusion coefficient g |∇u| is to control the smoothing adaptively.When g ≡ 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 |∇u| div ∇u/|∇u| and |∇u| ∂ 2 u/∂ξ 2 g |∇u| ∂ 2 u/∂η 2 , respectively, in 2.4 and 2.7 result in different interpolation effects.In fact, |∇u| div ∇u/|∇u| |∇u| ∂ 2 u/∂ξ 2 is the secondorder directional derivative in the direction that is orthogonal to the gradient |∇u|, and ∂ 2 u/∂η 2 is the second-order directional derivative in the direction of the gradient |∇u|.
From the viewpoint of geometry, the evolution processes in the artificial time t given by these models are seen as energy dissipation processes in two orthogonal directions η and ξ 9 .The diffusion process of u x, t 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 g |∇u| 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.

Nonlocal p-Laplacian Image Interpolation
In this section, we adopt nonlocal p-Laplacian evolution to overcome the local limit.Our proposed energy functional for regularized image interpolation is given by

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 E u is 3.2 where P J p u Ω J x, y u y, t − u x, t p−2 u y, t − u x, t dy.

3.3
The kernel J : Ω → R is assumed to be nonnegative, bounded continuous radial function, with sup J ⊂ B 0, d and Ω J z dz 1.
The nonlocal energy dissipation is implemented mostly by P J p u in our model.It is necessary to investigate the relation between the heat diffusion equation related to 2.7 and the p-Laplacian equation.The p-Laplacian evolution equation u t div |∇u| p−2 ∇u is well studied in image processing, which can be represented as 18

3.4
When p 1, it is the TV flow keeping the edges but suffering from the staircase effect.When p 2, u t Δu, this is isotropic diffusion because of the same diffusion coefficients.This model can smooth image, while bluring sharp edges.When 1 < p < 2, the grey values of u x, t in 3.4 diffuse along the directions η and ξ, respectively, as in the following equation:  and J x, y is thought of as the probability distribution of jumping from location y to location x, then P J p u is the rate at which individuals are arriving to position x from all other places.In image interpolation, P J p u is also the rate at which individuals are devoting to interpolated pixel x from all other pixels.The evolution process in the artificial time t given by 3.6 is seen as an anisotropic energy dissipation process.The direction of anisotropic diffusion is indicated by the magnitude of |u y, t − u x, t | p−2 in a larger neighborhood.It approximates to the direction of edge curve more accurate than the direction indicated by gradient.When 1 < p < 2, the energy dissipation process is adaptively controlled by |u y, t − u x, t | p−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 |u y, t − u x, t | p−2 , and the diffusion along the orthogonal direction is enhanced for larger |u y, t − u x, t | p−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.

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 h and Δt be the space and time steps, respectively, and let x 1i ; x 2j ih; jh be the grid points.Let u n i; j be an approximation of the function u nΔt; x 1i ; x 2j , with n ≥ 0. Equation 3.2 can be discretized as follows: In all numerical experiments, we choose the following kernel function:

4.2
The constant C > 0 is selected such that Ω J x dx 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 p 1.8 for the proposed interpolation method, k 0.0001 for the heat diffusion model, and time step Δt 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 |u n 1 − u n | 2 < 10 −6 , normally within a few decade iterations.
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 p-Laplacian evolution works better than other methods.

Conclusion
In this paper, a new image interpolation model based on TV and nonlocal p-Laplacian regularization is proposed.It combines the advantages of TV regularizer and nonlocal p-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 p-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.
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 p-Laplacian if p > 1 19 .The nonlocal p-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 x and time t depends on all the values of u in a larger neighborhood of x.More precisely, if u t, x is thought of as the density at the point x at time t,

Table 1 :
Comparison of different interpolation algorithms using PSNR values.

Table 2 :
Comparison of different interpolation algorithms using MSSIM values.