# A Bidirectional Flow Joint Sobolev Gradient for Image Interpolation

Authors: Yi Zhan; Sheng Jie Li; Meng Li
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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
asL2 gradient flows on integral functionals and then modify the inner product from L2 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.


---

## Body

## 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 discretizedL2 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 L2 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 edgeinstead of along the gradient direction to the isophote line (edge) and the tangent direction as in [7]. We consider the diffusion equations as L2 gradient flows on integral functionals and then modify the inner product from L2 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 imageu0 from a high-resolution image u can be formulated as
(1)u0=Hu,
where H is a sparse matrix that combines both filtering and downsampling processes. The image interpolation is to solve unknown image u 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]:
(2)minuE(u)=Jd(u,u0)+λJr(u),
where λ is a regularization parameter that controls the tradeoff between Jd and Jr. The data fidelity function Jd is generally formulated in the classical least-squares sense as Jd(u,u0)=(1/2)|Hu-u0|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)Jr(u)=∫Ω‍φ(|∇u(x)|)dx.
Using Euler’s equation, the minimizer of problem (3) is the steady-state solution of the nonlinear PDE given by
(4)ut=div(φ′(|∇u(x)|)|∇u(x)|∇u(x)).
A direct calculation verifies
(5)ut=φ′(|∇u(x)|)div(∇u|∇u|)+φ′′(|∇u|)∂2u∂η2.
The two second-order directional derivatives are calculated as follows:
(6)∂2u∂ξ2=ξ(uxxuxyuxyuyy)ξT=uxxux2-2uxyuxuy+uyyuy2|∇u|2,∂2u∂η2=η(uxxuxyuxyuyy)ηT=uxxux2+2uxyuxuy+uyyuy2|∇u|2.A simple and effective functionφ is φ(x)=x, which corresponds to TV regularization. This derived the image interpolation based on the TV regularization under the level set method [4]
(7)ut=λ|∇u|div(∇u|∇u|)+HT(Hu-u0).An analogue regularization of (5) was proposed by Belahmidi and Guichard based on the classical heat diffusion model [5], formulated as follows:
(8)ut=|∇u|div(∇u|∇u|)+g(|∇u|)∂2u∂η2-HT(Hu-u0),
where the function g(s) is typically defined as
(9)g(s)=11+(s/K)2
with K>0 a constant to be tuned for a particular application. The role of the diffusion coefficient g(|∇u|) 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:
(10)ut=-α|∇u|th(k|∇(Gσ*u)|(Gσ*u)ηη)+β1+l1uξξ2uξξ+l2|∇u|21+l2|∇u|2(u0-u),
where α, β, k, l1, l2 are positive constant numbers, th(x) is a hyperbolic tangent function, and Gσ is a Gaussian function with standard deviation σ.In all these models,∂2u/∂η2 is the second-order directional derivative in the direction of the gradient |∇u| and |∇u|div(∇u/|∇u|)=∂2u/∂ξ2 in the direction orthogonal to the gradient |∇u|. From the viewpoint of geometry, the evolution process in the artificial time t given by these models is seen as energy dissipation process in two orthogonal directions η and ξ [4]. The diffusion process of u(x,t) 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 HT, and it will introduce some blur. In (8), the forward diffusion process of the grey values along η conserves the edges by the stop function g(s). 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, φ(x)=arctan(x). Equation (3) is rewritten as
(11)Jr(u)=∫Ω‍arctan(|∇u(x)|k)dx.Thus, our energy functional for image interpolation is defined as(12)E(u)=∫Ω‍arctan(|∇u(x)|k)dx+12|Hu-u0|2.
The problem now is to minimize or at least find critical points of the energy functional E(u) for which the Euler-Lagrange equation is a nonlinear PDE. There are three approaches to minimize E(u): Newton’s method, evolving the parabolic PDE ut=-∇E(u) 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 E(u) similar to [15, 17].Define a differential operatorD:H1,2(Ω)→L2(Ω)3 by
(13)D(u)=(u∇u)
and let
(14)〈g,h〉H1,2(Ω)=∫Ω‍gh+〈∇g,∇h〉=〈Dg,Dh〉L2(Ω)3.
From (12) the Fréchet derivative of E is
(15)E′(u)h=limα→01α[E(u+ah)-E(u)]=∫Ω‍μk〈∇u,∇h〉|∇u|(k2+|∇u|2)+λHT(Hu-u0)h.
For u∈H1,2(Ω), E′(u) is a bounded linear functional on H1,2(Ω), and hence, by the Riesz representation theorem, there exists ∇SE(u)∈H1,2(Ω), termed the Sobolev gradient, such that
(16)E′(u)h=〈∇SE(u),h〉H1,2(Ω)=〈D∇SE(u),Dh〉L2(Ω)3=〈D*D∇SE(u),h〉L2(Ω)
for all h∈H1,2(Ω), where D* denotes the adjoint of D. For u∈H2,2(Ω), E′(u) is also bounded on L2(Ω) and is represented by an L2 gradient
(17)E′(u)h=〈∇E(u),h〉L2(Ω)
for all h∈H1,2(Ω). The two gradients are thus related by
(18)∇SE(u)=(D*D)-1∇E(u),
where
(19)D*D=(I-∇)(I∇)=I-Δ.
Since H2,2(Ω) is a dense subspace of H1,2(Ω), (18) can be extended to u∈H1,2(Ω) by continuity. In order to obtain an expression for ∇E(u), we integrate by parts in (15):(20)E′(u)h=∫Ω‍[-μ(k2+3|∇u|2|∇u|(k2+|∇u|2)2div(∇u|∇u|)-2|∇u|(k2+|∇u|2)2Δu)-λHT(Hu-u0)]h,with ∇u=0 on the boundary of Ω. Hence
(21)-∇E(u)=μ(k2+3|∇u|2|∇u|(k2+|∇u|2)2div(∇u|∇u|)-2|∇u|(k2+|∇u|2)2Δu)+λHT(Hu-u0).
We improve the variant coefficients of the divergence and Laplace terms in (21) to obtain terse formula and better interpolation effects as follows:
(22)-∇E(u)=α|∇u|div(∇u|∇u|)-β|∇u|k2+|∇u|2Δu+λHT(Hu-u0),
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, Δu>0 and the second term in (22) is negative, and the dissipation process given by (22) diffuses u backward. While in darker lateral, Δu<0 and the second term in (22) is positive, the dissipation process given by (22) diffuses u 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 h and Δt be the space and time steps, respectively, and let (x1i;x2j)=(ih;jh) be the grid points. Let un(i;j) be an approximation of the function u(nΔt;x1i;x2j), with n≥0. Then for each n>0 (given un),(1)
discretize by finite differences and computeGn:=-∇E(u) from (22),(2)
introduce the notationw (will correspond to w=(un+1-un)/Δt),(3)
solve(I-Δ)w=Gn with (∂w/∂n→)|∂Ω in w using the semi-implicit scheme as in [16]: start with w0=0 (or the previous w), iterate for l=1,2,… until reaching the steady-state solution w:
(23)wi,jl+1-{wi+1,jl-2wi,jl+1+wi-1,jlΔx2+wi,j+1l-2wi,jl+1+wi,j-1lΔy2}=Gi,jn,(4)
updateun+1=un+Δt·w. The parameters are fixed by hand at some fixed values. In general, parameter K 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.Koala image expanded by1.5×1.5. (a) LSS, (b) FBF (k=500,α=0.04,β=0.01,l1=0.005,l2=2×10-4); (c) BG (λ=0.3,K=1); (d) proposed (k=100,λ=45,α=5,β=0.5).
(a)
(b)
(c)
(d)Barbara (portion) expanded by5×5. (a) LSS, (b) FBF (k=5,α=0.05,β=0.05,l1=5×10-4,l2=2×10-4); (c) BG (λ=1,K=1); (d) proposed (k=1,λ=25,α=10,β=1.5).
(a)
(b)
(c)
(d)Old man image reduced and expanded by3×3. Top row: original, LR image, LSS. Bottom row: FBF (k=500, α=0.01, β=0.05, l1=5×10-5, l2=2×10-5); BG (λ=0.5, K=0.2); proposed (k=70,λ=25,α=5,β=3).
(a)
(b)
(c)
(d)
(e)
(f)Flower1 image reduced and expanded by2×2. Top row: original, NEDI, LSS. Bottom row: FBF (k=500,α=0.01,β=0.05,l1=5×10-5,l2=2×10-5); BG (λ=0.5,K=0.2); proposed (k=70,λ=25,α=5,β=3).
(a)
(b)
(c)
(d)
(e)
(f)Flower2 image reduced and expanded by2×2. Top row: original, NEDI, LSS. Bottom row: FBF (k=500,α=0.01,β=0.05,l1=5×10-5,l2=2×10-5); BG (λ=0.5,K=0.2); proposed (k=70,λ=25,α=5,β=3).
(a)
(b)
(c)
(d)
(e)
(f)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 (Figure1), 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 t given by (22). The diffusion process of u(x,t) 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 of2×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.
Image
NEDI
LSS
FBF
BG
The proposed
Parrot (×3)
—
22.9104
25.9583
27.6972
27.8102
Old man (×3)
—
23.9840
29.2953
33.8563
33.9642
Flower1 (×2)
26.7876
28.4274
27.9264
30.4265
30.3702
Flower1 (×2)
25.9601
27.1878
26.4071
28.6159
28.7703


## 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.



---

*Source: 571052-2013-05-25.xml*