Structured Pursuits for Geometric Super-Resolution

Document Sample
Structured Pursuits for Geometric Super-Resolution Powered By Docstoc
					 Structured Pursuits for Geometric Super-Resolution
                                                St´ phane M ALLAT and Guoshen Y U
                                      CMAP, Ecole Polytechnique, 91128 Palaiseau Cedex, France

   Abstract—Super-resolution image zooming is possible when the                         II. D IRECTIONAL I NTERPOLATIONS
image has some geometric regularity. We introduce a general
class of non-linear inverse estimators, which combines linear
                                                                             Let us consider a subsampling operator U by a factor 2.
estimators with mixing weights in a frame providing a sparse               The measured low-resolution image y is
representation. Mixing weights are computed with a block decom-
position, which minimizes a Tikhonov energy penalized by an l1
                                                                                             y[n1 , n2 ] = f [2n1 , 2n2 ] + w[n1 , n2 ]              (2)
norm of the mixing weights. A fast orthogonal matching pursuit             and if we neglect the noise w its Fourier transform is
algorithm computes the mixing weights. Adaptive directional
image interpolations are calculated with mixing weights in a                               1 ˆ             ˆ
wavelet frame.
                                                                             y (ω1 , ω2 ) = f (ω1 , ω2 ) + f (ω1 + π, ω2 )
                    I. I NTRODUCTION                                                         ˆ                  ˆ
                                                                                           + f (ω1 , ω2 + π) + f (ω1 + π, ω2 + π) .                  (3)
  We consider linear inverse problems where the measure-
                                                                           Fig. 1(b) shows such an aliased Fourier transform.
ments y of a signal f are obtained through a non-invertible
operator U to which is added a noise w:
                         y = U f + w.
Image interpolation is an example of such ill-posed inverse
                                                                    −pi                                                                   −pi
problems, where U is a subsampling operator. This paper
studies non-linear inverse estimators, obtained as a mixture
of linear inverse estimators, for image interpolation.               0                                                                     0

   Increasing the size of images is needed for many image
or video screens. Most high definition televisions perform            pi
                                                                     −pi         0      pi
                                                                                                                                           −pi   0         pi

a scaling of input standard definition images, with a spatial                     a                     b                   c                     d
interpolation. When images are aliased, linear interpolations
introduce artifacts such as Gibbs oscillations or zigzag along             Fig. 1.    (a): a high-resolution edge image f [n] and its Fourier
                                                                           transform. (b): a low-resolution image y[n] = f [2n] and its aliased
edges. Better estimations can be obtained with non-linear                  Fourier transform. (c): 2D cubic-spline interpolation (50.45 dB). (d):
estimators which take into account local image properties.                 super-resolution estimation with directional interpolation along the
   Suppose that we are given a family of linear estimators                 contour direction (72.72 dB). The ellipse in (b) gives the frequency
 ˜                                                                         domain preserved by the directional interpolator.
fθ = Iθ y that are optimized with specific priors on f . For                               ˆ
                                                                             To recover f (ω1 , ω2 ), a filter h is adjusted to reduce the
image interpolation, Iθ is a directional interpolation along
                                                                           impact of the three aliased components:
a direction θ, which is precise if the image has regular
                                                                             f (n) = Iy(n) =        y(q) h(n − 2q) for n = (n1 , n2 ).
variations along θ. Mixing linear estimators is very effective
for noise removal [3]. However, these approaches do not apply
to ill-posed inverse problems because they rely on a Stein                 Interpolation filters satisfy h(0) = 1 and h(2n) = 0 for n = 0,
                                                                           so that f (2q) = f (2q) if w = 0. If f is uniformly regular
unbiased empirical estimator of the risk, and there is no such
estimator for ill-posed inverse problems. In the context of                then its Fourier transform has an energy mostly concentrated
image interpolation, this paper introduces a new approach to               in [−π/2, π/2]2 and h is chosen to be a low-pass filter.
define an adaptive inverse estimator which combines linear                  Among all linear filters, a cubic spline interpolation filter
estimators. The measurement y is decomposed in a basis or                  h = hs defines an interpolator Is which yield nearly the
a frame {φp }p∈G providing a sparse representation, and a                  best PSNR for most images [5]. When the frequency support
                                                                           of f exceeds [−π/2, π/2]2 then a low frequency interpolator
family of linear estimators {Iγ }γ∈Γ are combined with mixing
weights a(γ, p):                                                           produces Gibbs oscillations, zigzag jaggies along contours.
                                                                         Fig. 1(c) gives an example with a cubic spline filter.
                                                                              The optimization of h can be adapted to prior information
            f=         Iγ                     ˜
                                    a(γ, p) y, φp φp  .      (1)          on f with a quadratic Tikhonov regularization. Suppose that f
                 γ∈Γ          p∈G                                          has some form of regularity, which is measured by a regularity
Mixing weights are computed for directional image interpo-                 norm Rf , where Rf is typically a high-pass filtering.
lations with an adaptive block transform in a wavelet frame.               Sobolev norms are particular instances of such regularity
Directional blocks are adapted to the directional regularity of            norms. Let σ 2 be the variance of the noise w. A Tikhonov
wavelet coefficients, with a block matching pursuit algorithm.                                    ˜
                                                                           estimator computes f = Iy by
                          2                             2
  minimizing        RIy       subject to    U Iy − y≤ σ 2 . (4)
We shall neglect the noise for image interpolation, so U f = y.
Minimizing Rθ f                                          ˜
                   ˜ where Rθ computes variations of f in the
direction θ is a prior regularity condition which corresponds to
images where f (ω1 , ω2 ) has a frequency support concentrated
in a narrow band around a line of angle θ+π/2, as illustrated in
Fig. 1(a). The resulting Tikhonov interpolator Iθ in Fig. 1(d),           Fig. 2.  A block Bθ,q of orientation θ regroups wavelet coefficients
                                                                          that are aligned with the directional image regularity in a neighbor-
performs an interpolation in the direction θ with a filter                 hood of a point q, in the 3 wavelet coefficient images corresponding
hθ which eliminates the aliased frequency components while                to the 3 wavelet directions d = 1, 2, 3. Three blocks are shown in
preserving the support of f (ω1 , ω2 ).                                   this figure.
   The direction of image regularity changes in space and it is
therefore necessary to adapt locally the interpolation angle θ.           function of a block, which does not depend on d = 1, 2, 3.
Many researchers and engineers in industry have studied such              The block position q is sampled along a grid Gθ ⊂ G so
adaptive interpolators. The interpolation direction is adapted            that G = ∪q∈Gθ Bθ,q and each block overlaps with only two
along edges, using ad-hoc edge detectors [2] or gradient based            other blocks of same orientation. Among all Bθ,q , the blocks
orientation estimators [4], [6]. These algorithms can provide             selected by interpolator Iθ are located where the image is
better visual quality then a linear interpolation but the PSNR            regular in the direction θ, as illustrated in Fig. 2. This selection
is hardly improved.                                                       is performed by coefficients (θ, q). They should be close to
                                                                          1 when the image is regular in the direction θ and otherwise
     III. B LOCK A DAPTIVE T IKHONOV E STIMATIONS                         close to 0. Selection coefficients decompose yh according to
    Adaptive interpolation amounts to combine a family of                 its directional regularity plus a residue:
directional interpolators Iθ and the cubic spline interpolator
Is . This section introduces an algorithm which computes                                       yh =                    (θ, q) yθ,q + yr .                (6)
                                                                                                        θ∈Θ q∈Gθ
mixing weights in a wavelet frame, according to the general
formulation (1), with an adaptive block decomposition.                    The residue yr has no directional regularity and is therefore
    The subsampled image y[n] for n ∈ G is decomposed in a                interpolated with a cubic spline interpolator Is . The resulting
                                            d                             estimator is:
translation invariant wavelet tight frame {ψp }0≤d≤3,p∈G on 1
                                                                                                                      
scale. It separates low frequencies yl over the low frequency
                  0                                                                 ˜
scaling filters {ψp }p∈G and high frequencies yh over fine scale                      f = Is (yl ) +            Iθ           (θ, q) yθ,q  + Is (yr ).    (7)
wavelets {ψp }1≤d≤3,p∈G :                                                                              θ∈Θ           q∈Gθ

                                            3                             This interpolator is a mixture of the interpolators {Iθ , Is }θ∈Θ
                             0                                d
       yl =         c(0, p) ψp and yh =              c(d, p) ψp ,   (5)   where the image low frequency and the high frequency residue
              p∈G                          d=1 p∈G                        are interpolated with a cubic spline:
                                               ˜d                                                                   
where the wavelet coefficients c(d, p) = y, ψp are computed                                                3
with the dual wavelet frame. The low frequency image yl has                  ˜
                                                                             f=            Iθ                a(θ, p) c(d, p) ψp 
little aliasing and can thus be precisely interpolated with a                        θ∈Θ        p∈G d=1
cubic spline interpolator Is . The high frequency image yh is                                                                                       
interpolated by combining directional interpolators and a cubic                      + Is                     0
                                                                                                      c(0, p) ψp +                a(s, p) c(d, p) ψp  .
spline interpolator.                                                                           p∈G                     p∈G d=1
    Let Θ be a set of angles θ uniformly discretized between
0 and π. An adaptive interpolator is implemented over fine                 The mixing weights are given by
scale wavelets ψp , by mixing the interpolators {Iθ , Is }θ∈Θ             a(θ, p) =             (θ, q)1Bθ,q (p) and a(s, p) = 1−                   a(θ, p) .
with weights which depend on the image regularity. The ap-                              q∈Gθ                                                 θ∈Θ
pendix specifies the directional interpolators used in numerical                                                                     (9)
experiments.                                                                 Selection coefficients (θ, q) are computed by minimizing
    To construct a mixing operator, wavelets are regrouped in             a criterion which multiplies these coefficients by a Tikhonov
blocks for which a specific interpolator is chosen. Let us                 regularity term in the corresponding blocks, while minimizing
denote                                                                    the norm of the residue yr . The norm of the residue is
                   yθ,q =           c(d, p) ψp                            controlled by their wavelet coefficient energy. The regularity
                              (d,p)∈Bθ,q                                  in the direction θ is calculated with an operator Rθ applied
                                                                          to wavelet coefficients, which computes their variation in the
the image projection over wavelets in a spatial block Bθ,q .
                                                                          direction θ:
This block is a rectangle elongated in the direction θ, and
located in the spatial neighborhood of q for the three directions                      Rθ c    Bθ,q   =                |c(d, p) − Aθ c(d, p)|2 ,
d = 1, 2, 3. We write 1Bθ,q (d, p) = 1Bθ,p (p) the indicator                                              (d,p)∈Bθ,q
where Aθ c(d, p) is constant over each line of angle θ included                           The overall algorithm thus proceeds as follow:
in Bθ,q and equal to the average of c(d, p) over each of these                            • Compute the wavelet transform of y.
lines, for each d = 1, 2, 3. Selection coefficients are thus                               • Compute selection coefficients        with the orthogonal
computed by minimizing                                                                      pursuit.
                 3                                                                        • Apply the reconstruction (8) with the mixing weights (9).
            1                                                                        2
L( )   =                    |c(d, p)|2 1 −                         (θ, q) 1Bθ,q (p)       Fig. 3 compares a separable cubic spline interpolation with
                d=1 p∈G                          θ∈Θ q∈Gθ                              this adaptive interpolation computed by mixing directional
           +λ                       (θ, q) Rθ c      2
                                                     Bθ,q   ,                     (10) interpolations. The mixing algorithm improves the PSNR and
                 θ∈Θ q∈Gθ                                                              the visual image quality where the image is geometrically
                                                                                       regular. This clearly appears in the straws, the hat border, the
                         2                                                             cloth and the hairs of various directions.
                     c   Bθ,q   =                |c(d, p)|2 .
                                    (d,p)∈Bθ,q                                                                     A PPENDIX
The energy L( ) has a quadratic term in , which is penalized                              We describe a directional interpolator Iθ for upscaling
by a weighted l1 norm. It is minimized by standard numerical                           images by a factor 2, which is illustrated in Fig. 4. Origi-
algorithms such as an iterated thresholding [1].                                       nal samples are shown as crosses. Step 1 computes a one-
                                                                                       dimensional interpolations in the direction θ. We consider all
       IV. O RTHOGONAL B LOCK M ATCHING P URSUIT                                       lines of angle θ that intersect original image samples (crosses
   An l1 minimization is computationally too expensive for                             in Fig. 4) and we compute mid-points (circles) with a cubic
real time image processing applications. In the following, we                          spline interpolation. It oversamples by a factor two either
develop a fast orthogonal matching pursuit algorithm, which                            the image rows, or the image columns or the diagonals of
computes an approximate solution.                                                      angle ±π/4. The missing coefficients are shown as squares
   An orthogonal block matching pursuit is a greedy algorithm                          in Fig. 4. Step 2 computes new samples (dots) with a cubic
which selects progressively non-overlapping blocks which                               spline interpolation along these oversampled rows, columns or
minimize the energy (10). If a single block Bθ,q is selected                           diagonals. The position of these new samples are chosen so
then one can verify that (10) is minimized for                                         that any missing coefficient (square) is a mid-point between
                                             Rθ c        2                             two dots on a line of angle θ. Step 3 computes missing samples
               (θ, q) = max 1 − λ                    2          ,0 .            (11)   (squares) with a cubic spline linear interpolation along the
                                                 c   Bθ,q                              direction θ from the previously calculated new samples (dots).
Setting (θ, q) to this optimal value reduces L( ) by
                     2                                      2
                 c   Bθ,q                         Rθ c      Bθ,q        2
   e(θ, q) =                    max 1 − λ                2         ,0       .   (12)
                     2                               c   Bθ,q

The best block Bθ,q is the one which minimizes L( ) and
hence which maximizes e(θ, q).
   An orthogonal matching pursuit algorithm selects blocks                             Fig. 4.  Directional interpolation in three steps: a one-dimensional
that do not intersect and which reduce as much as possible                             interpolations along the angle θ, a vertical interpolation and another
L( ). At each iteration, among all blocks that do not intersect                        interpolation along θ.
with the previously selected blocks, it finds the block which
                                                                                                                  R EFERENCES
maximizes e(θ, q). It proceeds as follows.
   1) Initialize Ω0 = {(θ, q) : θ ∈ Θ , q ∈ Gθ }, (θ, q) = 0                           [1] I. Daubechies. and M. Defrise, and C. De Mol, “An Iterative
                                                                                           Thresholding Algorithm for Linear Inverse Problems with a
       for all (θ, q) ∈ Ω0 and l = 0.                                                      Sparsity Constraint”, Comm. on Pure Appl. Math., 57:1413–
   2) Find (θl , ql ) ∈ Ωl which maximizes e(θ, q) in (12). and                            1457, 2004.
       set (θl , ql ) with (11).                                                       [2] K. Jensen and D. Anastassiou, “Subpixel edge localization
   3) Set Ωl+1 = {(θ, q) ∈ Ωl : Bθ,q ∩ Bθl ,ql = 0}   /                                    and the interpolation of still images”, IEEE Trans. on Image
   4) If e(θl , ql ) > 0 set l = l + 1 and go back to (2).                                 Processing, 4(3):285–295, 1995.
                                                                                       [3] G. Leung and Barron, A.R., “Information theory and mixing
       Otherwise stop.                                                                     least-squares regressions”, IEEE Trans. on Information Theory,
   This algorithm computes selection coefficients (θ, q) for                                52(8):3396–3410, 2006.
all θ ∈ Θ and q ∈ Gθ . A translation invariant estimator                               [4] X. Li and M. T. Orchard, “New edge-directed interpolation”,
computes selection coefficients (θ, q) for all q ∈ G. This                                  IEEE Trans. on Image Processing, 10(10):1521–1527, 2001.
                                                                                       [5] M. Unser, “Splines: A perfect fit for signal and image process-
is done by covering the original image grid G by a union                                   ing”, IEEE Signal Processing Magazine, 16(6):22–38, 1999.
of translated grids Gθ,τ = Gθ + τ . The translation invariant                          [6] Q. Wang and R.K. Ward, “A New Orientation-Adaptive Interpo-
selection coefficients (θ, q) are the average of the selection                              lation Method” IEEE Trans. on Image Processing, 16(4):889–
coefficients calculated over each translated grid Gθ,τ , with the                           900, 2007.
orthogonal matching pursuit algorithm.
     High-resolution image                            Cubic spline 21.37 dB                      Adaptive interpolation 23.82 dB

     High-resolution image                            Cubic spline 34.02 dB                      Adaptive interpolation 35.84 dB

     High-resolution image                            Cubic spline 29.89 dB                      Adaptive interpolation 30.49 dB

     High-resolution image                            Cubic spline 28.66 dB                      Adaptive interpolation 29.27 dB

Fig. 3.   From left to right: high-resolution image, cubic-spline interpolation, adaptive directional interpolation with mixing weights.