A Bidirectional Flow Joint Sobolev Gradient for Image Interpolation

Abstract

An energy functional with bidirectional flow is presented to sharpen image by reducing its edge width, which performs a forward diffusion in brighter lateral on edge ramp and backward diffusion that proceeds in darker lateral. We first consider the diffusion equations as 𝐿₂ gradient flows on integral functionals and then modify the inner product from 𝐿₂ to a Sobolev inner product. The experimental results demonstrate that our model efficiently reconstructs the real image, leading to a natural interpolation with reduced blurring, staircase artifacts and preserving better the texture features of image.

1. Introduction

Digital image interpolation is an inverse process of imaging process which samples higher resolution scene to lower resolution lattices. It is an important technology in digital image processing, artificial intelligence, pattern recognition, and modern industries. Producing visually pleasing resulting image is the key task in image interpolation. Many methods do their service for these problems, usually known as edge directed, level-set based, or isophote oriented, varying considerably.

The methods based on partial differential equations (PDEs) or energy variational usually have better effects in deblurring. Morse and Schwartzwald [1, 2] reduced the zigzagging artifacts by restricting curvature of the interpolated isophotes (equi-intensity contours). Minimum curvature is required on isophotes of the interpolated images. Malgouyres and Guichard [3] considered the preservation of one-dimensional structure to lead to the resulting images without blurring effects. Aly and Dubois [4] proposed a model-based TV regularization image up-sampling methods. Image acquisition process is modelled 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 oversmoothed homogeneous regions. Belahmidi and Guichard [5] have improved the TV-based interpolation by developing a nonlinear anisotropic PDE, hereafter referred to as BG interpolation method. Similar to [6], the nonlinear diffusion coefficient is locally adjusted according to image features such as edges, textures, and moments, which performs forward diffusion near homogeneous regions while backward diffusion near edges with strength and orientation adapted to image structures. To sharpen edges further, bidirectional diffusion (flow) was proposed [7–9]. In [7], a bidirectional flow is presented to sharpen image by reducing its edge width, which performs a fuzzy backward (inverse) diffusion along the gradient direction to the isophote line (edge), while it does a certain forward diffusion along the tangent direction on the contrary. A fractional-order bidirectional flow was also present in [9].

Generally, the solution of the energy functional is obtained by a steepest descent iteration with the discretized 𝐿₂ gradient. Richardson [10] demonstrated substantial benefits of Sobolev gradient method to a variety of PDEs that arise in image processing in computational efficiency. Sundaramoorthi et al. [11] presented Sobolev active contours with a Sobolev gradient in a space of smooth curves. Kazemi and Danaila [12] used Sobolev gradient method in conjunction with the steady-state solution of the Navier-Stokes equation in order to fill in missing pieces in a digital image. Yuan and He [13] proposed variational level set methods for image segmentation based on both 𝐿₂ and Sobolev gradients.

In this paper we focus on a new method for image interpolation based on variational regularization and Sobolev gradient. The variational regularization is a bidirectional flow in which forward diffusion and backward diffusion occur on both sides of image edge instead of along the gradient direction to the isophote line (edge) and the tangent direction as in [7]. We consider the diffusion equations as 𝐿₂ gradient flows on integral functionals and then modify the inner product from 𝐿₂ to a Sobolev inner product. The rest of the paper is organized as follows. In Section 2, we review the method based on variational method, PDEs method, and bidirectional diffusion method. The proposed variational regularization and Sobolev gradient for image interpolation are presented in Section 3. In Section 4 we demonstrate the experimental results on natural images which suggest that this method outperforms previous work, and Section 5 is for conclusion.

2. Background

The imaging processing of a digital lower resolution image 𝑢₀ from a high-resolution image 𝑢 can be formulated as 𝑢₀=𝐻𝑢, (1) where 𝐻 is a sparse matrix that combines both filtering and downsampling processes. The image interpolation is to solve unknown image 𝑢 from the ill-posed inverse problem (1). It is an efficient approach to reduce it to a minimization problem in a regularization-based framework as follows [4]: min_𝑢𝐸(𝑢)=𝐽_𝑑𝑢,𝑢₀+𝜆𝐽_𝑟(𝑢), (2) where 𝜆 is a regularization parameter that controls the tradeoff between 𝐽_𝑑 and 𝐽_𝑟. The data fidelity function 𝐽_𝑑 is generally formulated in the classical least-squares sense as 𝐽_𝑑(𝑢,𝑢₀)=(1/2)|𝐻𝑢−𝑢₀|². The variational regularization has been widely used in many areas of image processing [4, 7, 14]. The main advantage of the variational regularization is that it allows incorporating various prior knowledge, such as shape and intensity distribution, smoothing image edges while preserving important features of image. All of these formulations are different forms of a unified functional form of the variational regularizer given by 𝐽_𝑟(𝑢)=_Ω𝜑(|∇𝑢(𝑥)|)𝑑𝑥. (3) Using Euler’s equation, the minimizer of problem (3) is the steady-state solution of the nonlinear PDE given by 𝑢_𝑡𝜑=div^′(|∇𝑢(𝑥)|)|(𝑥(𝑥).∇𝑢)|∇𝑢 (4) A direct calculation verifies 𝑢_𝑡=𝜑^′(|∇𝑢(𝑥)|)div∇𝑢|∇𝑢|+𝜑^′′𝜕(|∇𝑢|)²𝑢𝜕𝜂². (5) The two second-order directional derivatives are calculated as follows: 𝜕²𝑢𝜕𝜉²𝑢=𝜉^𝑥𝑥𝑢_𝑥𝑦𝑢_𝑥𝑦𝑢_𝑦𝑦𝜉^𝑇=𝑢_𝑥𝑥𝑢²_𝑥−2𝑢_𝑥𝑦𝑢_𝑥𝑢_𝑦+𝑢_𝑦𝑦𝑢²_𝑦|∇𝑢|²,𝜕²𝑢𝜕𝜂²𝑢=𝜂^𝑥𝑥𝑢_𝑥𝑦𝑢_𝑥𝑦𝑢_𝑦𝑦𝜂^𝑇=𝑢_𝑥𝑥𝑢²_𝑥+2𝑢_𝑥𝑦𝑢_𝑥𝑢_𝑦+𝑢_𝑦𝑦𝑢²_𝑦|∇𝑢|². (6)

A simple and effective function 𝜑 is 𝜑(𝑥)=𝑥, which corresponds to TV regularization. This derived the image interpolation based on the TV regularization under the level set method [4] 𝑢_𝑡=𝜆|∇𝑢|div∇𝑢|∇𝑢|+𝐻^𝑇𝐻𝑢−𝑢₀. (7)

An analogue regularization of (5) was proposed by Belahmidi and Guichard based on the classical heat diffusion model [5], formulated as follows: 𝑢_𝑡=|∇𝑢|div∇𝑢||𝜕∇𝑢+𝑔(|∇𝑢|)²𝑢𝜕𝜂²−𝐻^𝑇𝐻𝑢−𝑢₀, (8) where the function 𝑔(𝑠) is typically defined as 𝑔(𝑠)=11+(𝑠/𝐾)² (9) with 𝐾>0 a constant to be tuned for a particular application. The role of the diffusion coefficient 𝑔(|∇𝑢|) is to control the smoothing adaptively.

In order to sharpen image by reducing its edge width, a fuzzy bidirectional flow framework based on generalized fuzzy set has been presented for image denoising by [7] as follows: 𝑢_𝑡=−𝛼|∇𝑢|th𝑘∇𝐺_𝜎∗𝑢𝐺_𝜎∗𝑢_𝜂𝜂+𝛽1+𝑙₁𝑢²_𝜉𝜉𝑢_𝜉𝜉+𝑙₂|∇𝑢|²1+𝑙₂|∇𝑢|²𝑢₀−𝑢, (10) where 𝛼, 𝛽, 𝑘, 𝑙₁, 𝑙₂ are positive constant numbers, th(𝑥) is a hyperbolic tangent function, and 𝐺_𝜎 is a Gaussian function with standard deviation 𝜎.

In all these models, 𝜕²𝑢/𝜕𝜂² is the second-order directional derivative in the direction of the gradient |∇𝑢| and |∇𝑢|div(∇𝑢/|∇𝑢|)=𝜕²𝑢/𝜕𝜉² in the direction orthogonal to the gradient |∇𝑢|. From the viewpoint of geometry, the evolution process in the artificial time 𝑡 given by these models is seen as energy dissipation process in two orthogonal directions 𝜂 and 𝜉 [4]. The diffusion process of 𝑢(𝑥,𝑡) along with 𝜉 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 [4, 5]. But this process cannot decrease the width of edges caused by the upsampling process 𝐻^𝑇, and it will introduce some blur. In (8), the forward diffusion process of the grey values along 𝜂 conserves the edges by the stop function 𝑔(𝑠). Unfortunately, the ability of the forward diffusion term in (8) is limited since it is a unilateral diffusion. The feature-dependent fuzzy bidirectional flow (10) has better effects in decreasing the width of edges since it is a bilateral diffusion across the image counters. However, the second-order directional derivative in the direction of gradient locate inaccurately the center of the “slop” of image edges, which reduces the ability of decreasing the width of edges. In the next section, we chose an appropriate function 𝜑 and improved its Sobolev gradient to decrease the width of edges.

3. The Variational Regularization and Sobolev Gradient

If we consider the gray-level value at a pixel to be analogous to the amount of particles, each having one unit of “mass,” stacked at the pixel, then in order to decrease the width of edges, we would like to move particles from the lower part of a “slope” upwards and from the higher part downwards. This process enhances the ratio of gray level on both sides of image edges. It can be achieved by enlarging the small and middle gradients in interpolation process. Since the arc tangent function rapidly increases near zero while very slowly on larger arguments, we take the arc tangent function for 𝜑 in (3); that is, 𝜑(𝑥)=arctan(𝑥). Equation (3) is rewritten as 𝐽_𝑟(𝑢)=_Ω|(𝑥arctan∇𝑢)|𝑘𝑑𝑥. (11)

Thus, our energy functional for image interpolation is defined as 𝐸(𝑢)=_Ω|(𝑥arctan∇𝑢)|𝑘1𝑑𝑥+2𝐻𝑢−𝑢₀². (12) The problem now is to minimize or at least find critical points of the energy functional 𝐸(𝑢) for which the Euler-Lagrange equation is a nonlinear PDE. There are three approaches to minimize 𝐸(𝑢): Newton’s method, evolving the parabolic PDE 𝑢_𝑡=−∇𝐸(𝑢) to steady state, and steepest descent iteration [15]. As pointed in [15], these methods all have some drawbacks. With the Sobolev gradient, however, steepest descent is often very effective [16]. Next we derive the Sobolev gradient of 𝐸(𝑢) similar to [15, 17].

Define a differential operator 𝐷∶𝐻^1,2(Ω)→𝐿₂(Ω)³ by 𝐷(𝑢)=𝑢∇𝑢 (13) and let 𝑔,ℎ_𝐻1,2(Ω)=_Ω𝑔ℎ+∇𝑔,∇ℎ=𝐷𝑔,𝐷ℎ_𝐿2(Ω)³. (14) From (12) the Fréchet derivative of 𝐸 is 𝐸^′(𝑢)ℎ=lim_𝛼→01𝛼[𝐸()(𝑢=𝑢+𝑎ℎ−𝐸)]_Ω𝜇𝑘⟨∇𝑢,∇ℎ⟩|∇𝑢|𝑘²+|∇𝑢|²+𝜆𝐻^𝑇𝐻𝑢−𝑢₀ℎ. (15) For 𝑢∈𝐻^1,2(Ω), 𝐸^′(𝑢) is a bounded linear functional on 𝐻^1,2(Ω), and hence, by the Riesz representation theorem, there exists ∇_𝑆𝐸(𝑢)∈𝐻^1,2(Ω), termed the Sobolev gradient, such that 𝐸^′(𝑢)ℎ=∇_𝑆𝐸(𝑢),ℎ_𝐻1,2(Ω)=𝐷∇_𝑆𝐸(𝑢),𝐷ℎ_𝐿2(Ω)³=𝐷^∗𝐷∇_𝑆𝐸(𝑢),ℎ_𝐿2(Ω) (16) for all ℎ∈𝐻^1,2(Ω), where 𝐷^∗ denotes the adjoint of 𝐷. For 𝑢∈𝐻^2,2(Ω), 𝐸^′(𝑢) is also bounded on 𝐿₂(Ω) and is represented by an 𝐿₂ gradient 𝐸^′(𝑢)ℎ=⟨∇𝐸(𝑢),ℎ⟩_𝐿2(Ω) (17) for all ℎ∈𝐻^1,2(Ω). The two gradients are thus related by ∇_𝑆𝐸(𝑢)=𝐷^∗𝐷⁻¹∇𝐸(𝑢), (18) where 𝐷^∗()𝐼∇𝐷=𝐼−∇=𝐼−Δ. (19) Since 𝐻^2,2(Ω) is a dense subspace of 𝐻^1,2(Ω), (18) can be extended to 𝑢∈𝐻^1,2(Ω) by continuity. In order to obtain an expression for ∇𝐸(𝑢), we integrate by parts in (15): 𝐸^′(𝑢)ℎ=∫_Ω𝑘−𝜇²+3|∇𝑢|²|∇𝑢|𝑘²+|∇𝑢|²²div∇𝑢||−2||∇𝑢∇𝑢𝑘²+|∇𝑢|²²Δ𝑢−𝜆𝐻^𝑇𝐻𝑢−𝑢₀ℎ, (20) with ∇𝑢=0 on the boundary of Ω. Hence (𝑢)𝑘−∇𝐸=𝜇²+3|∇𝑢|²|∇𝑢|𝑘²+|∇𝑢|²²div∇𝑢||−2||∇𝑢∇𝑢𝑘²+|∇𝑢|²²Δ𝑢+𝜆𝐻^𝑇𝐻𝑢−𝑢₀. (21) We improve the variant coefficients of the divergence and Laplace terms in (21) to obtain terse formula and better interpolation effects as follows: −∇𝐸(𝑢)=𝛼|∇𝑢|div∇𝑢||||∇𝑢−𝛽∇𝑢𝑘²+|∇𝑢|²Δ𝑢+𝜆𝐻^𝑇𝐻𝑢−𝑢₀, (22) where the parameter 𝛼, 𝛽 controls the tradeoff between the diffusion along the direction orthogonal to the gradient and the bidirectional diffusion.

Now we illustrate the action of the regularization terms in (22). As mentioned at the outset of this section, in order to decrease the width of edges, we would like to move particles from the lower part of a “slope” upwards and from the higher part downwards. It can be achieved by an energy dissipation process in two opposite directions, forward diffusion in brighter lateral on edge ramp and backward in darker lateral. In brighter lateral on edge ramp, Δ𝑢>0 and the second term in (22) is negative, and the dissipation process given by (22) diffuses 𝑢 backward. While in darker lateral, Δ𝑢<0 and the second term in (22) is positive, the dissipation process given by (22) diffuses 𝑢 forward.

4. Numerical Algorithm and Experimental Results

In this section, we develop a fully discrete numerical method to approximate problem (22). 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. Then for each 𝑛>0 (given 𝑢^𝑛),(1)discretize by finite differences and compute 𝐺^𝑛∶=−∇𝐸(𝑢) from (22),(2)introduce the notation 𝑤 (will correspond to 𝑤=(𝑢^𝑛+1−𝑢^𝑛)/Δ𝑡),(3)solve (𝐼−Δ)𝑤=𝐺^𝑛 with (𝜕𝑤/𝜕𝑛)|_𝜕Ω in 𝑤 using the semi-implicit scheme as in [16]: start with 𝑤⁰=0 (or the previous 𝑤), iterate for 𝑙=1,2,… until reaching the steady-state solution 𝑤: 𝑤^𝑙+1_𝑖,𝑗−𝑤^𝑙_𝑖+1,𝑗−2𝑤^𝑙+1_𝑖,𝑗+𝑤^𝑙_{𝑖−1,𝑗}Δ𝑥²+𝑤^𝑙_𝑖,𝑗+1−2𝑤^𝑙+1_𝑖,𝑗+𝑤^𝑙_{𝑖,𝑗−1}Δ𝑦²=𝐺^𝑛_𝑖,𝑗, (23)(4)update 𝑢^𝑛+1=𝑢^𝑛+Δ𝑡⋅𝑤. The parameters are fixed by hand at some fixed values. In general, parameter 𝐾 varies largely according to image features while 𝛼,𝛽,𝜆 in a certain range ([2,10], [0.5,5], [20,50]) are better.

We evaluate the performance of the proposed interpolation method in comparison with some representative work in the literatures: new edge-directed interpolation (NEDI) [18], the BG [5], the fuzzy bidirectional flow method (FBF) [7], and the local self-similarity algorithm (LSS) [19]. Experiments results are shown in Figures 1–4. 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. The first experiment directly interpolates a koala image and a Barbara image (portion) by a factor of 1.5×1.5 and 5×5, respectively, shown in Figures 1 and 2. In the second experiment, the images were first lowpass filtered and subsampled by a factor of 2×2 or 3×3; then the subsampled image was interpolated to the original image size, shown in Figures 3, 4, and 5.

(a)

(b)

(c)

(d)

Figure 1 

Koala image expanded by 1.5×1.5. (a) LSS, (b) FBF (𝑘=500,𝛼=0.04,𝛽=0.01,𝑙1=0.005,𝑙2=2×10−4); (c) BG (𝜆=0.3,𝐾=1); (d) proposed (𝑘=100,𝜆=45,𝛼=5,𝛽=0.5).

(a)

(b)

(c)

(d)

Figure 2 

Barbara (portion) expanded by 5×5. (a) LSS, (b) FBF (𝑘=5,𝛼=0.05,𝛽=0.05,𝑙1=5×10−4,𝑙2=2×10−4); (c) BG (𝜆=1,𝐾=1); (d) proposed (𝑘=1,𝜆=25,𝛼=10,𝛽=1.5).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Old man image reduced and expanded by 3×3. Top row: original, LR image, LSS. Bottom row: FBF (𝑘=500, 𝛼=0.01, 𝛽=0.05, 𝑙1=5×10−5, 𝑙2=2×10−5); BG (𝜆=0.5, 𝐾=0.2); proposed (𝑘=70,𝜆=25,𝛼=5,𝛽=3).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Flower1 image reduced and expanded by 2×2. Top row: original, NEDI, LSS. Bottom row: FBF (𝑘=500,𝛼=0.01,𝛽=0.05,𝑙1=5×10−5,𝑙2=2×10−5); BG (𝜆=0.5,𝐾=0.2); proposed (𝑘=70,𝜆=25,𝛼=5,𝛽=3).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

Flower2 image reduced and expanded by 2×2. Top row: original, NEDI, LSS. Bottom row: FBF (𝑘=500,𝛼=0.01,𝛽=0.05,𝑙1=5×10−5,𝑙2=2×10−5); BG (𝜆=0.5,𝐾=0.2); proposed (𝑘=70,𝜆=25,𝛼=5,𝛽=3).

From these figures, we can see that the LSS method sharpens the edges of image but smoothes unduly the image features, which produces the cartoon effects in result images. For example, the blocks or lumps of irregular shape appear in the fur of koala in face and forelimb (Figure 1), the strip (Figure 2), in face of old man and on the edge of hat (Figure 3), and in the petal (the detail texture disappeared in Figures 4 and 5). It is remarkable especially for larger expanded factor as shown in Figures 2 and 3, which results from using several small filters for larger factor. For example, to achieve a magnification of 5 it uses 5 : 4, 4 : 3, 4 : 3, 3 : 2, and 3 : 2. However, the effects appear even if only one filter is adopted, as shown in Figure 1. The FBF sharpens the image edges but cannot obtain smooth edges, unduly smooth image homogeneous areas, as shown in Figures 1 and 2. The BG produces better results than the LSS and FBF; the blur visual effect is still remarkable. Our approach not only sharpens the edges obtaining crisp and clear image edges not only, but also preserves the image detail feature resulting in the best visual effects. It results from the evolution process in the artificial time 𝑡 given by (22). The diffusion process of 𝑢(𝑥,𝑡) along the orthogonal direction to the gradient smooth image contours along them while the bidirectional flow decreases the width of edges maintaining their crispness. For example, the fur of koala, the fuzz on the hat of the old man, the stamen, and the stria in petals are seen clearly.

We use the classic PSNR to characterize the difference between the reference image and the output of a method. We use several test images including parrot, old man, and two flower images. To show the true power of the interpolation algorithms, we first downsampled each image by a factor of 2×2 or 3×3 and then interpolated the result back to its original size. The PSNR is shown in Table 1. From the table, the proposed method yields improved PSNR in all the experiments. This improvement may be attributed to the fact that our approach works better than other methods.

Table 1 

Comparison of different interpolation algorithms using PSNR values.

5. Conclusion

In this paper, a new image interpolation model based on energy variational is proposed. It diffuses forward and backward simultaneously on both sides of the image’s edge. A Sobolev gradient method is adopted to minimize the functional. We have shown improvement over our method on a subjective scale and in many cases with an improvement in PSNR.

Acknowledgments

This work was supported by Natural Science Foundation Project of CQ CSTC (Grant no. cstcjjA40012), Education Committee Project Research Foundation of Chongqing (no. KJ120709, no. KJ131209), the National Natural Science Foundation of China under Grant nos. 61202349 and 11226268, and the Research Foundation of Chongqing University of Arts and Sciences (no. R2012SC20).