Image Interpolation with Contour Stencils

Pascal Getreuer

doi:10.5201/ipol.2011.g_iics

Pascal Getreuer

→ BibTeX

@article{ipol.2011.g_iics,
    title   = {{Image Interpolation with Contour Stencils}},
    author  = {Getreuer, Pascal},
    journal = {{Image Processing On Line}},
    volume  = {1},
    year    = {2011},
    doi     = {10.5201/ipol.2011.g_iics},
}
% if your bibliography style doesn't support doi fields:
    note    = {\url{http://dx.doi.org/10.5201/ipol.2011.g_iics}}

published: 2011-08-01
reference: Pascal Getreuer, Image Interpolation with Contour Stencils, Image Processing On Line, 1 (2011). http://dx.doi.org/10.5201/ipol.2011.g_iics

Communicated by François Malgouyres
Demo edited by Pascal Getreuer

Pascal Getreuer pascal.getreuer@cmla.ens-cachan.fr, CMLA, ENS Cachan

Overview

Image interpolation is the problem of increasing the resolution of an image. Linear methods have traditionally been preferred, for example, the popular bilinear and bicubic interpolations are linear methods. However, a linear method must compromise between artifacts like jagged edges, blurring, and overshoot (halo) artifacts. These artifacts cannot all be eliminated simultaneously while maintaining linearity.

More recent works consider nonlinear methods, especially to improve interpolation of edges and textures. An important aspect of nonlinear interpolation is accurate estimation of edge orientations. For this purpose we apply contour stencils, a new method for estimating the image contours based on total variation along curves. This estimation is then used to construct a fast edge-adaptive interpolation.

References

F. Malgouyres and F. Guichard. “Edge direction preserving image zooming: A mathematical and numerical analysis.” SIAM J. Numerical Analysis, 39 (2002), pp. 1–37.
D. Muresan. “Fast edge directed polynomial interpolation.” ICIP (2), 2005, pp. 990–993.
A. Roussos and P. Maragos. “Reversible interpolation of vectorial images by an anisotropic diffusion-projection PDE.” Int. J. Comput. Vision, 2008.
P. Getreuer. “Contour stencils for edge-adaptive image interpolation.” Proc. SPIE, vol. 7257, 2009.
P. Getreuer. “Image zooming with contour stencils.” Proc. SPIE, vol. 7246, 2009.
onOne software. “Genuine Fractals.” www.ononesoftware.com

Online Demo

An online demo of this algorithm is available.

Contour Stencils

The idea in contour stencils is to estimate the image contours by measuring the total variation of the image along curves. Define the total variation (TV) along curve C

$\|u\|_{\mathrm{TV}(C)} = \int_0^T \bigl|\tfrac{\partial}{\partial t} u\bigl(\gamma(t)\bigr)\bigr| \,dt, \quad \gamma:[0,T]\rightarrow C$

where γ is a smooth parameterization of C. The quantity ||u||_TV(C) can be used to estimate the image contours. If ||u||_TV(C) is small, it suggests that C is close a contour. The contour stencils strategy is to estimate the image contours by testing the TV along a set of candidate curves, the curves with small ||u||_TV(C) are then identified as approximate contours.

Contour stencils are a discretization of TV along contours. As described in [4], [5], a ''contour stencil'' is a function S : Z^2 x Z^2 -> R^+ describing weighted edges between the pixels of v. Stencil is applied to v at pixel k ∈ Z^2 as

$(\mathcal{S} \star [v])(k) := \sum_{m,n\in\mathbb{Z}^2} \mathcal{S}(m,n) |v_{k+m} - v_{k+n}|.$

Defining [v](m,n) := |v_m − v_n|, this quantity is (with an abuse of notation) a cross-correlation over Z^2 x Z^2 evaluated at (k,k). The stencil edges are used to approximate a curve C so that the quantity $(S \star [v])(k)$ approximates ||u||_TV(C+k) (where C + k := {x+k : x ∈ C}). The image contours are estimated by finding a stencil with small TV. The best-fitting stencil at pixel k is

$\mathcal{S}^\star(k) = \operatorname*{arg\,min}_{\mathcal{S}\in\Sigma} (\mathcal{S} \star [v])(k)$

where Σ is a set of candidate stencils. It is possible that the minimizer is not unique, for example in a locally constant region of the image. For simplicity, we do not treat this situation specially and always choose a minimizer even if it is not unique. This best-fitting stencil provides a model of the image contours in the neighborhood of pixel k. The stencils used in this work are shown below.

The proposed stencil set Σ. The edge weights are denoted by superscript α, β, δ, γ.

For the set of candidate stencils Σ, we use 8 line-shaped stencils that were designed to distinguish between the functions

$f^j(x) = x_1 \sin \tfrac{\pi}{8}j - x_2 \cos \tfrac{\pi}{8}j, \quad j = 0, \ldots, 7.$

The edge weights α, β, δ, γ are selected so that on the function f(x) = x₁sinθ − x₂cosθ,

$|theta - pi/16 j| < pi/16 implies S^{pi/8 j} = argmin_{S in Sigma} (S star [f]).$

In this way, the stencils can fairly distinguish 8 different orientations. An estimate of the local contour orientation at point k is obtained by noting which stencil is the best-fitting stencil .

Normalized stencil total variations $(\mathcal{S}^{\frac{\pi}{8}j} \star [f])$ vs. θ.
Left: The first three stencils, j = 0, 1, 2. Right: All eight stencils.

For a color image, the image is converted from RGB to a luma+chroma space

[Y U V]^T = [2 4 1; -1 -2 3; 4 -3 -1] * [R G B]^T

and the stencil TV is computed as the sum of $(S \star [v])(k)$ applied to each color channel.

Interpolation

Given image v known on Z^2 , we seek to construct an image u on R^2 such that

$v(x) = (h \ast u)(x), \quad \text{for all } x\in\mathbb{Z}^2.$

where h is the (assumed known) point spread function and ∗ denotes convolution.

The goal is to incorporate deconvolution yet maintain computational efficiency. To achieve this, the global operation of deconvolution is approximated as a local one, such that pixels only interact within a small window.

Local Reconstructions

For every pixel k in the input image, we begin by forming a local reconstruction

$u_k(x) = v_k + \sum_{n\in\mathcal{N}} c_n \,\varphi_{\mathcal{S}^\star(k)}(x - n),$

where ⊂ Z^2 is a neighborhood of the origin and $phi^n_{S*(k)}$ is a Gaussian oriented with the contour modeled by the best-fitting stencil .

Local reconstruction u_k with oriented functions functions .

The c_n are chosen such that u_k satisfies the discretization model locally,

$(h \ast u_k)(m) = v_{k+m}, \quad m\in\mathcal{N}.$

This condition implies that the c_n satisfy the linear system

$\sum_n (A_\mathcal{S})_{m,n} c_n = v_{k+m} - v_k, \quad n \in \mathcal{N},$

where is a matrix with elements $(A_S)_{m,n} = (h * phi_S)(m - n)$ . By defining the functions

$\psi_\mathcal{S}^n(x) := \sum_{m\in\mathcal{N}} (A_\mathcal{S}^{-1})_{m,n} \varphi_\mathcal{S}(x-m),$

u_k can be expressed directly in terms of the samples of v,

$u_k(x) = v_k + \sum_{n\in\mathcal{N}} (v_{k+n} - v_k)\psi_{\mathcal{S}^\star(k)}^n(x).$

Global Reconstruction

The u_k are combined with overlapping windows to produce the interpolated image,

$u(x) = \sum_{k \in \mathbb{Z}^2} w(x-k) u_k(x-k).$

The window should satisfy ∑_k w(x − k) = 1 for all x ∈ R^2 and w(k) = 0 for k ∈ Z^2 \ .

Iterative Refinement

This global reconstruction satisfies the discretization model approximately, ↓(h ∗ u) ≈ v. The accuracy may be improved using the method of iterative refinement. Let denote the global reconstruction formula,

$(\mathcal{R}v)(x) := \sum_{k\in\mathbb{Z}^2} w(x - k)\Bigl[ v_k + \sum_{n\in\mathcal{N}} (v_{k+n} - v_k) \psi^n_{\mathcal{S}^\star(k)}(x - k) \Bigr]$

such that u = v (where we consider as fixed parameters so that is a linear operator). Then the deconvolution accuracy is improved by the iteration

$u^{0} = 0, \quad u^{i+1} = u^{i} + \mathcal{R}\bigl(v - \operatorname{sample}(h*u^{i})\bigr).$

Each iteration should reduce the residual in satisfying the discretization model,

$r^{i} = v - \operatorname{sample}(h*u^{i}).$

The residual reduces quickly in practice, usually three or four iterations is sufficient for accurate results.

Parameters

The following parameters are fixed in the experiments:

h is a Gaussian with standard deviation 0.5,
= {-1,0,1}×{-1,0,1},
w is the cubic B-spline,
$w(x,y) = B(x)B(y), \quad B(t) = \bigl(1-|t| + \tfrac{1}{6} |t|^3 - \tfrac{1}{3}\bigl|1-|t|\bigr|^3\bigr)^+$
is an oriented Gaussian,
$\varphi_\mathcal{S}(x,y) = \exp\Bigl( -\frac{\tau^2}{2\sigma_\tau^2} -\frac{\nu^2}{2\sigma_\nu^2} \Bigr), \quad \begin{pmatrix} \tau \\ \nu \end{pmatrix} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}$
with σ_τ = 1.2 and σ_ν = 0.6, and θ is the orientation modeled by ,

and three iterations of iterative refinement are applied (one initial interpolation and two correction passes).

For sake of demonstration, the examples below use a PSF with a substantial amount of blur, σ_h = 0.5. The default value for σ_h is 0.35 in the online demo associated with this article, which better models the blurriness of typical images.

Algorithm

The interpolation is computationally efficient. We first consider the complexity without iterative refinement.

The matrices can be precomputed for each stencil ∈Σ, allowing the c_n coefficients to be computed in 6² + 3 operations per (color) input pixel. Furthermore, since w has compact support, u only depends on the small number of u_k where w(x − k) is nonzero. Let W be a bound on the number of nonzero terms,

$\#\{k\in\mathbb{Z}^2 : w(x-k) \ne 0\} \le W, \quad \text{for all } x\in\mathbb{R}^2.$

We suppose that W is O(). Given the c_n, each evaluation of u(x) costs O(²) operations. So for factor-d scaling, the total computational cost is O(²d²) operations per input pixel. For scaling by rational d, samples of w and can also be precomputed, and scaling costs 6Wd² operations per input pixel. For the settings used in the examples, this is 864d² operations per input pixel.

With iterative refinement, the previous cost is multiplied by the number of steps and there is the additional cost of computing the residual. If h is quickly decaying, then it is accurately approximated by an FIR filter with O(d²) taps and the residual

$r^{i} = v - \operatorname{sample}(h*u^{i})$

can be computed in O(d²) operations per input pixel.

Implementation

This software is distributed under the terms of the simplified BSD license.

source code zip tar.gz
online documentation

Please see the readme.html file or the online documentation for details.

Implementation notes:

Fixed-point arithmetic is used to accelerate the main computations.
For efficiency in the correction passes of iterative refinement, the u_k for which |v_k| is small are not added (so that they do not need to be computed),

$u(x) = \sum_{k \in \mathbb{Z}^2:|v_k| \ge T} w(x-k) u_k(x-k).$

Examples

Here we perform an interpolation experiment to test the performance of the proposed interpolation strategy.

First, a high-resolution image u_o is smoothed and downsampled by factor 4 to obtain a coarsened image v = ↓(h ∗ u_o) where h is a Gaussian with standard deviation 0.5 in units of input pixels, σ_h = 0.5. This amount of smoothing is somewhat weak anti-aliasing, so the input data is slightly aliased.

The value of σ_h should estimate the blurriness of the PSF used to sample the input image. It is better to underestimate σ_h rather than overestimate: if σ_h is smaller than the true standard deviation of the PSF, the result is merely blurrier, but using σ_h slightly to large creates ripple artifacts.

The method works well for 0 ≤ σ_h ≤ 0.7. For σ_h above 0.7, the method produces visible ringing artifacts (even if the true PSF used to sample the input image has standard deviation σ_h). One could expect this effect, since there is no kind of regularization in the deconvolution. In the online demo, the default value for σ_h is 0.35, which reasonably models the blurriness of typical images.

Interpolation is then performed on v to produce u approximating the original image u_o. The interpolation and the original image are compared with the peak signal-to-noise ratio (PSNR) and mean structural similarity (MSSIM) metrics (How are these computed?).

Comparison with Other Methods

Original Image (332×300)	Input Image (83×75)

Estimated Contour Orientations	Contour Stencil Interpolation PSNR 25.77, MSSIM 0.7165, CPU time 0.109s

The following table shows the convergence of the residual rⁱ = v − ↓(h ∗ u) where the image intensity range is [0,1].

Iteration i	\|\|rⁱ\|\|_∞
1	0.05409007
2	0.01677390
3	0.00661765

For comparison, the same experiment is performed with standard bicubic interpolation, Muresan's AQua-2 edge-directed interpolation [2], Genuine Fractals fractal zooming [6], Fourier zero-padding with deconvolution, Malgouyres' TV minimization [1], and Roussos and Maragos' tensor-driven diffusion [3]. The first three of these methods do not take advantage of knowledge about the point spread function, while the later three do (notice their sharper appearance).

Bicubic PSNR 24.36, MSSIM 0.6311, CPU time 0.012s	AQua-2 [2] PSNR 23.97, MSSIM 0.6062, CPU time 0.016s

Fractal Zooming [6] PSNR 24.50, MSSIM 0.6317	Fourier Zero-Padding with Deconvolution PSNR 25.70, MSSIM 0.7104, CPU time 0.049s

TV Minimization [1] PSNR 25.87, MSSIM 0.7181, CPU Time 2.72s	Tensor-Driven Diffusion [3] PSNR 26.00, MSSIM 0.7297, CPU Time 5.11s

The contour stencil interpolation has good quality similar to tensor-driven diffusion but with an order magnitude lower computation time.

Geometric Features

The following experiment on a synthetic image tests the method's ability to handle different geometric features.

Original Image (320×240)	Input Image (80×60)

Estimated Contour Orientations	Contour Stencil Interpolation PSNR 21.23, MSSIM 0.8548, CPU time 0.078s

Textures

Because the method is sensitive to the image contours, oriented textures like hair can be reconstructed to some extent. Interpolation of rough textures with turbulent contours is less successful.

Original Image (392×304)	Original Image (332×304)

Input Image (98×76)	Input Image (83×76)

Contour Stencil Interpolation PSNR 33.24, MSSIM 0.7762, CPU time 0.129s	Contour Stencil Interpolation PSNR 22.47, MSSIM 0.6051, CPU time 0.115s

Noisy Images

A limitation of the method is the design assumption that noise in the input image is negligible. If noise is present, it is amplified by the deconvolution. The sensitivity to noise increases with the PSF standard deviation σ_h, which controls the deconvolution strength. Similarly, if σ_h is larger than the standard deviation of the true PSF that sampled the image, then the method produces significant oscillation artifacts because the deconvolution exaggerates the high frequencies.

Original Image (256×256)

The top row shows the input images and the bottom row shows their interpolations.

Clean Input	JPEG Compressed	Quantized Colors

Contour Stencil Interpolation PSNR 26.48, MSSIM 0.8196	Contour Stencil Interpolation PSNR 20.38, MSSIM 0.5244	Contour Stencil Interpolation PSNR 18.30, MSSIM 0.3393

This material is based upon work supported by the National Science Foundation under Award No. DMS-1004694. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Work partially supported by the MISS project of Centre National d'Etudes Spatiales, the Office of Naval Research under grant N00014-97-1-0839 and by the European Research Council, advanced grant “Twelve labours.”

image credits

John D. Willson, USGS Amphibian Research and Monitoring Initiative
Pascal Getreuer
Pascal Getreuer
Peppers standard test image
Gary Kramer, U.S. Fish and Wildlife Service

Maximum absolute difference	$:= \\|A - B\\|_\infty$ ,
Mean squared error (MSE)	$:= \tfrac{1}{3N}\\|A - B\\|_2^2$ ,
Root mean squared error (RMSE)	$:= \sqrt{\text{MSE}}$
	$\hphantom{:}= \tfrac{1}{\sqrt{3N}}\\|A - B\\|_2$ ,
Peak signal-to-noise ratio (PSNR)	$:= 10 \log_{10}\frac{255^2}{\text{MSE}}$
	$\hphantom{:}= 20 \log_{10}\frac{255}{\text{RMSE}}$ .

Overview

References

Online Demo

Contour Stencils

Interpolation

Local Reconstructions

Global Reconstruction

Iterative Refinement

Parameters

Algorithm

Implementation

Examples

Image Quality Metrics

L^p Metrics

MSSIM

Computation Time

Comparison with Other Methods

Geometric Features

Textures

Noisy Images

Overview

References

Online Demo

Contour Stencils

Interpolation

Local Reconstructions

Global Reconstruction

Iterative Refinement

Parameters

Algorithm

Implementation

Examples

Image Quality Metrics

Lp Metrics

MSSIM

Computation Time

Comparison with Other Methods

Geometric Features

Textures

Noisy Images

L^p Metrics