image processing and Analysis
This page intentionally left blank
image
processing
and Analysis
variational PDE, ...

This content was uploaded by our users and we assume good faith they have the permission to share this book. If you own the copyright to this book and it is wrongfully on our website, we offer a simple DMCA procedure to remove your content from our site. Start by pressing the button below!

image processing and Analysis

This page intentionally left blank

image

processing

and Analysis

variational PDE, wavelet, and stochastic Methods

Tony F. chan

university OF california, LOS Anggles LOS AnGeLeS, california

jianHonG [jacKie] shen

university OF Minnesota MinneapoLis, Minnesota

Society for Industrial and Applied Mathematics Philadelphia

Copyright © 2005 by the Society for Industrial and Applied Mathematics. 1098765432 1 All rights reserved. Printed in the United States of America. No part of this book may be reproduced, stored, or transmitted in any manner without the written permission of the publisher. For information, write to the Society for Industrial and Applied Mathematics, 3600 University City Science Center, Philadelphia, PA 19104-2688. MATLAB® is a trademark of The MathWorks, Inc. and is used with permission. The MathWorks does not warrant the accuracy of the text or exercises in this book. This book's use or discussion of MATLAB® software or related products does not constitute endorsement or sponsorship by The MathWorks of a particular pedagogical approach or particular use of the MATLAB® software. For MATLAB® product information, please contact: The MathWorks, Inc., 3 Apple Hill Drive, Natick, MA 01760-2098 USA, 508-647-7000, Fax: 508-647-7101, [email protected], www.mathworks.com

Library of Congress Cataloging-in-Publication Data Chan, Tony F. Image processing and analysis : variational, PDE, wavelet, and stochastic methods / Tony F. Chan, Jianhong (Jackie) Shen. p. cm. Includes bibliographical references and index. ISBN 0-89871-589-X (pbk.) 1. Image processing—Mathematical models. I. Shen, Jianhong, 1971- II. Title. TA1637.C4775 2005 621.36'7-dc22

2005048960

• 51£LJTL is a registered trademark.

T0 Monica, Micfiettej Ityan, andCtaudia

This page intentionally left blank

Contents List of Figures

Xlll

Preface

xix

1

Introduction 1.1 Dawning of the Era of Imaging Sciences 1.1.1 Image Acquisition 1.1.2 Image Processing 1.1.3 Image Interpretation and Visual Intelligence 1.2 Image Processing by Examples 1.2.1 Image Contrast Enhancement 1.2.2 Image Denoising 1.2.3 Image Deblurring 1.2.4 Image Inpainting 1.2.5 Image Segmentation 1.3 An Overview of Methodologies in Image Processing .3.1 Morphological Approach .3.2 Fourier and Spectral Analysis .3.3 Wavelet and Space-Scale Analysis .3.4 Stochastic Modeling .3.5 Variational Methods .3.6 Partial Differential Equations (PDEs) .3.7 Different Approaches Are Intrinsically Interconnected 1.4 Organization of the Book 1.5 How to Read the Book

1 1 1 5 6 6 6 8 9 9 11 12 12 14 15 16 17 19 21 24 26

2

Some Modern Image Analysis Tools 2.1 Geometry of Curves and Surfaces 2.1.1 Geometry of Curves 2.1.2 Geometry of Surfaces in Three Dimensions 2.1.3 Hausdorff Measures and Dimensions 2.2 Functions with Bounded Variations 2.2.1 Total Variation as a Radon Measure 2.2.2 Basic Properties of BV Functions

31 31 31 36 44 45 46 49

VII

viii

Contents

2.3

2.4

2.5

2.6

3

2.2.3 The Co-Area Formula Elements of Thermodynamics and Statistical Mechanics 2.3.1 Essentials of Thermodynamics 2.3.2 Entropy and Potentials 2.3.3 Statistical Mechanics of Ensembles Bayesian Statistical Inference 2.4.1 Image Processing or Visual Perception as Inference 2.4.2 Bayesian Inference: Bias Due to Prior Knowledge 2.4.3 Bayesian Method in Image Processing Linear and Nonlinear Filtering and Diffusion 2.5.1 Point Spreading and Markov Transition 2.5.2 Linear Filtering and Diffusion 2.5.3 Nonlinear Filtering and Diffusion Wavelets and Multiresolution Analysis 2.6.1 Quest for New Image Analysis Tools 2.6.2 Early Edge Theory and Marr's Wavelets 2.6.3 Windowed Frequency Analysis and Gabor Wavelets 2.6.4 Frequency-Window Coupling: Malvar-Wilson Wavelets 2.6.5 The Framework of Multiresolution Analysis (MRA) 2.6.6 Fast Image Analysis and Synthesis via Filter Banks

Image Modeling and Representation 3.1 Modeling and Representation: What, Why, and How 3.2 Deterministic Image Models 3.2.1 Images as Distributions (Generalized Functions) 3.2.2 Lp Images 3.2.3 Sobolev Images H n (W) 3.2.4 BV Images 3.3 Wavelets and Multiscale Representation 3.3.1 Construction of 2-D Wavelets 3.3.2 Wavelet Responses to Typical Image Features 3.3.3 Besov Images and Sparse Wavelet Representation 3.4 Lattice and Random Field Representation 3.4.1 Natural Images of Mother Nature 3.4.2 Images as Ensembles and Distributions 3.4.3 Images as Gibbs'Ensembles 3.4.4 Images as Markov Random Fields 3.4.5 Visual Filters and Filter Banks 3.4.6 Entropy-Based Learning of Image Patterns 3.5 Level-Set Representation 3.5.1 Classical Level Sets 3.5.2 Cumulative Level Sets 3.5.3 Level-Set Synthesis 3.5.4 An Example: Level Sets of Piecewise Constant Images 3.5.5 High Order Regularity of Level Sets 3.5.6 Statistics of Level Sets of Natural Images

52 54 55 56 58 61 61 62 64 65 65 67 70 73 73 75 76 77 80 86 91 91 93 93 96 98 98 99 99 104 107 115 115 116 117 119 122 124 126 127 127 129 129 130 131

Contents

ix

3.6

The Mumford-Shah Free Boundary Image Model 132 3.6.1 Piecewise Constant 1-D Images: Analysis and Synthesis 132 3.6.2 Piecewise Smooth 1-D Images: First Order Representation . . . 134 3.6.3 Piecewise Smooth 1-D Images: Poisson Representation 135 3.6.4 Piecewise Smooth 2-D Images 136 3.6.5 The Mumford-Shah Model 138 3.6.6 The Role of Special BV Images 140

4

Image Denoising 4.1 Noise: Origins, Physics, and Models 4.1.1 Origins and Physics of Noise 4.1.2 A Brief Overview of 1-D Stochastic Signals 4.1.3 Stochastic Models of Noises 4.1.4 Analog White Noises as Random Generalized Functions 4.1.5 Random Signals from Stochastic Differential Equations 4.1.6 2-D Stochastic Spatial Signals: Random Fields 4.2 Linear Denoising: Lowpass Filtering 4.2.1 Signal vs. Noise 4.2.2 Denoising via Linear Filters and Diffusion 4.3 Data-Driven Optimal Filtering: Wiener Filters 4.4 Wavelet Shrinkage Denoising 4.4.1 Shrinkage: Quasi-statistical Estimation of Singletons 4.4.2 Shrinkage: Variational Estimation of Singletons 4.4.3 Denoising via Shrinking Noisy Wavelet Components 4.4.4 Variational Denoising of Noisy Besov Images 4.5 Variational Denoising Based on BV Image Model 4.5.1 TV, Robust Statistics, and Median 4.5.2 The Role of TV and BV Image Model 4.5.3 Biased Iterated Median Filtering 4.5.4 Rudin, Osher, and Fatemi's TV Denoising Model 4.5.5 Computational Approaches to TV Denoising 4.5.6 Duality for the TV Denoising Model 4.5.7 Solution Structures of the TV Denoising Model 4.6 Denoising via Nonlinear Diffusion and Scale-Space Theory 4.6.1 Perona and Malik's Nonlinear Diffusion Model 4.6.2 Axiomatic Scale-Space Theory 4.7 Denoising Salt-and-Pepper Noise 4.8 Multichannel TV Denoising 4.8.1 Variational TV Denoising of Multichannel Images 4.8.2 Three Versions of TV[w]

145 145 145 147 150 151 153 155 156 156 157 159 160 160 163 165 171 174 174 175 175 177 178 183 185 191 191 194 198 203 203 204

5

Image Deblurring 5.1 Blur: Physical Origins and Mathematical Models 5.1.1 Physical Origins 5.1.2 Mathematical Models of Blurs 5.1.3 Linear vs. Nonlinear Blurs

207 207 207 208 214

x

6

Contents

5.2 III-posedness and Regularization 5.3 Deblurring with Wiener Filters 5.3.1 Intuition on Filter-Based Deblurring 5.3.2 Wiener Filtering 5.4 Deblurring of BV Images with Known PSF 5.4.1 The Variational Model 5.4.2 Existence and Uniqueness 5.4.3 Computation 5.5 Variational Blind Deblurring with Unknown PSF 5.5.1 Parametric Blind Deblurring 5.5.2 Parametric-Field-Based Blind Deblurring 5.5.3 Nonparametric Blind Deblurring

216 217 217 218 220 220 222 224 226 226 230 233

Image Inpainting 6.1 A Brief Review on Classical Interpolation Schemes 6.1.1 Polynomial Interpolation 6.1.2 Trigonometric Polynomial Interpolation 6.1.3 Spline Interpolation 6.1.4 Shannon's Sampling Theorem 6.1.5 Radial Basis Functions and Thin-Plate Splines 6.2 Challenges and Guidelines for 2-D Image Inpainting 6.2.1 Main Challenges for Image Inpainting 6.2.2 General Guidelines for Image Inpainting 6.3 Inpainting of Sobolev Images: Green's Formulae 6.4 Geometric Modeling of Curves and Images 6.4.1 Geometric Curve Models 6.4.2 2-, 3-Point Accumulative Energies, Length, and Curvature . . . . 6.4.3 Image Models via Functionalizing Curve Models 6.4.4 Image Models with Embedded Edge Models 6.5 Inpainting BV Images (via the TV Radon Measure) 6.5.1 Formulation of the TV Inpainting Model 6.5.2 Justification of TV Inpainting by Visual Perception 6.5.3 Computation of TV Inpainting 6.5.4 Digital Zooming Based on TV Inpainting 6.5.5 Edge-Based Image Coding via Inpainting 6.5.6 More Examples and Applications of TV Inpainting 6.6 Error Analysis for Image Inpainting 6.7 Inpainting Piecewise Smooth Images via Mumford and Shah 6.8 Image Inpainting via Euler's Elasticas and Curvatures 6.8.1 Inpainting Based on the Elastica Image Model 6.8.2 Inpainting via Mumford-Shah-Euler Image Model 6.9 Inpainting of Meyer's Texture 6.10 Image Inpainting with Missing Wavelet Coefficients 6.11 PDE Inpainting: Transport, Diffusion, and Navier-Stokes 6.11.1 Second Order Interpolation Models 6.11.2 A Third Order PDE Inpainting Model and Navier-Stokes . . . .

245 246 246 248 249 251 253 256 256 257 258 263 264 265 268 269 270 270 272 274 274 276 277 279 282 284 284 287 289 291 295 295 299

Contents 6.11.3 TV Inpainting Revisited: Anisotropic Diffusion 6.11.4 CDD Inpainting: Curvature Driven Diffusion 6.11.5 A Quasi-axiomatic Approach to Third Order Inpainting 6.12 Inpainting of Gibbs/Markov Random Fields 7

xi 301 302 303 307

Image Segmentation 309 7.1 Synthetic Images: Monoids of Occlusive Preimages 309 7.1.1 Introduction and Motivation 309 7.1.2 Monoids of Occlusive Preimages 310 7.1.3 Mimimal and Prime (or Atomic) Generators 315 7.2 Edges and Active Contours 318 7.2.1 Pixelwise Characterization of Edges: David Marr's Edges . . . .318 7.2.2 Edge-Regulated Data Models for Image Gray Values 320 7.2.3 Geometry-Regulated Prior Models for Edges 322 7.2.4 Active Contours: Combining Both Prior and Data Models . . . . 325 7.2.5 Curve Evolutions via Gradient Descent 327 7.2.6 F-Convergence Approximation of Active Contours 329 7.2.7 Region-Based Active Contours Driven by Gradients 331 7.2.8 Region-Based Active Contours Driven by Stochastic Features 332 7.3 Geman and Geman's Intensity-Edge Mixture Model 338 7.3.1 Topological Pixel Domains, Graphs, and Cliques 338 7.3.2 Edges as Hidden Markov Random Fields 339 7.3.3 Intensities as Edge-Regulated Markov Random Fields 342 7.3.4 Gibbs' Fields for Joint Bayesian Estimation of u and F 343 7.4 The Mumford-Shah Free-Boundary Segmentation Model 344 7.4.1 The Mumford-Shah Segmentation Model 344 7.4.2 Asymptotic M.-S. Model I: Sobolev Smoothing 345 7.4.3 Asymptotic M.-S. Model II: Piecewise Constant 347 7.4.4 Asymptotic M.-S. Model III: Geodesic Active Contours 351 7.4.5 Nonuniqueness of M.-S. Segmentation: A 1-D Example 355 7.4.6 Existence of M.-S. Segmentation 355 7.4.7 How to Segment Sierpinski Islands 359 7.4.8 Hidden Symmetries of M.-S. Segmentation 362 7.4.9 Computational Method I: F-Convergence Approximation . . . .364 7.4.10 Computational Method II: Level-Set Method 366 7.5 Multichannel Logical Segmentation 369

Bibliography

373

Index

393

This page intentionally left blank

List of Figures 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

1.9 1.10 1.11

1.12 1.13 1.14 1.15 1.16

An ideal image: noiseless, complete, and in good contrast A low-contrast version of the ideal image Degradation by additive Gaussian noise (Chapter 4) Salt-and-pepper noise with 10% spatial density (Chapter 4) Degradation by out-of-focus blur: the digital camera focuses on a fingertip about 1 inch away while the target scene is about 1 foot away (Chapter 5) Degradation by motion blur due to horizontal hand jittering during a single exposure (Chapter 5) 150 8-by-8 packets are randomly lost during transmission (Chapter 6). The goal of error concealment (or more generally, inpainting) is to develop models and algorithms that can automatically fill in the blanks [24, 67, 166]. Such cartoonish segmentation seems trivial to human vision but still remains the most fundamental and challenging problem in image processing and low-level computer vision (Chapter 7). (The background segment Wo is not shown explicitly here.) A binary set A and a structure element S (for dilation and erosion). ... Dilation D s ( A ) (left) and erosion E s ( A ) (right): dilation closes up small holes or gaps, while erosion opens them up A digital image and its discrete Fourier transform. The example reveals a couple of salient features of Fourier image analysis: (1) most highamplitude coefficients concentrate on the low-frequency band; (2) dominant directional information in the original image is easily recognizable in the Fourier domain; and (3) the coefficients, however, decay slowly for Heaviside-type directional edges (i.e., a jump line with distinct constant values along its two shoulders) An example of a (mother) wavelet by Daubechies' design [96]. Localization and oscillation are characteristic to all wavelets A noisy 1-D signal and its optimal restoration according to (1.7) A trefoil evolves under mean-curvature motion (1.11) or (1.12) Different approaches are intrinsically connected: an example Organization of the book

XIII

7 7 8 9 10

10 11

11 13 14

15 16 19 21 22 24

xiv

List of Figures

2.1 2.2 2.3 2.4

2.5 2.6 2.7

2.8

2.9

2.10

2.11 2.12 2.13

2.14 2.15 2.16 2.17

A planar curve: tangent t, normal n, and curvature K Second order local geometry: surface normal N and two perpendicular principle curvature circles A surface patch z — h(u, v] — cos(u)cos(u) (top), its mean curvature field H (lower left), and Gauss curvature field K (lower right) A set E (left) and a covering A (right). Covering elements like the huge box with dashed borders (right) fail to faithfully capture small-scale details. The definition of Hausdorff measures forces covering scales to tend to zero The Hausdorff dimension dimH(E) of a set E is the critical d, above which E appears too "thin," while below it appears too "fat." Three 1-D images with TV[f] = TV[g] = TV[h] = 2 An example of L 1 -lower semicontinuity. The sequence (un) of 1-D images on [0, 1 ] converges to u = O i n L ' since ||u,n+1 — u||L1 2-n. Notice that TV(u) = 0 while TV(u n ) = 2, which is consistent with the property of lower semicontinuity: TV(u) lim inf n TV(u n ). In particular, strict inequality is indeed realizable Geometric meaning of TV: the co-area formula. For smooth images, TV[w] is to sum up the lengths of all the level curves, weighted by the Lebesgue element dk. Plotted here is a discrete approximation: TV[«] ~ £„ length(« = X,,) AA Gibbs' CE: higher temperature T corresponds to smaller /3 and more uniform distribution; on the other hand, when T — 0, the system exclusively remains at the ground state (leading to superfluids or superconductance in physics) Gibbs' entropies for two imaginary 4-state systems (with K set to 1). Entropy generally measures the degrees of freedom of a system. A lower entropy therefore means that the target system is more restrained. (This observation led Shannon [272] to define negative entropies as information metrics, since less randomness implies more information.) An example of linear anisotropic diffusion by (2.38) with diagonal diffusivity matrix D = diag(D v , D v ) and Dy : D v = 10 : 1. The image thus diffuses much faster along the vertical y-direction The median of a 5-component segment from an imaginary 1-D signal (*„). An example of using the median filter (with 7 x 7 centered square window) to denoise an image with severe salt-and-pepper noise. Notice the outstanding feature of median filtering: the edges in the restored image are not blurry Two examples of Marr's wavelets ("Mexican hats") as in (2.50) A generic approach to designing the window template w(x) via symmetrically constructing the profile of its square w2(x) MRA as (Hilbert) space decompositions: the finer resolution space V-± is decomposed to the detail (or wavelet) space W\ and coarser resolution space V\. The same process applies to V\ and all other V/s [290] A pair of compactly supported scaling function 4>(x) and mother wavelet iAU) by Daubechies' design [96]

32 38 40

45 45 47

49

54

59

60 68 70

71 75 79 84 87

List of Figures 2.18 3.1

3.2 3.3

3.4

3.5 3.6

3.7

4.1 4.2

4.3 4.4

4.5

Fast wavelet transform via filter banks: the two-channel analysis (or decomposition) and synthesis (or reconstruction) banks

xv

88

Images as distributions or generalized functions. Test functions then model various biological or digital sensors, such as retinal photoreceptors of human vision or coupled charge devices in CCD cameras 94 Three Haar mother wavelets in two dimensions via tensor products. . . . 103 Besov norms in B®(LpYs measure the strength of signals in the spacescale plane: Lp for intrascale variations, Lq for inter- or cross-scale variations (in terms of dX = dh/h), while h~a for comparison with Holder continuity 108 Ingredients of Markov random fields by examples: a neighborhood Na (left), two doubleton cliques C e C (middle), and locality of conditional inference p(ua \ «n\ a ) = p(ua \ uNu) (right) 120 An example of branches LX,AA.'S with A A, = h and A. = nh. The branch L^,h contains two leaflets 132 The(xn, bn, gn) -representation of a (compactly supported) piecewise smooth signal u. On smooth regions where the signal u varies slowly, gn is often small and only a few bits suffice to code them 135 Poisson representation of a piecewise smooth signal u. Shown here is only a single representative interval. The signal u is represented by the two boundary values «+ and u ~+, and its second order derivative / = u" (corresponding to the source distribution in electromagnetism). Reconstructing u on the interval then amounts to solving the Poisson equation (3.62). The advantage of such representation is that / is often small (just like wavelet coefficients) for smooth signals and demands fewer bits compared with u 136 A sample 2-D Brownian path W(t) with 0 < / < 4, and the reference circle with radius 2 = V4, the standard deviation of W (4) A sample random signal generated by the SDE dX = 2Xdt + XdW with X(Q) = 1. The smooth dashed curve denotes the mean curve x(t) = E X ( t ) = e2'. The random signal is clearly not stationary since both the means and variances of X ( t ) evolve Hard thresholding T\(t) and soft shrinkage 5^(0 The shrinkage operator a* = Sff(a0) achieves minimum shrinkage (Theorem 4.4) among all the estimators a that satisfy the uniform shinkage condition (4.25) Two singleton error functions e\(t\ao) with A. = JJL = 1 and ao = 1,2. Notice that the optimal estimators (i.e., the valleys) a = 0, 1 are precisely the shrinkages Sa(ao) with a = fji/'k = 1

152

155 161

163

165

xvi

List of Figures

4.6

4.7 4.8

4.9 4.10 4.11

4.12 4.13

5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9

6.1 6.2

6.3

An example of shrinkage-based image denoising, with Gaussian white noise level a =0.1, and four levels of wavelet decompositions based on one of Daubechies' orthonormal wavelets [96]. The three bottom panels demonstrate the effect of the threshold parameter A on the shrinkage operator Sx '• a too large X causes overshrinkage and oversmoothing, while a too small one leads to the opposite 171 Same example as in Figure 4.6 with increased Gaussian noise level a =0.2 172 Denoising effects of moving means vs. moving medians. Both employ the same symmetric (moving) window of four neighbors on either side. As expected, the median filter preserves sharp edges better than the mean filter. (The dashed curve denotes the ideal clean signal.) 175 An example of Rudin, Osher, and Fatemi's TV denoising 179 A target pixel O and its neighbors 181 An example of applying Perona and Malik's anisotropic diffusion to a noisy image. The nonlinear diffusivity coefficient D(\ VM|) for this example is the one in (4.93) and (4.94). Notice the remarkable feature of edge preservation, as contrast to linear heat diffusions that invariably blur sharp edges while removing noises 194 Two examples of salt-and-pepper noises and the denoising performance of median filtering. As consistent with theoretical analysis, median filtering works very efficiently for low spatial densities but poorly for high ones. 202 An example of applying the TV inpainting technique (4.103) (also see Chan and Shen [67] or Chapter 6) to the denoising of high-density saltand-pepper noises 203 An example of motion blur (due to camera jittering) Motion blur of the image of a point source An example of out-of-focus blur Geometric optics of out-of-focus imaging An example of Wiener filtering for both denoising and deblurring Deblurring an out-of-focus image (with known PSF) Deblurring an image blurred by horizontal hand jittering Another example of restoring an image blurred by motion A computational example for the double-BV blind deblurring model (by Chan and Wong [77]). Left: the blurry image; Right: the deblurred image The goal of inpainting is to reconstruct the ideal image u on the entire domain based on the incomplete and often degraded data u° available outside the missing (or inpainting) domain D The effect of local and global pattern recognition on image inpainting: the black color seems to be a reasonable inpainting solution in the left panel, while for the right chessboard pattern, the white color becomes more plausible (Chan and Shen [67]) The effect of aspect ratios on image inpainting [67]

209 210 212 213 219 225 226 226 243

246

257 258

List of Figures

6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17 6.18 6.19 6.20 6.21 7.1

7.2 7.3 7.4 7.5 7.6

xvii

Harmonic inpainting of a smooth image (u = r = y/x 2 + j2) and an ideal step edge 263 Can the TV inpainting model explain Kanizsa's Entangled Man? 272 The model for Kanizsa's Entangled Man 273 Inpainting a noisy edge 277 Inpainting occluded bars 277 Inpainting a noisy face 278 Inpainting for text removal 278 Digital harmonic zoom-in and TV zoom-in (Chan and Shen [67]) 279 Edge decoding by TV inpainting (Chan and Shen [67]; also see the recent pattern-theoretic approach by Guo, Zhu, and Wu [148]) 280 The aspect ratios of inpainting domains influence the performance of inpainting schemes (Chan and Kang [59]) 281 Inpainting based on the F-convergence approximation (6.42) and its associated elliptic system (6.47) 284 Text erasing by inpainting based on the Mumford-Shah image model. . .284 Two examples of elastica inpainting, as compared with TV inpainting. In the case of large aspect ratios [68], the TV inpainting model fails to comply to the connectivity principle 286 Inpainting based on the Mumford-Shah-Euler image model can satisfactorily restore a smooth edge as expected 289 An example of inpainting a noisy and incomplete BV image with missing wavelet components using the model (6.64) (Chan, Shen, and Zhou [74]). 294 Another example of inpainting a noisy and incomplete BV image with missing wavelet components using the model (6.64) (Chan, Shen, and Zhou [74]) 294 TV fails to realize the connectivity principle in inpainting problems with large scales (or aspect ratios) 301 An example of CDD inpainting for scratch removal (Chan and Shen [68]) 302 Role of occlusion in visual perception. Left panel: a gray vase against a bright background or two opposite human faces against a gray background? This is a typical illusion problem caused by the lack of occlusion cues; Right panel (a trefoil): the role of occlusion in knot theory [253]. An example of two preimages aY and b$ and their occlusion ay H bs. Complex image patterns in the real world often originate from simple objects via occlusion Left: a diskette preimage; Right: a planar preimage On an edge pixel x (with y ( x ) — 1), p — \ Vu\(x) is more likely to be large according to (7.17) An example of geodesic active contouring: the given image and three different stages of the contour evolution by (7.32) An example of Chan and Vese's Active Contour model (7.52) without gradients (i.e., region based and driven by mean fields): three different evolution stages for a sectional brain image (from Chan and Vese [75]).

.310 312 316 322 329 . 337

xviii 7.7

7.8 7.9

7.10

7.11 7.12

7.13

7.14 7.15 7.16 7.17

List of Figures An example of texture segmentation by combining Chan and Vese's Active Contour model (7.52) and Gabor filters (7.53) (shown at increasing times; see [65, 260]) A pixel lattice £1 = [a,b, ...} and an edge lattice E = [a, @ , . . . } Six binary edgel configurations on a maximal edge clique and their empirical potentials V^(F)'s assigned by Geman and Geman [130] (where A is a constant) The total curvatures (in the interior of the pixel squares) of closed tours around the six edge configurations qualitatively match Geman and Geman's empirical choice of potentials in Figure 7.9 Visualization of the example: the set S consisting of the circles at all scale levels has the weird property H* (S) < oo but ~Hl (S) = oo, since S = Q. Left panel: the complement G of the classical Sierpinski gasket (Mandelbrot [206]) consists of all the shaded open triangles (excluding their sides). Right panel: Sierpinski Islands G(p) is the open set obtained from G by contracting each triangle of G with some suitable rate controlled by some p e (0, 3/2) (see the text). And the target image up discussed in this section is the indicator function of G(p) An example of F-convergence-based approximation to the M.-S. segmentation model (7.97): the optimal image estimation u = u £ ( x ) and its associated edge "canyon" function z — z£(x). Computation has been based on the alternating minimization scheme (7.98) An example of M.-S. segmentation using the level-set algorithm by Chan and Vese [75, 76] A synthetic example of an object in two different channels. Notice that the lower left corner of A\ and the upper corner of A 2 are missing Different logical combinations for the sample image: the union, the intersection, and the differentiation Region-based logical model on a medical image. In the first channel A \, the noisy image has a "brain tumor," while channel AI does not. The goal is to spot the tumor that is in channel A\, but not in AT, i.e., the differentiation A\ D —>A2. In the right column, we observe that the tumor has been successfully captured

338 340

341

342 357

360

366 369 369 370

372

Preface No time in human history has ever witnessed such explosive influence and impact of image processing on modern society, sciences, and technologies. From nanotechnologies, astronomy, medicine, vision psychology, remote sensoring, security screening, and the entertainment industry to the digital communication technologies, images have helped mankind to see objects in various environments and scales, to sense and communicate distinct spatial or temporal patterns of the physical world, as well as to make optimal decisions and take right actions. Image processing and understanding are therefore turning into a critical component in contemporary sciences and technologies with many important applications. As a branch of signal processing, image processing has traditionally been built upon the machinery of Fourier and spectral analysis. In the past few decades, there have emerged numerous novel competing methods and tools for successful image processing. They include, for example, stochastic approaches based upon Gibbs/Markov random fields and Bayesian inference theory, variational methods incorporating various geometric regularities, linear or nonlinear partial differential equations, as well as applied harmonic analysis centered around wavelets. These diversified approaches are apparently distinct but in fact intrinsically connected. Each method excels from certain interesting angles or levels of approximations but is also inevitably subject to its limitations and applicabilities. On the other hand, at some deeper levels, they share common grounds and roots, from which more efficient hybrid tools or methods can be developed. This highlights the necessity of integrating this diversity of approaches. The present book takes a concerted step towards this integration goal by synergistically covering all the aforementioned modern image processing approaches. We strive to reveal the few key common threads connecting all these major methodologies in contemporary image processing and analysis, as well as to highlight some emergent integration efforts that have been proven very successful and enlightening. However, we emphasize that we have made no attempt to be comprehensive in covering each subarea. In addition to the efforts of organizing the vast contemporary literature into a coherent logical structure, the present book also provides some in-depth analysis for those relatively newer areas. Since very few books have attempted this integrative approach, we hope ours will fill a need in the field. Let u denote an observed image function, and T an image processor, which can be either deterministic or stochastic, as well as linear or nonlinear. Then a typical image

XIX

xx

Preface

processing problem can be expressed by the flow chart

where F represents important image or visual features of interest. In the present book, we explore all three key aspects of image processing and analysis. • Modeling: What are the suitable mathematical models for u and Tl What are the fundamental principles governing the constructions of such models? What are the key features that have to be properly respected and incorporated? • Model Analysis: Are the two models for u and T compatible? Is T stable and robust to noises or general perturbations? Does F = T[u] exist, and if so, is it unique? What are the fine properties or structures of the solutions? In many applications, image processors are often formulated as inverse problem solvers, and as a result, issues like stability, existence, and uniqueness become very important. • Computation and Simulation: How can the models be efficiently computed or simulated? Which numerical regularization techniques should be introduced to ensure stability and convergence? And how should the targeted entities be properly represented? This view governs the structure and organization of the entire book. The first chapter briefly summarizes the emerging novel field of imaging science, as well as outlines the main tasks and topics of the book. In the next two chapters, we introduce and analyze several universal modern ways for image modeling and representation (for u), which include wavelets, random fields, level sets, etc. Based on this foundation, we then in the subsequent four chapters develop and analyze four specific and significant processing models (for T) including image denoising, image deblurring, inpainting or image interpolation, and image segmentation. Embedded within various image processing models are their computational algorithms, numerical examples, or typical applications. As the whole spectra of image processing spread so vastly, in this book we can only focus on several most representative problems which emerge frequently from applications. In terms of computer vision and artificial intelligence, these are often loosely categorized as low-level vision problems. We do not intend to cover high-level vision problems which often involve pattern learning, identification, and representation. We are enormously grateful to Linda Thiel, Alexa Epstein, Kathleen LeBlanc, Michelle Montgomery, David Riegelhaupt, and Sara Murphy of the SI AM Publisher for their constant encouragement and care throughout the project. It has been such a wonderful experience of planning, communication, and envisaging. We also owe profound gratitude to the following colleagues whose published works and personal discussions have greatly influenced and shaped the contents and structures of the current book (in alphabetical order): Antonin Chambolle, Ron Coifman, Ingrid Daubechies, Rachid Deriche, Ron DeVore, David Donoho, Stu Geman, Brad Lucier, Jitendra Malik, Yves Meyer, Jean-Michel Morel, David Mumford, Stan Osher, Pietro Perona, Guillermo Sapiro, Jayant Shah, James Sethian, Harry Shum, Steve Smale, Gilbert Strang, Curt Vogel, Yingnian Wu, Alan Yuille, and Song-Chun Zhu, and the list further expands. The book would be impossible without the generous support and help of many friends: Doug Arnold, Andrea Bertozzi, Carme Calderon, Charles Chui, Babette Dalton, Bjorn

Preface

xxi

Enquist, Bob Gulliver, Xiaoming Huo, Kate Houser, Yoon-Mo Jung, Dan Kersten, Paul Schrater, Mitch Luskin, Riccardo March, Willard Miller, Peter Olver, Hans Othmer, Fadil Santosa, Zuowei Shen, Chi-Wang Shu, Luminita Vese, Kevin Vixie, Tony Yezzi, Hong-Kai Zhao, and Ding-Xuan Zhou. Tony Chan would like to personally thank his colleagues for numerous inspirations and interesting discussions: Emmanuel Candes, Raymond Chan, Li-Tien Cheng, Ron Fedkiw, Mark Green, Michael Ng, Stefano Soatto, Xue-Cheng Tai, Richard Tsai, and Wing Wong. Tony Chan is also deeply grateful to the following students and postdoctoral fellows for the privilege of working with them in this new and exciting field: Peter Blomgren, Jamylle Carter, Selim Esedoglu, Sung-Ha Kang, Mark Moelich, Pep Mulct, Fred Park, Berta Sandberg, Jackie Shen, Bing Song, David Strong, Justin Wan, Yalin Wang, Luminita Vese, Chiu-Kwong Wong, Andy Yip, Hao-Min Zhou, and Wei Zhu, Jackie Shen wishes to personally thank Professors Gilbert Strang, Tony Chan, Stan Osher, David Mumford, and Stu Geman for their profound influence on his scholastic growth in the field, as well as numerous personal friends for warming and shining up each ordinary day during the project: Tianxi Cai, Shanhui Fan, Chuan He, Huiqiang Jiang, Ming Li, Tian-Jun Li, William Li, Chun Liu, Hailiang Liu, Huazhang (Andy) Luo, Conan Leung, Mila Nikolova, Jiaping Wang, Chao Xu, Jianling Yuan, Wen Zhang, and Qiang Zhu. We are very grateful to Jean-Michel Morel and Kevin Vixie specifically for their insightful and constructive comments on an early version of the manuscript for this book, which have significantly improved the quality. During this book project, the authors are enormously grateful for the support from the National Science Foundations (NSF-USA), the Office of Naval Research (ONR-USA), as well as the National Institute of Health (NIH-USA). In particular, Tony Chan would like to thank Wen Master (at ONR) for the continuous support to this novel area in applied and computational mathematics. We also acknowledge the tremendous benefits from the participation of numerous imaging sciences related workshops at the Institute of Pure and Applied Mathematics (IPAM) at UCLA, the Institute of Mathematics and Its Applications (IMA) at the University of Minnesota, the Institute of Mathematical Sciences (IMS) at the National University of Singapore, the Mathematical Sciences Research Institute (MSRI) at Berkeley, as well as the Center of Mathematical Sciences (CMS) at Zhejiang University, China. Finally, this book project, like all the others in our life, is an intellectual product under numerous constraints, including our busy working schedules and many other scholastic duties. Its contents and structures as presented herein are therefore only optimal subject to such inevitable conditions. All errata and suggestions for improvements will be received gratefully by the authors. This book is absolutely impossible without the pioneering works of numerous insightful mathematicians and computer scientists and engineers. It would be our great pleasure to see that the book can faithfully reflect many major aspects of contemporary image analysis and processing. But unintentional biases are inevitable due to the limited views and experiences of the authors, and we are happy to hear any criticisms from our dear readers.

This page intentionally left blank

Chapter 1

Introduction

The last few decades have witnessed a new era of imaging sciences. From satellite imaging and X-ray imaging to modern computer aided CT, MRI, and PET imaging in medical sciences, mankind's naked vision is no longer restricted by the conventional notions of space, time, scales, and visibility. Consequently, just as Mother Nature has achieved for biological vision, there is tremendous need for developing the virtual brain structures for all these imaging areas, i.e., models and algorithms for making these mechanical or artificial vision systems more efficient, accurate, and stable, as well as artificial intelligence that can successfully process, communicate, and interpret the resultant images. The current book attempts to offer a well-organized view of this rapidly evolving field. This introductory chapter first gives an overview of this exciting emerging field, and then introduces the scopes and structures of the current book. A detailed guidance on how to efficiently read and use the book is also speculated for different groups of potential readers.

1.1

Dawning of the Era of Imaging Sciences

Imaging sciences consist of three relatively independent components: image acquisition, image processing, and image interpretation. This book mainly focuses on the second subject, though we shall first give a quick overview of all three.

1.1.1

Image Acquisition

Image acquisition (orformation) studies the physical mechanisms by which particular types of imaging devices generate image observations, and it also investigates the associated mathematical models and algorithms employed by the computers integrated into such imaging devices. We now briefly describe several important image acquisition problems. 1

2

Chapter 1. Introduction

Human Vision The human vision system is Mother Nature's biological device for carrying out intelligent visual perception. The imaging formation system mainly consists of optical peripherals including the cornea, pupil, iris, and lens, as well as two different types of photoreceptors of the retinas, cones and rods, which are responsible for capturing colors and gray value intensities separately [168, 227]. The left and right retinas make up the frontal end of a complex vision system that transforms light or photons into electrochemical signals and prepares for further image processing and interpretation in the inner visual cortices through optical neuron fibers. The perceptible light spectra for average naked vision are about between 400 and 700nm (nanometer, or 10~9 meter) in terms of wavelengths. On the other hand, human vision achieves astounding sensitivity to a large dynamic range of light intensities due to its delicate adaptive mechanism described by the celebrated Weber's law in vision psychology, psychophysics, and retinal physiology [121, 168, 275, 315]. Night Vision and Thermal Imaging Wavelengths longer than 700nm belong to the infrared (IR) range. Mother Nature forces human beings to sleep and replenish during the night by making IR light imperceptible to normal naked eyes. However, assisted by proper signal amplification and enhancement devices (e.g., night vision goggles), human beings are capable of perceiving IR images ranging from just above 700nm to a few microns (10~6 meter). Furthermore, the majority of light waves in the IR spectra (ranging from 3 to 30 microns) are not reflected by objects as in the visible spectra but emitted by objects due to their intrinsic thermal and statistical activities at atomic levels. For such weak photons, thermal imaging devices can still collect, enhance, and process them, producing the socalled thermograms that map the detailed distribution of temperatures of a target scene. Assisted with proper signal processing hardware and software, one can then visualize the otherwise imperceptible images. Radar Imaging and SAR Radar is a device that employs RAdio waves to make Detection And Ranging of targets. The signals are in the microwave range between 1cm and 1m wavelengths and polarized. Radar antennas emit these polarized radio wave pulses and receive their echoes reflected from objects. The waiting time between sending and receiving, as well as the strength of the echoes or backscalters, reveal important information such as ranges, local configuration, and material differentiation of the targets. Applied to imaging, radar becomes a powerful tool for remote sensoring. To acquire ground images, a radar device is usually attached to a flight vehicle such as a space shuttle and continuously sends, receives, and processes signalling microwave pulses. The radar image is then constructed along its flight path. The resolution of the acquired image is inversely proportional to the size of the radar; that is, larger antennas result in higher resolution images. This hardware restriction on image resolution can, however, be effectively loosened using clever models, algorithms, and software that can extract extra spatial information

1.1. Dawning of the Era of Imaging Sciences

3

from the correlation among the echoes. This amounts to increasing the effective size of the real antenna being employed, and thus the resultant imaging technique is called synthetic aperture radar (SAR) imaging. Aperture refers to the opening area of an antenna that collects echoed pulses. Ultrasound Imaging Like ultraviolet light, ultrasound refers to acoustic waves whose frequencies are higher than the audible range of normal human ears, i.e., about between 20Hz and 20KHz. Medical ultrasound is often above 1MHz and below 20MHz. The different responses of the various tissues and structures to ultrasound pulses form the physiological foundation for ultrasound imaging. By measuring the time intervals between emitting and receiving sound pulses, as well as their strengths, the range and local topological information of target physiological entities can be inferred and displayed using computers. Furthermore, using Doppler shifts, ultrasound imaging can be employed to study organs in real-time motion (e.g., heart valves and blood vessels). Recall that the Doppler shift refers to the phenomenon that the reflected echoes from a moving object have higher or lower frequencies than the incident waves depending on whether the motion is toward or away from the wave source. Computed (Axial) Tomography (CAT or CT) Ever since the legendary discovery of X-rays by the German physicist Wilhelm Rontgen in 1895, who was then awarded the first Nobel Prize of physics in 1901, X-rays have become the most useful probing signals in medical imaging and diagnosis. Unlike the visible spectra, IR range for night vision or thermal imaging and the microwaves for radar, X-rays have wavelengths in the nanometer (10~ 9 ) or picometer (10~ 12 ) scales, which are energetic enough to penetrate most materials including biological tissues. This makes X-rays ideal for noninvasive medical probing and imaging. Different tissues and organs in the human body absorb X-rays at different rates. Denser structures like the bones absorb more X-rays than softer tissues. Thus conventional X-ray imaging is very similar to ordinary photographing, with a film directly exposed to the penetrated X-rays, and the different intensity patterns on the film differentiate different organs and structures. Computed (axial) tomography, often referred to as CAT or CT scanning, improves the above conventional X-ray imaging method by more delicate hardware design as well as computer software for reconstructing the final images. CT can be called rotational and localized X-ray imaging technology since it consists of a localized X-ray source and a detector which are put opposite to each other and rotate around the human body's longitudinal central axis. A two-dimensional (2-D) slice image is then synthesized from the collected X-ray signals along all the directions. Since the mid 1970s, CT scanning has witnessed broad applications in medical imaging and diagnosis, thanks to its fast speed and patient-friendly nature, as well as its unique capability in imaging a wide variety of soft tissues, hard bones, blood vessels, and so on.

4

Chapter 1. Introduction

Magnetic Resonance Imaging (MRI) The 2003 Nobel Prize in Physiology or Medicine was awarded to Paul Lauterbur and Peter Mansfield for their contributions on magnetic resonance imaging (MRI). The MRI imaging technology has been built upon the remarkable physical theory of nuclear magnetic resonance (NMR), discovered and developed by two American physicists, Felix Bloch and Edward Purcell, who won the Nobel Prize in 1952. A nucleon such as a proton or neutron possesses spin, which comes at a multiple of 1/2 with plus or minus differentiation. Paired opposite spins can effectively cancel each other out, leading to no evident observables. In an external magnetic field, unpaired or net spins tend to align with the orientation of the field, and statistical mechanics describes the distribution of two opposite alignments. The two alignments in different energy levels can transit to each other by either absorbing or releasing a photon energy proportional to the external magnetic field. For hydrogen atoms in clinical applications, for example, the frequencies of such photons are in the radio frequency (RF) range, i.e., MHz. Simply speaking, MRI imaging works as follows. First the atomic spins (e.g., hydrogen nuclei) of the human body are aligned by applying a strong external magnetic field. Then high frequency RF pulses are emitted into a slice plane perpendicular to the external field. By absorbing the photons in the RF waves, the aligned spins can transit from the lower energy alignment to the higher (or excited) one. After the RF waves are turned off, the excited nucleus ensemble starts to gradually resume to its original distribution, a process known as relaxation. It is the relaxation time that reveals or differentiates different tissues or body structures. MRI imaging hardware and software measure the distribution of relaxation times and construct the associated spatial images. MRI technology has many new important developments, including functional MRI (fMRI) for brain mapping [36] and diffusion tensor imaging (DTI) for studying neural fibers [314, 319]. Computer Graphics and Synthesized Images All the above image acquisition processes concern how to design hardware and combine it with powerful software to construct targeted images in the real three-dimensional (3-D) world. Computer graphics, on the other hand, is a purely so/rvvare-based image acquisition approach, or rather, an art that creates images of imaginary 3-D scenes. As a result, computer graphics has become increasingly important in the movie industry (e.g., popular movies like Shrek, Finding Nemo, Titanic, The Matrix, and The Lord of the Rings, to name only a few). Computer graphics can be formulated as follows. With an imaginary scene in mind, a designer first creates the 3-D coordinates G for shapes and geometric information of imaginary objects, as well as their optical properties such as the degree of transparency and surface reflectance R. The designer then designs a light source / for illuminance, which could be a point source or scattered sources such as in a normal daytime environment, and also picks a viewpoint 9 of an imaginary observer. The final imaginary image u is then synthesized based on some proper optical relation F:

1.1. Dawning of the Era of Imaging Sciences

5

If motion is involved, such as in movie production, time t has to be included as well:

which then becomes a dynamic process of image sequence generation. Human dreaming can be considered as Mother Nature's subconscious art of graphics, not by computers but by human visual neurons and cortices. In addition to all the aforementioned optical, electromagnetic, or nuclear imaging areas, there also exist other imaging techniques such as underwater or geophysical acoustic imaging (see, e.g., Papanicolaou [244]).

1.1.2

Image Processing

Any image processor can be abstractly formulated by an input-output system:

The input UQ denotes an already acquired image, which could be degraded due to either poor imaging conditions or problems during storage and communication. Mathematically speaking, an image processor T could be any linear or nonlinear operator that operates on the inputs and produces targeted features or patterns F. The design and development of image processors are often driven by specific applications or tasks. For example, in situations where observed images are degraded, one wishes to restore or enhance them to obtain high-quality images, based upon which other important visual decisions can be made more robustly and safely. Tumor detection from lowly contrasted CT images is such an example, while detecting celestial details such as binary stars from a blurry Hubble telescope image is another one. The main challenge for designing and computing T is that most image processing problems are ill-posed inverse problems. That is, there often exists a. forward problem that generates an image observation UQ from the targeted feature variable F if F is known or given. Most image processing problems are, however, to recover or detect F from UQ, and are therefore inverse problems. Ill-posedness often results from the nonuniqueness or instability of direct inversion. The output F = T[UQ] can be generally any collection of features that are visually meaningful. For example, F could denote the location of corners, the layout of object boundaries, the disjoint regions associated with distinct objects, or the relative order or depth of different objects in a scene. When the input UQ is a degraded image, the output F could also simply be its enhanced ideal version, often denoted by u, The input UQ can contain more than one image, such as in surveillance cameras, video or movie processing, and continuous medical monitoring of a single patient. For video or movie processing in particular, the time variable can lead to new features such as velocity distribution (or optical flows), dynamic evolution of the shapes and areas of interesting targets (or target tracking), and image registration. This book mainly focuses on image processing, and more specifically, on four classes of processing tasks which frequently arise from numerous applications: denoising, deblurring, inpainting or image interpolation, and segmentation. We shall further elaborate on them in the next major section.

Chapter 1. Introduction

6

By far it seems that image processing can be separated apart from image acquisition. This is, however, only partially true and is dependent on the real tasks at hand. A doctor who tries to decipher the image patterns of a patient certainly does not need to know nuclear physics or the physics details of MRI. But knowledge of an image acquisition process can often help develop more powerful image processing models or algorithms. One also must be aware that many image acquisition devices have their built-in image processors in order to construct high-quality images during the image formation stage.

1.1.3

Image Interpretation and Visual Intelligence

The ultimate purpose of imaging sciences is to be able to "see," monitor, and interpret the targeted portion of the world being imaged, whether it is the surface of Mars, the newly discovered tenth planet of the solar system, Sedna, or an MRI image slice of the brain. Therefore image interpretation is of fundamental importance for imaging sciences and is intrinsically connected to vision-based artificial intelligence. Classical image processing is often loosely identified with low-level vision, while image interpretation generally belongs to high-level vision. An ideal artificial vision system shall be able to perform all the major functions of the human vision system, including depth perception, motion estimation, and object differentiation and recognition. A combination of image processing and image interpretation is often necessary for these tasks. Image interpretation is partially the inverse problem of image acquisition. The latter studies how to form 2-D images from 3-D structures and patterns, which are arranged or positioned with relative orders, depths, opaqueness or transparency, etc. Image interpretation, on the other hand, attempts to reconstruct or interpret the 3-D world from 2-D images. For example, a well-trained radiologist can identify some abnormal image patterns in an MRI image with tumor tissues in a real 3-D body. Mathematically speaking, the foundation of image interpretation is a rapidly growing field called pattern theory, which mainly covers pattern modeling and representation, learning theory, and pattern detection and identification. This developing field in mathematics has been pioneered by, for example, S. Geman and D. Geman [130], Grenander [143], Mumford [224], Poggio and Smale [252], and Vapnik [307].

1.2

Image Processing by Examples

This section explains some typical image processing tasks through tangible examples. Some will be explored in later chapters with great details, while others can be found in more classical image processing books [140, 194]. Throughout this section, the working example is the ideal image displayed in Figure 1.1.

1.2.1

Image Contrast Enhancement

Figure 1.2 shows an image with very low contrast, which makes it difficult for human vision to clearly perceive important visual features. Assume that u is a gray-scale image in the range of 0 (black) and 1 (white). Let M denote the maximum value of u and m the minimum value. Then low contrast often means

1.2. Image Processing by Examples

7

Figure 1.1. An ideal image: noiseless, complete, and in good contrast. that M and m are not separated far enough. It may appear that a simple linear transform like can renormalize the range and suffice to enhance the image. But imagine the following scenario. Let w be a normal image with natural contrast and w(x) e [0, 1] for any pixel x. Suppose it is distorted by some nonlinear process to a low-contrast image u by

Figure 1.2. A low-contrast version of the ideal image.

Chapter 1. Introduction

8

Then for u, 0.4 < m < M < 0.5, and indeed the contrast is quite low. Furthermore, due to the nonlinearity of the distortion, the simple linear transform u —> i> proposed above cannot work faithfully to approximate the original image w. In real scenarios, there are infinitely many transforms like (1.1) that could lead to lowcontrast images. Further complexities may arise when the transforms are pixel dependent. It thus starts to become an interesting and nontrivial problem: without any specific knowledge on low-contrast transforms like (1.1), how could one design a scheme that can optimally restore the contrast, and perhaps even more importantly, what are the proper criteria for optimality? A very important solution to this problem in classical image processing is the histogram equalization technique, which can be read in many image processing books (e.g., [ 140,194]).

1.2.2

Image Denoising

Figure 1.3 shows a noisy image UQ\

where n(x) denotes additive Gaussian white noise and u the ideal "clean" image. To extract M from its noisy version MO is the denoising problem. The above additive Gaussian white noise model is popular for testing the performance of denoising models and algorithms. Another common noise model is the salt-and-pepper noise, which is nonlinear:

where u(x) e [0, 1] stands for the ideal clean image, c(x) a random field of independent 0-1 (with probabilities 50%:50%) binary random variables, and s(x) a random field of independent 0-1 binary random variables with p — Pmb(s(x) — 1) fixed for all .v e £i. c

Figure 1.3. Degradation by additive Gaussian noise (Chapter 4).

1.2. Image Processing by Examples

9

Figure 1.4. Salt-and-pepper noise with 10% spatial density (Chapter 4). and s are assumed to be independent. A sample of such noises can be simulated as follows. Take any two imaginary coins, one biased for s so that the chance of heads (s = 1) is p and the other unbiased for c. At each pixel x, first flip the s-coin. If s is a tail (s = 0), leave u(x) as it is and move to an unvisited new pixel. Otherwise (i.e., 5 = 1), further flip the c-coin, and change u(x) to the outcome of c. Repeat this game until each pixel has been visited exactly once. Then the resultant distorted image gives a sample of MO- Figure 1.4 shows an image degraded by the salt-and-pepper noise with p = 10%. Chapter 4 will further elaborate on noises and various approaches to suppress them.

1.2.3

Image Deblurring

Figure 1.5 shows a real blurry image out of focus. Blur is a common problem in satellite imaging, remote sensoring, and cosmological observations. For instance, before the Hubble telescope was launched into space in 1990, astronomic image observations had always been made by telescopes directly mounted on the surface of the Earth. As a result, images qualities were often inevitably blurred by atmospheric turbulence. Image blur occurs not only in ground telescope observations but also in other scenarios including fast motion or camera jittering (see Figure 1.6). Chapter 5 further explains in details different types of blurs, and models and algorithms to correct them, which is the problem of image deblurring.

1.2.4

Image Inpainting

Figure 1.7 shows a simulated image with 150 8-by-8 blocks randomly missing. The phenomenon often occurs in image communication where the packets are uncontrollably lost, for example, during wireless transmission. To restore the missing image information is in essence an interpolation problem, and is professionally called the error concealment problem by the communication community [166, 266].

10

Chapter 1. Introduction

Figure 1.5. Degradation by out-of-focus blur: the digital camera focuses on a fingertip about 1 inch away while the target scene is about 1 foot away (Chapter 5).

Figure 1.6. Degradation by motion blur due to horizontal hand jittering during a single exposure (Chapter 5).

In Chapter 6, readers see that error concealment is just one particular example of general image interpolation problems, which arise from numerous important applications including restoration of ancient paintings, movie production, image compression and construction, zooming and superresolution, and visual interpolation. The word "inpainting" is the artistic synonym for image interpolation and has been circulating among museum restoration artists for a long time [112, 312]. It was first introduced into digital image processing in the work of Bertalmio et al. [24] and later further popularized in many other works [23, 67, 61, 116].

1.2. Image Processing by Examples

11

Figure 1.7. 150 8-£>j-8 packets are randomly lost during transmission (Chapter 6). The goal of error concealment (or more generally, inpainting) is to develop models and algorithms that can automatically fill in the blanks [24, 67, 166].

1.2.5

Image Segmentation

Image segmentation crucially bridges low-level and high-level computer visions. Figure 1.8 shows an ideal image and its visually meaningful partition (or segmentation) of the entire image domain into different regions and boundaries. Almost trivial to human vision, image segmentation still remains one of the most challenging and most studied problems in image processing, image understanding, and artificial intelligence. A general image segmentation can be formulated as follows. Given on a 2-D domain £1 an image observation UQ, which could be degraded by noise or blur, find a visually

Figure 1.8. Such cartoonish segmentation seems trivial to human vision but still remains the most fundamental and challenging problem in image processing and low-level computer vision (Chapter 1). (The background segment QQ is not shown explicitly here.)

12

Chapter 1. Introduction

meaningful partitioning of the image domain,

for some image-dependent N, such that each component £2n (n > 1) visually corresponds to an "object" (besides the "background" patch S70)- Beneath this intuitive description of the segmentation problem lie the unavoidable complications on how to properly characterize the notions of (1) "visually meaningful," (2) image-dependent N, and even (3) "objects." These questions have been partially answered, e.g., in the celebrated works of Geman and Geman [130] and Mumford and Shah [226].

1.3

An Overview of Methodologies in Image Processing

In this section, we give an overview of some major approaches to image processing. It is unfair to say which one is necessarily superior to all the others. The efficiency and advantages of a particular methodology often depend on the concrete tasks at hand, as well as the classes and data structures of the images provided.

1.3.1

Morphological Approach

Since images capture objects, image processing naturally studies operations of objects. 2-D objects can be modelled as sets or domains of the 2-D plane R2 in a continuum setting or of the 2-D lattice Z2 in the discrete or digital setting. Equivalently, an object A can be identified with its binary characteristic function:

A morphological transform T is a map among objects,

which is often local in the sense that whether or not a pixel x belongs to B can be completely determined by the local behavior of A in the vicinity of x. Two most fundamental morphological transforms are the dilation operator D and erosion operator E, both depending on a structure element 5, or equivalent!y, a local neighborhood template. For example in the discrete setting, one could take

which is a 3-by-3 square template. The associated dilation and erosion operators are defined by [140]

Here 0 denotes the empty set and y + S = {y + s \ s e S}. It is evident that

1.3. An Overview of Methodologies in Image Processing

13

as long as (0, 0) e S, from which come the names "erosion" and "dilation." It is evident that both operators are monotonic:

Dilation and erosion provide convenient ways to define the "edge" or boundary of objects. For instance the dilative edge of an object can be defined by

while the erosive edge can be defined by

Edges are very crucial in visual inference and perception [163, 210]. Figures 1.9 and 1.10 show the effect of dilation and erosion on a binary object. In applications, two other morphological transforms generated from dilation and erosion play even more important roles. The first is the closing operator C$ — E$ o DS, i.e., with dilation followed by erosion. When applied to shapes, C$ can close small holes or gaps. And the second is the opening operator O$ = DS o E$, i.e., with dilation following after erosion. The net effect is that O$ can often erase narrow connectors and open up large holes or caves. The above binary morphological transforms can be naturally extended to general gray-level images via their level sets. Take the dilation operator DS for example. Let u be a general gray-level image defined on the lattice Z2. For each gray-level A., define the cumulative A-level set of u by

Then one defines another gray-level image v = D$(u], a dilated version of u, by

Figure 1.9. A binary set A and a structure element S (for dilation and erosion).

Chapter 1. Introduction

14

Figure 1.10. Dilation DS(A) (left) and erosion E S ( A ) (right): dilation closes up small holes or gaps, while erosion opens them up. By the monotonicity condition (1.2), FA (v) defined in this way indeed satisfies the necessary condition for cumulative level sets:

Therefore v is uniquely determined by the reconstruction formula

1.3.2

Fourier and Spectral Analysis

Fourier or spectral analysis has been one of the most powerful and favored tools in classical signal and image processing. If an image u is considered as a continuous function on the canonical rectangular domain £2 = (0, I) 2 , with periodic extension, the information of u can be completely encoded into its Fourier coefficients:

The completeness is in the sense of L 2 . Furthermore, in the digital setting when the image domain is actually a finite lattice Q — (0 : N — 1) x (0 : Af — 1) and the image u = (uj) = («;,./,) with j = (ji, 72) e Q a square matrix of data, one resorts to the discrete Fourier transform (DFT):

DFT is actually an orthonormal transform after a multiplicative scalar renormalization. Furthermore, DFT allows fast implementation known as Ihefast Fourier transform (FFT), which greatly facilitates numerous computational tasks in signal and image processing [237] (see the example in Figure 1.11). As in the spectral analysis of linear differential equations, Fourier transforms and their variations (e.g., cosine or sine transforms, and Hilbert transforms [237]) have been

1.3. An Overview of Methodologies in Image Processing

15

Figure 1.11. A digital image and its discrete Fourier transform. The example reveals a couple of salient features of Fourier image analysis: (1) most high-amplitude coefficients concentrate on the low-frequency band; (2) dominant directional information in the original image is easily recognizable in the Fourier domain; and (3) the coefficients, however, decay slowly for Heaviside-type directional edges (i.e., a jump line with distinct constant values along its two shoulders). universally applicable in many image processing tasks, including linear filtering and filter design, shift-invariant linear blurs, as well as classical image compression schemes such as the old JPEG protocol. In image analysis and processing, the Fourier approach has, however, been greatly challenged since the 1980s by another even more powerful tool—wavelet analysis.

1.3.3 Wavelet and Space-Scale Analysis Fourier transform mixes long-range spatial information, which makes it less ideal for handling localized visual features like edges. One extreme example is Dirac's delta image:

Then its Fourier coefficients are

That is, all the Fourier coefficients respond indiscriminately despite such a simple image. The harmonics are therefore inefficient in representing and coding localized image information [122]. On the other hand, cognitive and anatomic evidences in vision research have shown that human vision neurons are cleverly organized to be able to resolve localized features more efficiently [122, 152]. Wavelets in a certain sense ideally embody the idea of locality by resorting to localized bases, which are organized according to different scales or resolutions [96, 203]. Simply put, the wavelet representation of a given image u is defined by where i/ra's are wavelets, with the index set A resolving scale-adapted spatial locations.

Chapter 1. Introduction

16

Figure 1.12. An example of a (mother) wavelet by Daubechies' design [96]. Localization and oscillation are characteristic to all wavelets. The pioneering works of Nobel laureates Hubel and Wiesel [ 152] revealed the remarkable physiological fact that the simple and complex cells of human vision system behave like differentiators. Likewise, each wavelet \(fa possesses the similar differentiation property

implying that it remains dormant to featureless constant images. More generally, wavelets often satisfy the so-called vanishing-moment condition:

for some positive integer m. Coupled with spatial localization (i.e., rapid decay at infinities or being compactly supported [96, 204, 290]), this annihilation property immediately implies that the wavelet coefficients of a generic piecewise smooth image are mostly negligible except for those along important visual cues such as jumps or edges. Thus wavelets are efficient tools for adaptive image representation and data compression. Figure 1.12 displays a typical example of Daubechies' wavelets.

1.3.4

Stochastic Modeling

For images with notable stochastic nature, statistical methods become more suitable than deterministic approaches. There are two different sources leading to the statistical nature of a given image observation UQ. (a) UQ is the composition of an ideal (and often deterministic) image u with some random effect X: where the function F can be either deterministic or stochastic. For example, F(u, X) = u + X, with X = n denoting the Gaussian white noise, represents a noisy image. Or, when X denotes a spatial Poisson process and F is defined by

MO becomes a stochastic image with salt-and-pepper noise.

1.3. An Overview of Methodologies in Image Processing

17

(b) Treat individual images as the typical samples of some random fields, as in Geman and Geman's celebrated work [130] on Gibbs and Markov random fields. Most natural images such as the pictures of trees, cloudy skies, sandy beaches, or other natural landscapes are often more properly handled by the stochastic framework. For images with statistical nature, stochastic methods are the most ideal and relevant tools for modeling image processing. Especially important are statistical pattern theory [143, 224], learning theory [252], Bayesian inference theory for signal and parameter estimations [176], and stochastic algorithms such as Monte-Carlo simulation, simulated annealing, and the EM algorithm [33]. Combined with filtering techniques, a stochastic approach can become even more powerful in transformed feature spaces than directly in pixel domains. The Bayesian framework is worth extra highlighting among all stochastic methods due to its fundamental role. Let Q denote some observed data and F some hidden features or patterns embedded within that have generated Q. Bayesian inference of F is to maximize a posteriori probability (MAP):

The feature distribution probability p ( F ) is called the prior model since it specifies the a priori bias among the targeted patterns and is independent of data observation. The conditional probability p(Q \ F) is called the (generative) data model since it describes how the observation Q is distributed once F is specified. As far as MAP estimation is concerned, the denominator p(Q) is simply a probability normalization constant and plays no essential role. Without any prior knowledge,

alone would also give an optimal feature estimation, and is called the maximum likelihood (ML) estimator. In image processing, however, due to the huge degrees of freedom of images as high-dimensional data sets, prior knowledge becomes crucial for effective signal or feature estimation. The Bayesian estimator combines both the prior and data knowledge

Such a Bayesian principle frequently emerges in numerous image processing models in later chapters.

1.3.5

Variational Methods

Variational approach could be formally considered as the deterministic reflection of the Bayesian framework in the mirror of Gibbs' formula in statistical mechanics [82, 131]:

where ft = \/(kT) denotes the reciprocal of temperature T multiplied by the Boltzmann constant k, and Z = Z^ denotes the partition function for probability normalization. The

18

Chapter 1. Introduction

formula directly expresses the likelihood p(F) of a feature configuration F in terms of its "energy" E[F]. Therefore under any given temperature, the Bayesian MAP approach in (1.5) amounts to the minimization of the posterior energy:

where an additive constant (or ground energy level) independent of F has been dropped. When F and Q belong to certain functional spaces, such as Sobolev or bounded variation (BV) spaces, the minimization of this posterior energy naturally leads to a variational model (see, e.g., the monograph by Aubert and Kornprobst [15]). We must mention that here F denotes a general feature or pattern variable, which could further contain multiple components, say F = ( F \ , F~L, ..., FN). According to the telescoping formula of conditional probabilities, one has

In the important case when the feature components have a natural Markov-chain structure, i.e., one has

In terms of energy formulation, it amounts to

In later chapters, readers will witness this universal structure in several celebrated models including Geman and Geman's mixture image model [130], Mumford and Shah's free boundary segmentation model [226], and Rudin, Osher, and Fatemi's total variation-based image restoration model [258]. As an example, consider the following additive noise model:

Assume that (a) n is a homogeneous field of Gaussian white noise with zero mean; and (b) V« is a homogeneous random field of isotropic Gaussian white vectors with zero means. By the variational formulation, estimating F = u from Q — UQ can be achieved by

where the two weights are inversely proportional to the variances. The preceding stochastic language is a helpful metaphor or rationale for such variational formulations [73, 223, 276]. Figure 1.13 shows an example of applying the above variational model to denoise a smooth signal polluted with noise.

1.3. An Overview of Methodologies in Image Processing

19

Figure 1.13. A noisy \-D signal and its optimal restoration according to (1.7). The second theoretical foundation for the variational approach is the Tikhonov regularization technique for solving inverse problems [298]. Many image processing tasks carry the nature of ill-posed inverse problems, and the Tikhonov regularization technique is a powerful tool to make processing tasks well-posed. In essence, Tikhonov regularization is equivalent to the above Bayesian variational formulation:

where now the prior model E[F] is understood as the regularization term and the data model E[Q | F] as the fitting or fidelity term.

1.3.6

Partial Differential Equations (PDEs)

The successful applications of partial differential equations (PDEs) in image processing can be credited to two main factors. First, many variational problems or their regularized approximations can often be effectively computed from their Euler-Lagrange equations. Second, PDEs as in classical mathematical physics are powerful tools to describe, model, and simulate many dynamic as well as equilibrium phenomena, including diffusion, advection or transport, reaction, and so on [15, 117, 146, 154, 172, 254, 262, 305, 317]. By calculus of variation, for example, the denoising model (1.7) amounts to solving the following elliptic boundary value problem:

20

Chapter 1. Introduction

or dynamically via gradient descent marching, the diffusion-reaction equation

with some suitable initial condition u(x, t = 0). On the other hand, PDE modeling in image processing does not have to always follow a variational model. This is well known in physics, with famous examples of Navier-Stokes equations in fluid dynamics, Schrodinger equations in quantum mechanics, and Maxwell equations in electromagnetics. One most remarkable example of PDEs in image processing is the mean curvature motion (MCM) equation [40, 119, 241]

under suitable boundary as well as initial conditions. The MCM equation shares the double features of diffusion and advection. First, note that by naively crossing out the two | V0|'s, one arrives at the more familiar diffusion equation 0, = A0. On the other hand, define N = V0/| V0| to be the normal direction of level curves. Then the MCM equation (1.10) can be rewritten as

which is in the form of an advection equation with characteristics (i.e., particle motion) given by the curvature motion [117]

where at each pixel a, K(a) is precisely the (signed) scalar curvature of the level curve {* | (*) = 0( fl )}. More interestingly, let s denote the (Euclidean) arc length parameter of a level set and T = xs its unit tangent vector. Then K N = Ts = xss, and the above particle motion driven by curvature simply becomes

which by itself is a system of two heat equations coupled by the arc length s. Therefore, it is conceivable that under the MCM equation (1.10), all the level curves tend to become more and more regular. As a result, the MCM equation can be employed for tasks such as image denoising. Figure 1.14 displays the evolution of a trefoil shape under the mean curvature motion. Another important example of PDE processing model is Perona and Malik's anisotropic diffusion for image denoising and enhancement [251, 317]:

1.3. An Overview of Methodologies in Image Processing

21

Figure 1.14. A trefoil evolves under mean-curvature motion (1.11) or (1.12).

where the diffusivity coefficient D depends on the gradient information of u. D often inversely responds to p = V«|, for example, D(p) = 1/(1 + p2). The Perona-Malik model improves the ordinary constant-coefficient diffusion equation in its adaptivity to edge information: faster diffusion in interior regions where image changes mildly and slower diffusion in the vicinities of object boundaries. Thus the model can maintain sharpness of edges while diminishing the effect of noises, which is desirable for image analysis and processing. Well-posedness of the Perona-Malik model is, however, not so trivial as its intuitive behavior might have suggested (e.g., see the work by Catte et al. [51]). PDEs involving important geometric features are often called geometric PDEs. As one shall see in later chapters, geometric PDEs can usually be constructed by either optimizing some global geometric quantities in the variational setting (e.g., lengths, areas, and total squared curvatures) or by geometric invariance under certain transform groups. High order PDE models also frequently emerge in modern image processing [24, 26, 68, 61, 184, 262].

1.3.7

Different Approaches Are Intrinsically Interconnected

The current book attempts to present most modern image processing approaches and reveal their qualitative or quantitative connections. Some of such efforts have already been made in the literature (see, e.g., [270, 274, 283, 330]). To illustrate how different approaches are intrinsically connected, Figure 1.15 gives a very specific example on a classical image denoising model, which has just been mentioned in the preceding sections. In the stochastic setting (lower left panel in the figure), the noise model is given by uQ = u + n,

with n denoting Gaussian white noise of zero mean.

22

Chapter 1. Introduction

Figure 1.15. Different approaches are intrinsically connected: an example.

Given a typical (in the information-theoretic sense [93]) single observation w 0 , one attempts to remove noise n and extract the ideal image u. By the Bayesian MAP estimation theory, one attempts to maximize the posterior likelihood

Since the generative data model is straightforward in the Gaussian white noise case, the performance of such a MAP estimator is therefore mostly determined by the image prior model p ( u ) . Using Gibbs' ensemble models for images as random fields [130], the prior p ( u ) is then completely determined by its "energy" E[u]. In the case when (i) a standard Cartesian graph (or topological) structure is endowed to the underlying pixel lattice Q and (ii) the energy E[u] is built upon dipolar (or doubleton cliques; see Chapter 3) quadratic potentials, the MAP estimator becomes the variational model in the upper left panel of Figure 1.15 (after being normalized by some multiplicative constant):

1.3. An Overview of Methodologies in Image Processing

23

where the integrals are understood as discrete summations if the pixel domain £2 is a Cartesian lattice. In the continuum limit, it is naturally required that the compatible images u belong to the Sobolev class //' (Q,). Applying the first variation u -> u + 8u to the energy functional E[u \ UQ], one then obtains the generalized differential (or Euler-Lagrange equation)

in the distributional sense. Here the boundary term is understood as an element in the Hilbert space of L 2 (3£2, %'), i.e., all boundary functions that are square integrable with respect to the one-dimensional (1 -D) Hausdorff measure of 9^2. Therefore, by either gradient descent time evolution or directly solving the equilibrium equation, one arrives at the two equations in the upper right panel of Figure 1.15. Notice that the boundary term gives rise to the Neumann boundary conditions for both equations. Finally, from the equilibrium equation

one has for the optimally denoised estimator u,

This solution can be understood as a new way to decompose the given image UQ into two components: u and if. u belongs to the Sobolev space //' and is therefore a smooth or regular component, w, on the other hand, is oscillatory since it is the distributional Laplacian of u: w — — Aw j\ and generally belongs only to L 2 (£2). w is large along edges where u experiences noticeable jumps. In combination, locality (since differentiation is local), oscillations, and strong responses to edges qualify the w-component as a generalized wavelet projection of UQ that encodes detailed features. By analogy, the smoother w-component behaves more like a coarse-scale projection in the multiresolution setting of wavelet theory [96, 204, 290]. (The two are, however, not exactly orthogonal.) This brings us to the lower right panel of Figure 1.15. Naive as it may have sounded, this approach indeed led Burt and Adelson [35] to introduce the famous Gauss/Laplace pyramid algorithm for image coding and compression, which has been acknowledged as one of the few earliest wavelets and multiresolution ideas [96, 204, 290]. In Burt and Adelson's pioneering work [35], one should make the following adjustment in (1.13):

where Ga — G\(x/a}/o2 with G\(x} — G\(x\, jc2) denoting the 2-D canonical Gaussian. A similar observation was also made earlier by Gabor in 1960 in the context of image deblurring, as pointed out by Guichard, Moisan, and Morel [146]. In Chapter 3, readers shall also be able to see that for Besov images, there indeed exist precise variational formulations that can lead to Donoho and Johnstone's thresholding-based

24

Chapter 1. Introduction

denoising method in wavelet domains [106, 109]. Such formulations are due to the works of DeVore, Lucier, and their collaborators [55, 101]. Notable efforts have also been made by Steidl et al. [286] to understand the intrinsic connections among wavelets, PDEs, and variational models.

1.4

Organization of the Book

The current book can be divided into three major parts (see Figure 1.16): (A) Chapter 2: general mathematical, physical, and statistical background for modern image analysis and processing; (B) Chapter 3: several generic ways to model and represent images; and (C) Chapters 4,5,6, and 7: models and computation of four most common image processing tasks: denoising, deblurring, inpainting or interpolation, and segmentation. Chapter 2 introduces several general topics significant for modern image analysis and processing, which can be very helpful for understanding the materials of later chapters. It also provides independent reference sources for most contemporary works in the field. These topics are: (a) differential geometry of curves and surfaces in two and three dimensions; (b) the space of functions with bounded variations (BV); (c) elements of statistical mechanics and their implications for image analysis; (d) an introduction to the general framework of Bayesian estimation theory; (e) a compact theory of filtering and diffusion; and (f) elements of wavelet theory.

Figure 1.16. Organization of the book.

1.4. Organization of the Book

25

On one hand, these topics reveal the broad scopes of modern image analysis and processing, and on the other, they make the current book more self-contained. Readers can find numerous other monographs for these six individual topics, but few have been exclusively and systematically devoted to image analysis and processing. Chapter 3 introduces several ways to model and represent images, which in [273] the second author has characterized as the fundamental problem of image processing. Proper mathematical models for images are crucial for designing efficient models and algorithms to process them. Some image models are more general and universal than others. It is, however, often difficult to say which one offers the best model or representation, since optimality is inseparable from the concrete tasks at hand. This is quite similar to the particle-wave duality in quantum mechanics. Particle theory and wave theory are just two different approaches to modeling the same phenomenon of microscopic particles. Particle theory is more convenient in explaining behavior like the mass of photons according to Einstein's equation E = me2, while wave theory becomes superior in other occasions such as to explain the diffraction phenomenon of electrons through slim slits. To conclude, each model or representation has its own strengths and limitations. In Chapter 3, we first discuss several deterministic image models in the context of real or functional analysis, including distributions or generalized functions, L^-functions, Sobolev functions, functions with BV, as well as Besov functions for which the multiscale nature of images can be conveniently revealed using wavelet bases. We then turn to stochastic models of images, as originated from Geman and Geman's celebrated work on Gibbs' image models [130]. It is also discussed how to learn image field characteristics by combining filtering techniques with the maximum entropy principle, as first proposed in the remarkable work of Zhu and Mumford [328] and Zhu, Wu, and Mumford [329]. This chapter concludes with two other geometry-based image models: the level-set representation and Mumford and Shah's free boundary image model. These two models play fundamental roles in designing and understanding numerous variational or PDE models in modern image processing. Based on these preparations, the remaining four chapters study individually four significant image processing tasks and their models: noise and denoising (Chapter 4), blur and deblurring (Chapter 5), inpainting or image interpolation (Chapter 6), and edge extraction and image segmentation (Chapter 7). All these four tasks are in essence inverse problems. We thus usually start with explaining the corresponding "forward" problems, such as the precise meaning or models of noises, the physical mechanisms leading to various blurs and their mathematical models, as well as how to design generative or synthetic image models leading to multiple segments in general images. We then explore various approaches to their inverse problems, i.e., suppressing embedded noises (denoising), reducing blurry effects (deblurring), completing lost image information (inpainting), and partitioning into different objects (segmentation). To sum up, the current book has been organized in the logical order of

and the central line of philosophy or conclusion is

26

Chapter 1. Introduction

These characteristics well echo the Bayesian inference theory introduced in the preceding section (see (1.5)):

That is, proper prior knowledge on the target features or patterns F can introduce necessary bias among a large pool of candidates and lead to efficient schemes for estimation and decision making.

1.5

How to Read the Book

A monograph like the present book has multiple purposes: (a) to demonstrate the intrinsic logical structure behind this nontraditional novel field of image analysis and processing in contemporary applied mathematics; (b) to offer a well-organized presentation of modern image processing for interested students or scientific researchers in many related disciplines; (c) to summarize and highlight the most significant contributions by mathematicians during the past few decades, especially for our colleagues in engineering, computer sciences, astronomy, medical imaging and diagnosis, and brain and cognitive sciences; (d) to explain the vast range of mathematical as well as computational techniques and challenges in contemporary image processing, especially for general applied mathematicians working on computation, modeling, and applied analysis; (e) to reveal to the mathematicians of other areas many fascinating applications of the knowledge of topology, geometry, real, functional, and harmonic analysis, invariant theory, probability theory, partial differential equations, calculus of variations, and so on, as well as the genuine impact of mathematical power on modern societies and humanity; (f) to highlight and emphasize the interconnections among different mathematical tools and methodologies in contemporary image processing, with more weight on the variational PDE approach, an area in which the authors have made many contributions; and (g) as a whole, to serve as an effective communication channel bridging different communities related to image processing, including computer and human vision, computer graphics, signal and data analysis, statistical inference and learning, and numerous aforementioned application fields in medicine, astronomy, and scientific probing and sensoring. With such a broad range of readers in mind, the authors thus by no means expect a uniformly optimal approach to reading the present book. Each individual reader should come up with his or her own comfortable usage of the book, subject to important factors such as one's background training in imaging sciences or mathematics, one's primary focuses and motivations, as well as potential time pressure, and so on. Below we try to envision several classes of readers and offer our unnecessarily optimal personal suggestions. These categories are by no means mutually exclusive.

1.5. How to Read the Book

27

Scientific Workers with Specific Tasks in Image Processing Some readers are being driven by certain specific image processing tasks at hand. Examples include an astronomy data analyst who is bothered by blurs in telescope images, a medical radiologist who attempts to remove noises in some CT images, or some security developer who engages in designing automatic tracking software for surveillance cameras. We then recommend such a reader to jump right into the image processing chapters (i.e., Chapters 4 to 7; see Figure 1.16) that are most pertinent to the problem at hand and skip over Chapters 2 and 3 in the beginning. That is, for example, the previously imagined security developer could start directly with Chapter 7 on image segmentation (since segmentation is the most essential ingredient for object tracking). After getting the general ideas and an integrated picture of those relevant chapters on image processing, the reader can then selectively consult related topics in Chapter 2 and some necessary materials in Chapter 3. Depending on his or her mathematical background, the reader could even consult some external monographs on specific topics there, for example, differential geometry [103], elements of statistical mechanics [82, 131], or partial differential equations [117, 291]. (We have made great efforts in making the current book more self-contained, but real situations often vary among different individuals.) With these investments, the reader can then frequent those image processing chapters that are directly related to the tasks at hand, gain more detailed understanding of the modeling ideas and techniques, and eventually start to develop his or her own novel versions of models or algorithms. Graduate Student Readers As far as the current book is concerned, there are mainly two "dual" classes of graduate student readers: (A) one already with some direct experience or knowledge on image processing but less mathematical training and (B) the other with more mathematical training but lacking previous exposure to imaging sciences. For example, an average graduate student in electrical or medical engineering or computer sciences falls into the first category, while a typical mathematics graduate student belongs to the second. For graduate students of class (A) who are working on some image processing project towards theses, we recommend to first read Chapter 2 on several important mathematical elements or techniques for modern image processing. Afterwards the reader can proceed directly to the image processing chapters (i.e., Chapters 4 to 7) that are intimately connected to his or her current research project. With adequate understanding of the main lines of those chapters, he or she can then further invest time and efforts in Chapter 3 in order to gain more fundamental and systematic understanding of image analysis and modeling, as well as their application in developing good image processing models. For mathematical graduate students of class (B), we recommend to follow the existent logical structure of the book. The reader can start with Chapter 2 and skip whatever he or she is familiar with and then proceed to Chapter 3 to understand how knowledge in real, functional, harmonic, and stochastic analysis can become very useful in describing various classes of images. Afterwards comes the freedom of order in reading the last four chapters on image processing since Chapters from 4 to 7 are relatively independent. The reader is then suggested to start with whatever topic he or she is more comfortable with. For example, if

28

Chapter 1. Introduction

one is familiar with function interpolation in numerical or approximation theory, the reader can first read Chapter 6 on image inpainting. However, please keep in mind that Chapter 4 on denoising contains many important techniques described earlier and is always a good starting spot if no preference of a specific reading order arises. General Mathematicians Interested in Image Processing There are a growing number of mathematicians in pure, applied, or computational mathematics who have no previous direct involvement in image processing research but wish to broaden their own horizon of view by understanding the main ideas of image processing. We then recommend to follow the logical structure of the book (Figure 1.16). They could first selectively read topics that are most recognizably connected to their own field of interest, based upon which gradually expand interest into other topics. A geometer, for example, can first read the geometric topics (Sections 1 and 2) in Chapter 2 and perhaps first skip the details of other sections in order to directly proceed to the similar topics in Chapter 3 with geometric flavors. He or she can then peruse some selected topics in Chapters from 4 to 7 to see how geometric ingredients are built into various image processing models and what the main resultant challenges are in modeling, analysis, and computation. A computational mathematician, for another instance, may be mainly thirsty for finding challenging models to compute or novel algorithms to analyze. He or she can then quickly go over both Chapters 2 and 3 and focus more on Chapters from 4 to 7 on image processing models and algorithms. To him or her, Chapter 2 or 3 can serve as either intellectual enlightenment on new areas or good reference materials on image processing. There are indeed numerous opportunities in image processing for further computational development. Other General Scientists Interested in Image Processing The group could include, for example, psychologists working on human visual perception, computer graphics designers, or scientists who have been mainly working on acoustic signal processing and just recently become curious as to how visual signals (or images) are being processed. These scientists have already been exposed to signal or image processing concepts previously but never yet had a good chance to adequately acquaint themselves with major image processing ideas and methods. They have no specific image processing tasks at hand and are willing to invest some of their time and efforts because of the kinship between their own fields and image processing. Some mathematical tools and topics presented in Chapter 2 or 3 may in the beginning sound intimidating, but these readers should always keep in mind that tools always follow the footprints of tasks. Without good understanding on the real tasks or missions, difficulty could multiply in comprehending the meaning and languages of their tools. Therefore, such readers can start from Chapters 2 and 3 to get basic tones of the book but should never be discouraged by some of the mathematical details, notations, or formulations that are temporarily beyond their comprehension. They can thus in the beginning spend more time focusing on Chapters 4 to 7 on four very specific but important image processing tasks. The tasks themselves are often relatively much easier to grasp.

1.5. How to Read the Book

29

When coming to the approaches to carrying out those tasks, the readers then become highly motivated in terms of the main challenges and key ideas. It is from this point on that our readers can start to use Chapters 2 and 3 and gradually understand in depth the real motivations, advantages or disadvantages, and the key ingredients of the mathematical entities. Depending on their mathematical background, these readers are also encouraged to read some external introductory materials on specific mathematical topics in Chapters 2 and 3. Finally, the authors wish that our readers can genuinely enjoy the numerous unique views and state-of-the-art methods presented in the book. Together we can further advance this exciting novel field which is benefiting profoundly not only contemporary sciences and technologies but also the society, civilization, and the very nature of humanity. Seeing is not passive watching or reception but active believing, inferencing, and decision making based upon conscious or subconscious computations in the human brain.

This page intentionally left blank

Chapter 2

Some Modern Image Analysis Tools

In this chapter, we help readers review several major tools in modern mathematical image and vision analysis, whose fundamental roles will manifest in later chapters. These topics reveal the geometric, stochastic, real analysis, and harmonic analysis aspects of this emerging field. In addition to providing useful background mathematical knowledge for later chapters, the present chapter can also serve as a valuable reference source for researchers in the field.

2.1

Geometry of Curves and Surfaces

Curves and surfaces are basic geometric elements in image and vision analysis, and computer graphics. For example, they define and reveal the information of an automobile (in automatic traffic control), a star or planet (in astronomic imaging), a human body (in video surveillance), and an internal organ (in medical imaging and 3-D reconstruction). In this section, we review some basic theories on curves and surfaces in two or three dimensions. Readers are also referred to the classical geometry book by do Carmo [103] or some other monographs on geometric image analysis and computation, e.g., Romeny [256] and Kimmel [173]. 2.1.1

Geometry of Curves

Local Geometry A parametric curve is often denoted by x(t), with t varying on an interval andx = (x\,X2) or x — (XI,XI,XT,) depending on the dimension. Differentiation with respect to t will be denoted by an overhead dot. We shall assume that x(t) is smooth enough, and that the parametrization is regular, meaning that x ( t ) is nonzero for all t. The local infinitesimal displacement is given by dx = xdt, leading to the Euclidean parameterization ds = \dx\ = \x\dt. Differentiation with respect to the arc length df/ds will be denoted by /' in this section for any function / defined on the curve. Then t = t(s) = x'(s) is the unit tangent vector. In two dimensions, rotation of the tangent by 90 degrees (see Figure 2.1) gives the normal vector n. More generally in 31

32

Chapter2. Some Modern Image Analysis Tools

Figure 2.1. A planar curve: tangent t, normal n, and curvature K. differential geometry, the normal is a second order feature, as manifest in the curvature formula In three dimensions, the displacement of n further introduces the third degree of local geometry, i.e., the binormal vector b and the torsion scalar T. Concisely, these combine into the well-known frame formula of Frenet [103]:

The skew-symmetry of the connection matrix results directly from the orthonormality of the moving frame (t,n,b). The torsion and binormal often appear in computer graphics and designs, while in 2-D image processing, curvature plays a more explicit role. In vision psychology, it has been observed that human subjects can easily observe and process convex or concave features of shapes and boundaries [163, 170]. That is to say, our biological vision networks are indeed able to evaluate and process the curvature information. The importance of curvature in image and vision analysis is perhaps rooted in the most fundamental laws of physics. From Newton's law in classical mechanics, Hamiltonian dynamics, and the Second Law in thermodynamics to Schrodinger equations in quantum mechanics, physics seems to be perfectly complete working with up to second order information. Take the Second Law for example; first order derivatives (of free energies) often define some key quantities such as temperature, pressure, and chemical potential, while the second order ones well explain their statistical fluctuations and introduce all the rest of the key quantities such as heat capacity and thermal compressibility. Though direct connection has not yet been established, we do believe that curvature-based vision cognition must result from the physics or biophysics of the vision neural networks (see, e.g., the work by Hubel andWiesel[152]). Under the arc length parametrization, assuming K(S) > 0 and applying the exterior product, we have

2.1. Geometry of Curves and Surfaces

33

Back to the general parametrization x(t), it gives

As an application, consider the graph curve of a function y = f ( x ) on the x-y plane. Take t = x as the natural parametrization. Then

Therefore, the curvature is given by

with positivity corresponding to convexity. For some applications in computer graphics or medical imaging (e.g., the action potential and its shock front on a normal heart), one is interested in the behavior of curves on curved 2-D surfaces instead of planes. One could generally study parametric curves x (t) on a smooth Riemannian manifold M [85, 103]. Then the tangent x belongs to the tangent space TXM, which is endowed with some local Riemannian inner product {• , •}. The arc length parameter is derived from

On a general Riemannaian manifold, high order derivatives are taken based on its LeviCivita connection D. For a local vector field X defined on a neighborhood of x e M, D X is a type (1,1) tensor in TXM, meaning that for each v e TXM, DVX e TXM. Thus the LeviCivita connection generalizes the ordinary (Euclidean) notion of directional derivative: DVX is the rate of change of X along the direction v. In particular, one could define the Levi-Civita derivative x"LC(s) = Dx>(S)x'(s). Like straight lines in Euclidean spaces, a curve with x"LC(s) ~ 0 for all s is interpreted as always heading straight forward on M, or simply, a geodesic. In image and vision analysis, most 2-D surfaces are naturally embedded in R3, and their Levi-Civita connections become more explicit. Let M denote a 2-D surface embedded in R3, and let x (t) be a parametric curve on M. Then x (t) has a natural representation under the canonical coordinates of R3: x = (x, y, z). Taking the ordinary first derivative in R3 gives the tangent x(t) & T X M . Under the arc length parameter s, the second order derivative in R3: x"(s) = (x"(s), y"(s), z " ( s ) ) often "sticks out" from the tangent plane TXM. To obtain the intrinsic or Levi-Civita second order derivative x"LC, one applies an extra step of orthogonal projection: where N is the normal of TXM in R3. In particular, a "straight line" with x"LC = 0 on M can still bear a hidden degree of bending, namely, along the direction of N. The big circles on a sphere are household examples. Furthermore, suppose that the surface M is defined as the zero level set of some function 0:

34

Chapter 2. Some Modern Image Analysis Tools

Let x(s) be a parametric curve on M with arc length s. Then (x' , V0) = 0, where V denotes the gradient. Taking differentiation again leads to

where D2 denotes the Hessian operator and the second term for the Hessian bilinear form. Notice that the surface normal is given by N = V0/| V0|. Therefore,

where the last step follows from (2.2). One clearly sees how the geometry of the surface and that of the curve interact. Variational Geometry In many situations, a collection of curves have to be dealt with simultaneously. In image segmentation for example, one attempts to select the "best" candidate among all possible boundary curves [75, 226, 299]; also for image interpolation or disocclusion in computer vision [67, 115, 116, 234], missing or incomplete boundary curves are usually fitted by some optimal candidates. In such situations, it is more ideal to treat a curve as a "point," out of an ensemble of candidate curves. In order to weigh the quality of each individual curve y, the first step is to establish some suitable measure E [ y ] . For curves generated by random walks, E can be specified by probability measures. For Brownian paths for example, the Wiener measure becomes natural and even unique under suitable assumptions. In this section, we shall deal with smooth curves, and regularity functionals are then the ideal candidates for E [ y ] . Let x(t) be a parametric curve with t e [a, b\. Then the first natural geometric measure is the length L:

At first sight, it seems necessary that x be C', or at least with L' integrable x. But the length L can be interpreted as the total variation of x(t), and it therefore suffices that x(t) belongs to the space of bounded variation (B V). Under variational processes, global geometry always points to local geometry. For length, consider a small perturbation of a given curve: x(t) —>• x(t) + Sx(t), a < t < b. By Taylor expansion, up to the linear order of Sx, the change of length 8L is given by

2.1. Geometry of Curves and Surfaces

35

Integration by parts further leads to

Thus away from the two ends, the rate of change of the tangent, —i, characterizes the sensitivity of length growth. Under the arc length parameter s,—t' = —Kn exactly defines the curvature information. A general second order (Euclidean) measure for a smooth curve y is given by

where /(/c) is some suitable function of curvature. In two dimensions for instance, /(/c) = K leads to

where 0 is the signed angle of t with respect to a fixed reference direction (such as the .x-axis) and changes continuously. For a closed smooth curve, E\ has to be an integer multiple of 2n, and is thus called the total circulation. In many situations, E\ is not a good measure since it does not penalize local oscillations as long as they cancel each other. If /(/c) = |K |, then the total absolute curvature

is exactly the total variation of the heading direction 9, It is not as memoryless as the circulation measure E\, in that local oscillations do not cancel each other. The total curvature measure £2, though capable of memorizing local ripples, still does not effectively discriminate local behavior and global trends. For instance, E2\y] equals 2n for any simple and closed convex curve y, including one that contains kinks or corners. In order to more effectively penalize local singularities, one weighs the arc length element ds by some w(|/c|), with w being nonnegative and increasing with respect to |/c|:

In the case when W(\K\) = \K\P~] for some p > 1, one obtains

A special but very useful case is Euler's elastica measure Ee\:

36

Chapter 2. Some Modern Image Analysis Tools

for some positive constant weights a and b. It first appeared in Euler's work in 1744 on modeling torsion-free elastic rods and has been reintroduced into computer vision by Mumford [222] for measuring the quality of interpolating curves in disocclusion. The ratio a/b indicates the relative importance of the total length versus total squared curvature. Variational calculus leads to the equilibrium equation of the energy Ee\:

which is nonlinear and could be solved explicitly using elliptic functions, as Mumford developed in [222].

2.1.2

Geometry of Surfaces in Three Dimensions

Local Geometry In abstract differential geometry, surfaces are 2-D Riemannian manifolds on which inner products are defined pointwise and smoothly. Their local differential structures are identical to R2. Klein bottle is a well-known example which, however, cannot be faithfully visualized infl 3 [103]. In computer graphics and 3-D image reconstructions, surfaces are always naturally embedded in their 3-D environment /?3, which makes their theory and computation relatively simpler. There are two alternative ways to define a surface embedded in /?3, parametric or implicit. A parametric surface is defined as a smooth map x(u, v) from a parameter domain ft c R2 to /?3: with the nondegeneracy condition xu x xv ^ 0. As a familiar example in image analysis, a color image could be considered as a parametric surface in the RGB color space, with the two Cartesian coordinates in the image plane as natural parameters. An implicit surface, on the other hand, is defined as a zero level set: (#) = 0 for some suitable function 0, often with the nondegeneracy condition V0 ^ 0 on the surface. By the implicit function theorem, an implicit surface could be locally parameterized. On the other hand, a parametric surface could also be made implicit locally using the distance function. Let £ c R3 denote the image of a parametric surface x(u, v). Define the distance function d^ by Then £ is exactly the zero level set of d-z. The distance function has many applications in image and vision analysis, especially in the level-set method of Osher and Sethian [241]. Consider a parametric surface x(u, v) : £1 —» R 3 . The tangent plane Tx is spanned by xu andxv. Define

Then the inner product on Tx is completely determined by the 2-by-2 structure matrix M — (E, F\ F, G), which contains all the first order information of the surface. A parametric

2.1. Geometry of Curves and Surfaces

37

curve on the surface is in the form of x(u(t), v ( t ) ) . Computing its arc length leads to the so-called first fundamental form /(-, •) of the surface:

Generally, I(u, v) is precisely the inner product for all tangent vectors u, v e TX*L. To study the second order geometric information of a surface E, one has to understand how to take derivatives for tangent vector fields. In general Riemannian geometry as already mentioned, this is realized by the Levi-Civita connection, or the covariant differential operator D. For a given tangent vector field X on Z and a direction v e TxE, DVX e 7^ E generalizes the ordinary notion of directional derivative of X along v. Due to the product rule of the covariant derivative,

and the linearity of Dv on v, it suffices to understand the following four basic actions:

If the surface is flat, these are simply the ordinary second order derivatives xuu, xuv, and X vv

For a general embedded surface, xuu may stick out from it. But creatures living on the surface cannot sense what is away. Therefore, xuu needs to be projected back onto the surface: where N denotes the unit normal in the direction of xu x xv. (If the codimension of the embedding is not 1, N has to be replaced by the entire normal space.) This discussion certainly applies to the other three as well. The quantities (x.. , N) completely measure the deviation of the covariant derivatives from the ordinary ones. Therefore they encode the local information of how far the surface deviates from being a plane, i.e., the bending and curving of the surface. Conventionally, they are denoted by

Since (x. , N) = 0, one has

In combination, they define the second order geometric information and lead to the second fundamental form //, which is defined as follows. By abusing the covariant derivative symbol D, for any v e Tx E, let us denote by Dv N the ordinary directional derivative of N along v in R3. Since along v (i.e., along a curve segment on E which is tangent to v at x), (N , N) = 1, one must have (DVN , N) = 0, implying that DVN e T^E! Thus D.N could be considered as a linear transformation in TXY. Indeed, —D.N is called the Weingarten linear map [103]. The second fundamental (bilinear) form in Tx E is then defined as

38

Chapter 2. Some Modern Image Analysis Tools

Figure 2.2. Second order local geometry: surface normal N and two perpendicular principle curvature circles. Noticing that DXu N = Nu and DXiN — Nv, one has

In particular, dx = xudu + xvdv and

The second form is thus symmetric as well. In order to understand how the second fundamental form encodes second order geometric information, consider a unit vector u e Tx S. The two orthonormal vectors u and N span a unique plane though jc, which locally intersects Z along a planar curve y. Applying Frenet's frame formula (2.1) to y at x, one has

In the current context, one has

Therefore, II(u, u) is exactly the curvature of y, which is denoted by KU and gives the explicit geometric meaning of the second fundamental form. Furthermore, one can introduce the Rayleigh quotient for the symmetric second form as in linear algebra [138, 289]. For any u e TXZ [103], define

Generically, its maximum K\ and minimum KI are called the principal curvatures (see Figure 2.2). Gaussian curvature is defined as K = K\K2, while mean curvature H = (K\ + K2)/2. From the perception point of view, it is unclear to which degree of precision the human vision system can detect these two curvatures. Psychological evidences show that human vision is at least good at detecting the sign of K, which corresponds to the visual differentiability between ellipticity (i.e., convexity) and hyperbolicity (i.e., saddle structures).

2.1. Geometry of Curves and Surfaces

39

Let M and L denote the two structure matrices of the first and second forms under the natural parametric basis (xu, xv):

Then from the computational point of view, the two principal curvatures are exactly the two eigenvalues of the generalized eigenvalue problem (L, M): Lu = KMu. Therefore, by linear algebra [139, 289],

or more explicitly,

The ubiquitous denominator EG — F2 = det M is the squared area of the parallelogram spanned by xu and xv. In particular, if xu and xv just happen to be orthonormal at x, then

The graph of a gray-scale image h(u, v) : Q -> R could be considered as a parametric surface in (u, v, h) e R3. Then it is easy to establish that

which is the determinant of the Hessian of h, normalized by the graph area element yi + hi + h\. Similarly for the mean curvature,

Figure 2.3 shows a surface patch and its associated Gauss and mean curvatures. Now consider an implicit surface £ defined by a potential function 0:

Working near a fixed target point :c, one may take any parametrization x (u, v) locally near x on E. In addition, assume that at this target point (only!), (xu, xv) is an orthonormal system of Tx E. Differentiating (f)(x(u, v)) — 0 gives the first order information:

Taking differentiation again to the first equation with respect to u gives

40

Chapter 2. Some Modern Image Analysis Tools

Figure 2.3. A surface patch z = h(u, v) = cos(w) cos(t') (top), its mean curvature field H (lower left), and Gauss curvature field K (lower right). where D20 denotes the Hessian matrix in /?3 and its associated bilinear form. Noticing that the surface normal /V = V(/>/| V|, we obtain the formula for e in the second form:

Similarly, one has

Since (xu, xv) is an orthonormal frame at the target point x, (2.5) leads to

These two beautiful formulae are independent of the parametrization and completely selfcontained. The mean curvature H has been frequently studied in dynamical processing of shapes and surfaces, mainly due to its variational meaning (see the coming subsection). In the literature, another more explicit formula for H is also commonly used based on the trace formula The first term is exactly the 3-D Laplacian

Therefore, —2H equals

and the mean curvature is thus the divergence of the normal field N (up to a constant)!

2.1. Geometry of Curves and Surfaces

41

Variational Geometry One common global feature of a given surface £ is its area A. Under the parametric form x(u, v) : £2 —»• H, the area element is given by da = \xu x xv\dudv. Thus

As for curve length, under a small perturbation,

the first order variation of A must be in the form of

where H is a vector field on E that takes values in R3. Its exact form can be worked out as follows. Let X = xu x xv. Then the normal is N = X/\X\. Notice that 8\X\ = (N , 8X) and

From the vector calculus of inner and exterior products, one has

Applying integration by parts and noticing that 8x(u, v) has been assumed compactly supported, one has

Compared with (2.10), it gives

Notice that all the four vectors on the right are in the tangent plane Tx E. Therefore the two exterior products are both in the normal direction N. Assume

for some scalar H. Then by (2.11),

42

Chapter 2. Some Modern Image Analysis Tools

Notice that the denominator is exactly EG — F2 from the first fundamental form. Applying the two rules in vector calculus

one establishes that

where e, f , and g are the coefficients of the second form. Referring to the mean curvature formula (2.4), we end up with Therefore, we have shown that

which is the global (or variational) meaning for the local mean curvature H. It also indicates that the most area-efficient displacements are along the normals, which well coincides with the daily life intuition. Consider, for example, the (top) surface area of a glass of still water. Smooth and slow rotation (i.e., 8x's are on the tangent planes) brings no noticeable change of surface area, but any wave ripples in the vertical direction do. A surface with mean curvature H = 0 is called a minimal surface, which has been a fascinating subject in geometry, geometric measure theory, and nonlinear PDEs ever since J. L. Lagrange [233]. Minimal surfaces are also very useful in computer graphics and geometric designs [92]. In image analysis, one often considers the graph surface Z/, of a given gray-scale image h(u,v) on an image domain £2. Then the area is given by

where the convenient notation |X| fl denotes ^/a2 + \X\2 for any vector X . In this situation, we are interested in the sensitivity of A(h) with respect to small changes of h: h —> h+8h. Noticing that 8\X\\ = ( X / \ X \ \ , 8X}, we have for the first order variation of A(h)

where the perturbation 8h has been assumed to be compactly supported. It was first derived by Lagrange in 1762 (see more detailed discussion in [233] for example). Notice that for image graphs, x = (u,v,h(u, v}}, 8x = (0, 0, 8h), and (N , 8x} = 8h/\Vh\\. Comparing (2.14) with the most general form (2.13), one concludes that the surface mean curvature is given by

2.1. Geometry of Curves and Surfaces

43

For smooth surfaces, we could consider high order global features as well. For example, the total curvature

the total absolute curvature

and the total squared curvature

More generally, let f ( K , H) = g(K\, ^2) be a proper function of the two principal curvatures. Then one could attempt to minimize the cost function

Many interesting results exist regarding these energies. For example, suppose the Gaussian curvature K does not change signs. Then from the Gauss map theorem [103]

one has

which is exactly the total parametric area of the normal map N : £2 -> S2. Suppose N is injective. Then this latter area could never be bigger than 4jr, the total area of S2, which implies that EI < 4rt in this situation. Sometimes one may also become interested in energies that are not intrinsically geometric, meaning that they do depend upon specific representations. For example, for a graph surface representation x(u, v) = (u, v, h(u, v)), one could attempt to minimize

By the total curvature formula (2.6),

which is not intrinsic due to the h-factor but is indeed a close approximation to the total squared curvature £3 when Vh is small. Finally, we would like to mention that besides curves and surfaces, high-dimensional differential or Riemannian geometry could also be very useful in image processing (see, e.g., [174]).

44

2.1.3

Chapter 2. Some Modern Image Analysis Tools

Hausdorff Measures and Dimensions

In image and vision analysis, one has to deal with curves and surfaces that are less regular, for example, fractal landscapes. Then differential geometry often demands too much, while the notions of Hausdorff measures and dimensions can become more powerful. Let E be a given set in R2 or /?3, or more generally any subset in a metric space with distance function. Suppose it is known a priori that £ is a d-dimensional object, with d = 1.618 say. Then a natural question is how to measure its J-dimensional "volume." Hausdorff measure Jid provides one good answer. A covering A of E is a (countable) collection of subsets A's whose union includes E as a subset. The scale ||^4|| of A is defined as

where diam(A) stands for the diameter of A. Furthermore, one defines

Then the J-dimensional Hausdorff measure 1-Ld(E) is defined as

where A's must be coverings of E. Notice that the e limit is well defined since the infimum as a function of € is nonincreasing. That is to say, the € limit process is in fact a supremum one. Thus the definition is a minimax statement. It is easy to see that all the ingredients in the definition are natural. (a) When d is an integer and A is an ordinary ^/-dimensional cube or ball, diam( A)d correctly measures the volume up to a multiplicative constant. (b) Taking infimum is to diminish the effect caused by inefficient overcovering or repetition. (c) Finally, by demanding 6 to vanish, a Hausdorff measure can potentially resolve details at all scales (see Figure 2.4). From the measure-theoretic point of view, T-Ld defined above is only an "outer" measure. Any outer measure, once restricted on the class of measurable sets, becomes a true measure. It is easy to see that for any d,t > 0,

from which one has, for any 6 > 0, Therefore, as long as T-Ld(E) is finite for some d, H.S(E} = 0 for all s > d. Therefore, there exists at most one d, so that l-Ld(E) is nonzero and finite, which leads to the notion of Hausdorff dimension dim//:

If dim//(E) > 0, then for any t e [0, dim//(E)), n'(E) = +00 (see Figure 2.5).

2.2. Functions with Bounded Variations

45

Figure 2.4. A set E (left) and a covering A (right). Covering elements like the huge box with dashed borders (right) fail to faithfully capture small-scale details. The definition of Hausdorjf measures forces covering scales to tend to zero.

Figure 2.5. The Hausdorff dimension dim//(£) of a set E is the critical d, above which E appears too "thin," while below it appears too "fat. "

2.2

Functions with Bounded Variations

Functions with bounded variations (B Vs) are ideal deterministic image models which allow the existence of jumps or edges, and are, however, still mathematically tractable. Theory of BV functions and its extensions [7, 9, 118, 195, 137] provide the necessary mathematical foundations for celebrated image processing models such as Rudin, Osher, and Fatemi's total variation denoising scheme [257, 258] and Mumford and Shah's segmentation model [226].

46

2.2.1

Chapter 2. Some Modern Image Analysis Tools

Total Variation as a Radon Measure

Let £2 c R2 denote a bounded open domain and u = u(x, y) belong to L 1 (£2). Motivated by both analysis and image processing, £2 is often assumed to be a Lipschitz domain. If u is smooth, then its total variation (TV) TV[w] is defined as

A function (or an image) u with TV[w] < oo is said to have B V. The notation B V(ft) denotes all functions in L 1 (ft) with BV. How smooth should a function be for the definition (2.18) to make sense? At first sight, it appears that at least u should belong to the Sobolev space WIA (ft), i.e., functions with integrable first order partial derivatives. But the power of TV and BV in image analysis exactly arises from the relaxation of such constraints. That is, the space BV(ft) is in fact much larger than W 1 - 1 (ft) to include many interesting image functions. To get acquainted with such relaxation, let us first look at a related but simpler example. For any function f(x, y) e L 1 (ft), we could define the total integral / by

The condition L 1 seems natural but not exclusively necessary for /[/] to be well defined. For example, assume that (0, 0) € ft, and denote by 8(x, y) the Dirac delta function. Then one can still define More generally, let \JL be a nonnegative measure on all Borel sets of ft with p.(K) < oo for all compact subsets K c ft. Such measures are examples of Radon measure [126, 118]. Then one could still define, as for Dirac's delta,

where /z could be generally a signed Radon measure. Signed Radon measures generalize the space of L1'oc(ft). Furthermore, in a certain proper sense (i.e., linear functionals) [126], it could be shown that signed Radon measures are indeed the most general class of mathematical entities for /[•] to be well defined. The same story applies to the TV functional. In some sense, TV[w] is the integrated version of /[/] by taking / = VM. Consider the case in Rl when /(*) = 8(x). Then u' = 8 would imply that u = H(x)—the Heaviside 0-1 function, provided that u(—00) = 0. Since H ( x ) contains a jump at the origin, it stays outside the locally Sobolev class H^,1. Inspired by the discussion on /, one could, however, still define

Generally, let ^t(jc) be any finite and signed Radon measure in /?', and let u = u(x) denote its cumulative distribution function (c.d.f. in probability theory). One could then define

2.2. Functions with Bounded Variations

47

Figure 2.6. Three \-D images with TV[/] = TV[g] = TV[/z] = 2. Here |/z| denotes the TV measure of /JL. That is, if /z = /JL+ — ^T denotes the Jordan decomposition into two mutually singular and nonnegative Radon measures [126], then \fji\ = ^L+ + IJL~. Figure 2.6 shows the TVs of three 1-D signals. With these insights in mind, we now formalize the definitions of TV[w] and BV(£2). Let U c Q c R2 denote any open set of £2. For any u e Ll (£2), define

where B2 denotes the open unit disk in /?2, and the admissible class C\(U, B2) denotes

If / | D*u | < oo for any open subset U with compact closure contained within £2, we say that locally u has BV and denote by BVioc(£2) the collection of all such functions in Ll(Q). For any u e BVioc(£2), by the standard practice in measure theory, one could construct an outer measure by defining, for any subset E c Q,

where U are open subsets. Restricted on measurable sets, it becomes a Radon measure on Q, which shall be denoted by / \Du\, and defines the TV of u (on £2) by

The notation BV(£2) denotes the collections of all functions with finite TVs in L 1 (fi). Assume that u e W M (n). Then for any g e C'(^, fi2),

48

Chapter 2. Some Modern Image Analysis Tools

by the Gauss-Green divergence theorem (or integration by parts). Since B2 is closed under radial reflection g -» — g, one has

On the other hand, by Lusin's continuity characterization, Uryson's extension lemma, and the mollification technique [126], for any specified precision e and -) = AH(x — 1/2), a product of an amplitude constant A > 0 and the shifted Heaviside function. It represents a binary black-white image of two adjacent vertical bands. Then

where D[ means 1-D TV. It is not a coincidence in this case that the TV measure is the product of the size d — c and strength A of the target edge. Example. Let £2 = /?2, and let aB2 denote the disk centered at the origin with radius a > 0. Let XaB2 = XaB2(x, y) be its indicator function. Then by definition,

where 3 (A) stands for the (topological) boundary of a set A, K 1 the 1-D Hausdorff measure, and n the outer normal. On the other hand, it is easy to construct a g e C\ (R2, B2), so that

2.2. Functions with Bounded Variations

49

Figure 2.7. An example of L 1 -lower semicontinuity. The sequence (un) of 1D images on [0, 1] converges to u = 0 in L 1 smce ||«n+i — u\\i\ < 2~n. Notice that TV(«) = OwhileTV&n) = 2, which is consistent with the property of lower semicontinuity: TV(w) < liming TV(« n ). In particular, strict inequality is indeed realizable. This result holds for any smooth domain E; namely, the TV of its indicator function XE is exactly its perimeter. In fact, this approach leads to a more general definition of perimeters for nonsmooth domains [137, 195].

2.2.2

Basic Properties of BV Functions

In this section, we discuss several major properties of BV functions that are important in image analysis. They are L 1 -lower semicontinuity, completeness (i.e., Banach), trace and divergence theorem, smoothness and regularity, and relative compactness in L'. Theorem 2.1 ( L1-Lower Semicontinuity). Suppose un (x, v) -> u(x, y) in L' (£2). Then

In particular, if(un)n is a bounded sequence in BV(£2), then u belongs to BV(fi) as well. The proof almost immediately follows from the definition of TV. For any g e Clc(Q,B2),

Taking supremum confirms Theorem 2.1 (see Figure 2.7). Among the many interesting results that Theorem 2.1 makes contribution to, perhaps the most fundamental one is that BV(£2) is a Banach space. First, it is easy to see that BV(£2) is a normed linear space under the BV norm:

which is stronger than the L 1 norm. To establish the completeness, let (un)n be a Cauchy sequence under || • ||BV- Then (un)n is also Cauchy in L 1 , and let M e L 1 be its limit. By

50

Chapter 2. Some Modern Image Analysis Tools

semicontinuity, u e BV(ft). Furthermore, applying semicontinuity again to (u,n — un}n with ra fixed, we have

which converges to 0 as m —> oo since (un)n is a Cauchy sequence. Thus un —>• u in BV(ft) under || • ||BV. The behavior of a BV function u 6 BV(ft) near the boundary 9ft is reflected in the notion of trace fu = T(u) e L ' ( 9 f t , Hl). In what follows, we shall assume that ft is a Lipschitz domain. Roughly speaking, the definability of /„ = T (u) indicates that a B V function u cannot oscillate too much near the boundary 9ft. Instead, it converges to some function (in a certain suitable sense) / along 9£2. To see clearly why this is the case, let K\ c K^ c • • • be a sequence of compact sets so that ft = [_J^° Kn. Then by the basic continuity property of a finite measure,

Therefore, roughness indeed gradually fades away on the "narrow rings" £2 \ Kn, and the function tends to converge. More explicitly, as far as a local segment of rectifiable boundary is concerned, one could "zoom in" and simply assume that ft = R x (0, oo). Take K€ = R x (e, oo), e > 0, and define f £ ( x ) = u(x, s) in the L 1 sense (by Fubini's theorem). For any 8 < e,

Therefore ( f s ) £ is a Cauchy sequence in L' (/?), and the trace fu = T(u) along the boundary is defined as its L 1 limit. The concept of trace leads to the Gauss-Green divergence theorem for BV images. Suppose that u € BV(£2), U d £2 (i.e., U is compact in ft) and is Lipschitz, and g e C,S(ft,IR 2 ). Then

where n denotes the outer normal of dU. As a Radon measure, for any compact set K c ft,

In image analysis and processing, of particular interest is the case when K is a compact piece of Lipschitz curve y. y is said to be singular for a given image u if

Notice that under the 2-D Lebesgue measure, y is always only a null set.

2.2. Functions with Bounded Variations

51

Let / + and /~ denote the corresponding traces of u along the two sides of y. Then one can show that the integration under the Hausdorff measure of the jumps along the curve! In image and vision analysis, intensity "jumps," usually referred to as "edges," have been considered as crucial visual cues. Equation (2.21) therefore partially explains why the TV measure and the B V space could have become so convenient and powerful in image analysis and processing. Another key question is how close W 1 ' 1 ^) is to BV(fi). By adaptive mollification based on the partition of unity [126], one could establish the following result. Theorem 2.2 (Mollification of BV). For any u e BV(£2), one can find a sequence of approximations (un)n such that 1. un e C°°(Q)forn = 1 , 2 , . . . ; 2. un -H» u in L' (Q) as n —> oo; 3. JR Dun | —> fft | Du \, and their traces gn = g for n — 1 , 2 , . . . . Theorem 2.2 is the best that one could wish for, since for a general BV image u it is infeasible to ask for a sequence of smooth functions (un)n so that

This could be easily explained by (2.21). Suppose u e BV(£2) and that a compact regular curve y c £2 is singular, i.e., f \Du\ = a > 0. Then for any smooth function «„,

That is, convergence under the TV seminorm is impossible for such an image. However, mollification provided by Theorem 2.2 has already been quite powerful. It allows one to conveniently transplant many results from the classical Sobolev space Wl'l(tt). Let u € BV(£2) and / = T(u) be its trace. A functional L[u, f] on BV(£2) is said to be L 1 -lower semicontinuous if, for any sequence (un}n e BV(£2) with un —> u in L 1 (£2) and /„ = / as n —> oo, one has

Corollary 2.3 (SoboIev-BV Transition). Suppose L[u, f] is an L[-lower semicontinuous functional in BV(£2) and that E = E(t) is a continuous function oft e [0, oo). If then the same inequality must hold for all u e BV(£2).

52

Chapter 2. Some Modern Image Analysis Tools

The proof follows readily from the semicontinuity property and Theorem 2.2 on mollification approximation. As an application, define in R" (with n > 1) p = ^-j-, and

Then the celebrated Sobolev inequality in W I A ( R " ) states that L < E for all compactly supported u e WIA(R") for some suitable constant c which depends only on n. By the preceding corollary, the same inequality must hold for any compactly supported u €. BV(£2) as well. The second most important application of the mollification theorem (Theorem 2.2) is the weak compactness property, which has become a fundamental tool for establishing the existence of solutions to various BV-based image processing models. Theorem 2.4 (Weak Compactness of BV). Let (un)n be a bounded sequence in BV(£2) where Q is a Lipschitz domain. There must exist a subsequence which converges in L 1 (Q). The proof is again immediate from the mollification theorem (Theorem 2.2) and the weak compactness of W L 1 (£2) in L'(£2). By Theorem 2.2, we could approximate un by wn € W{A(Q) such that

for each n. Thus (wn)n must be bounded in WL[(Q). By the weak compactness property of W L { ( £ 2 ) (or the Rellich theorem [3, 137, 193]), (u;,,),, contains a subsequence indexed by n(k), k = 1, 2 , . . . , which converges in L 1 (£2). Then the subsequence of (u,,)n with the same indices n(k) must also converge in L 1 (fi).

2.2.3

The Co-Area Formula

This beautiful formula was initially due to Fleming and Rishel [125], and De Giorgi gave the complete proof [134, 137]. It reveals the geometric essence of the TV measure, which has been fundamental for numerous applications in image analysis and processing. For smooth functions, it essentially results from a special change of variables, at least locally. Let u(x, >') denote a smooth image on a domain Q. Assume that Vw(jc 0 , >'o) / 0. Say uy(xQ, yo) 7^ 0, without loss of generality. Then by the implicit function theorem, locally near (JCQ, >'o), we could explicitly solve v:

For each given X, one could reparameterize the 1-level curve (x, y(x, 1)) using arc length 5, so that Since the reference points with 5 = 0 may also depend on A, one has

2.2. Functions with Bounded Variations

53

which is well defined at least locally. In summary, this change of variables is governed by

In the (s, A) system, taking partial derivatives with respect to s and A in the first equation gives The first identity, combined with the second one in (2.22), gives Vu = ±| Vu\(—ys,xs). Therefore, the second identity becomes

and with the Jacobian, we are able to conclude that

Assuming that this change of variables is global on £2, one has

where the second integral is over an appropriate domain in the (s, A) plane. Let y^ denote the A-level curve. Then

which is the celebrated co-area formula for regular functions. More generally, for any function

This page intentionally left blank

image

processing

and Analysis

variational PDE, wavelet, and stochastic Methods

Tony F. chan

university OF california, LOS Anggles LOS AnGeLeS, california

jianHonG [jacKie] shen

university OF Minnesota MinneapoLis, Minnesota

Society for Industrial and Applied Mathematics Philadelphia

Copyright © 2005 by the Society for Industrial and Applied Mathematics. 1098765432 1 All rights reserved. Printed in the United States of America. No part of this book may be reproduced, stored, or transmitted in any manner without the written permission of the publisher. For information, write to the Society for Industrial and Applied Mathematics, 3600 University City Science Center, Philadelphia, PA 19104-2688. MATLAB® is a trademark of The MathWorks, Inc. and is used with permission. The MathWorks does not warrant the accuracy of the text or exercises in this book. This book's use or discussion of MATLAB® software or related products does not constitute endorsement or sponsorship by The MathWorks of a particular pedagogical approach or particular use of the MATLAB® software. For MATLAB® product information, please contact: The MathWorks, Inc., 3 Apple Hill Drive, Natick, MA 01760-2098 USA, 508-647-7000, Fax: 508-647-7101, [email protected], www.mathworks.com

Library of Congress Cataloging-in-Publication Data Chan, Tony F. Image processing and analysis : variational, PDE, wavelet, and stochastic methods / Tony F. Chan, Jianhong (Jackie) Shen. p. cm. Includes bibliographical references and index. ISBN 0-89871-589-X (pbk.) 1. Image processing—Mathematical models. I. Shen, Jianhong, 1971- II. Title. TA1637.C4775 2005 621.36'7-dc22

2005048960

• 51£LJTL is a registered trademark.

T0 Monica, Micfiettej Ityan, andCtaudia

This page intentionally left blank

Contents List of Figures

Xlll

Preface

xix

1

Introduction 1.1 Dawning of the Era of Imaging Sciences 1.1.1 Image Acquisition 1.1.2 Image Processing 1.1.3 Image Interpretation and Visual Intelligence 1.2 Image Processing by Examples 1.2.1 Image Contrast Enhancement 1.2.2 Image Denoising 1.2.3 Image Deblurring 1.2.4 Image Inpainting 1.2.5 Image Segmentation 1.3 An Overview of Methodologies in Image Processing .3.1 Morphological Approach .3.2 Fourier and Spectral Analysis .3.3 Wavelet and Space-Scale Analysis .3.4 Stochastic Modeling .3.5 Variational Methods .3.6 Partial Differential Equations (PDEs) .3.7 Different Approaches Are Intrinsically Interconnected 1.4 Organization of the Book 1.5 How to Read the Book

1 1 1 5 6 6 6 8 9 9 11 12 12 14 15 16 17 19 21 24 26

2

Some Modern Image Analysis Tools 2.1 Geometry of Curves and Surfaces 2.1.1 Geometry of Curves 2.1.2 Geometry of Surfaces in Three Dimensions 2.1.3 Hausdorff Measures and Dimensions 2.2 Functions with Bounded Variations 2.2.1 Total Variation as a Radon Measure 2.2.2 Basic Properties of BV Functions

31 31 31 36 44 45 46 49

VII

viii

Contents

2.3

2.4

2.5

2.6

3

2.2.3 The Co-Area Formula Elements of Thermodynamics and Statistical Mechanics 2.3.1 Essentials of Thermodynamics 2.3.2 Entropy and Potentials 2.3.3 Statistical Mechanics of Ensembles Bayesian Statistical Inference 2.4.1 Image Processing or Visual Perception as Inference 2.4.2 Bayesian Inference: Bias Due to Prior Knowledge 2.4.3 Bayesian Method in Image Processing Linear and Nonlinear Filtering and Diffusion 2.5.1 Point Spreading and Markov Transition 2.5.2 Linear Filtering and Diffusion 2.5.3 Nonlinear Filtering and Diffusion Wavelets and Multiresolution Analysis 2.6.1 Quest for New Image Analysis Tools 2.6.2 Early Edge Theory and Marr's Wavelets 2.6.3 Windowed Frequency Analysis and Gabor Wavelets 2.6.4 Frequency-Window Coupling: Malvar-Wilson Wavelets 2.6.5 The Framework of Multiresolution Analysis (MRA) 2.6.6 Fast Image Analysis and Synthesis via Filter Banks

Image Modeling and Representation 3.1 Modeling and Representation: What, Why, and How 3.2 Deterministic Image Models 3.2.1 Images as Distributions (Generalized Functions) 3.2.2 Lp Images 3.2.3 Sobolev Images H n (W) 3.2.4 BV Images 3.3 Wavelets and Multiscale Representation 3.3.1 Construction of 2-D Wavelets 3.3.2 Wavelet Responses to Typical Image Features 3.3.3 Besov Images and Sparse Wavelet Representation 3.4 Lattice and Random Field Representation 3.4.1 Natural Images of Mother Nature 3.4.2 Images as Ensembles and Distributions 3.4.3 Images as Gibbs'Ensembles 3.4.4 Images as Markov Random Fields 3.4.5 Visual Filters and Filter Banks 3.4.6 Entropy-Based Learning of Image Patterns 3.5 Level-Set Representation 3.5.1 Classical Level Sets 3.5.2 Cumulative Level Sets 3.5.3 Level-Set Synthesis 3.5.4 An Example: Level Sets of Piecewise Constant Images 3.5.5 High Order Regularity of Level Sets 3.5.6 Statistics of Level Sets of Natural Images

52 54 55 56 58 61 61 62 64 65 65 67 70 73 73 75 76 77 80 86 91 91 93 93 96 98 98 99 99 104 107 115 115 116 117 119 122 124 126 127 127 129 129 130 131

Contents

ix

3.6

The Mumford-Shah Free Boundary Image Model 132 3.6.1 Piecewise Constant 1-D Images: Analysis and Synthesis 132 3.6.2 Piecewise Smooth 1-D Images: First Order Representation . . . 134 3.6.3 Piecewise Smooth 1-D Images: Poisson Representation 135 3.6.4 Piecewise Smooth 2-D Images 136 3.6.5 The Mumford-Shah Model 138 3.6.6 The Role of Special BV Images 140

4

Image Denoising 4.1 Noise: Origins, Physics, and Models 4.1.1 Origins and Physics of Noise 4.1.2 A Brief Overview of 1-D Stochastic Signals 4.1.3 Stochastic Models of Noises 4.1.4 Analog White Noises as Random Generalized Functions 4.1.5 Random Signals from Stochastic Differential Equations 4.1.6 2-D Stochastic Spatial Signals: Random Fields 4.2 Linear Denoising: Lowpass Filtering 4.2.1 Signal vs. Noise 4.2.2 Denoising via Linear Filters and Diffusion 4.3 Data-Driven Optimal Filtering: Wiener Filters 4.4 Wavelet Shrinkage Denoising 4.4.1 Shrinkage: Quasi-statistical Estimation of Singletons 4.4.2 Shrinkage: Variational Estimation of Singletons 4.4.3 Denoising via Shrinking Noisy Wavelet Components 4.4.4 Variational Denoising of Noisy Besov Images 4.5 Variational Denoising Based on BV Image Model 4.5.1 TV, Robust Statistics, and Median 4.5.2 The Role of TV and BV Image Model 4.5.3 Biased Iterated Median Filtering 4.5.4 Rudin, Osher, and Fatemi's TV Denoising Model 4.5.5 Computational Approaches to TV Denoising 4.5.6 Duality for the TV Denoising Model 4.5.7 Solution Structures of the TV Denoising Model 4.6 Denoising via Nonlinear Diffusion and Scale-Space Theory 4.6.1 Perona and Malik's Nonlinear Diffusion Model 4.6.2 Axiomatic Scale-Space Theory 4.7 Denoising Salt-and-Pepper Noise 4.8 Multichannel TV Denoising 4.8.1 Variational TV Denoising of Multichannel Images 4.8.2 Three Versions of TV[w]

145 145 145 147 150 151 153 155 156 156 157 159 160 160 163 165 171 174 174 175 175 177 178 183 185 191 191 194 198 203 203 204

5

Image Deblurring 5.1 Blur: Physical Origins and Mathematical Models 5.1.1 Physical Origins 5.1.2 Mathematical Models of Blurs 5.1.3 Linear vs. Nonlinear Blurs

207 207 207 208 214

x

6

Contents

5.2 III-posedness and Regularization 5.3 Deblurring with Wiener Filters 5.3.1 Intuition on Filter-Based Deblurring 5.3.2 Wiener Filtering 5.4 Deblurring of BV Images with Known PSF 5.4.1 The Variational Model 5.4.2 Existence and Uniqueness 5.4.3 Computation 5.5 Variational Blind Deblurring with Unknown PSF 5.5.1 Parametric Blind Deblurring 5.5.2 Parametric-Field-Based Blind Deblurring 5.5.3 Nonparametric Blind Deblurring

216 217 217 218 220 220 222 224 226 226 230 233

Image Inpainting 6.1 A Brief Review on Classical Interpolation Schemes 6.1.1 Polynomial Interpolation 6.1.2 Trigonometric Polynomial Interpolation 6.1.3 Spline Interpolation 6.1.4 Shannon's Sampling Theorem 6.1.5 Radial Basis Functions and Thin-Plate Splines 6.2 Challenges and Guidelines for 2-D Image Inpainting 6.2.1 Main Challenges for Image Inpainting 6.2.2 General Guidelines for Image Inpainting 6.3 Inpainting of Sobolev Images: Green's Formulae 6.4 Geometric Modeling of Curves and Images 6.4.1 Geometric Curve Models 6.4.2 2-, 3-Point Accumulative Energies, Length, and Curvature . . . . 6.4.3 Image Models via Functionalizing Curve Models 6.4.4 Image Models with Embedded Edge Models 6.5 Inpainting BV Images (via the TV Radon Measure) 6.5.1 Formulation of the TV Inpainting Model 6.5.2 Justification of TV Inpainting by Visual Perception 6.5.3 Computation of TV Inpainting 6.5.4 Digital Zooming Based on TV Inpainting 6.5.5 Edge-Based Image Coding via Inpainting 6.5.6 More Examples and Applications of TV Inpainting 6.6 Error Analysis for Image Inpainting 6.7 Inpainting Piecewise Smooth Images via Mumford and Shah 6.8 Image Inpainting via Euler's Elasticas and Curvatures 6.8.1 Inpainting Based on the Elastica Image Model 6.8.2 Inpainting via Mumford-Shah-Euler Image Model 6.9 Inpainting of Meyer's Texture 6.10 Image Inpainting with Missing Wavelet Coefficients 6.11 PDE Inpainting: Transport, Diffusion, and Navier-Stokes 6.11.1 Second Order Interpolation Models 6.11.2 A Third Order PDE Inpainting Model and Navier-Stokes . . . .

245 246 246 248 249 251 253 256 256 257 258 263 264 265 268 269 270 270 272 274 274 276 277 279 282 284 284 287 289 291 295 295 299

Contents 6.11.3 TV Inpainting Revisited: Anisotropic Diffusion 6.11.4 CDD Inpainting: Curvature Driven Diffusion 6.11.5 A Quasi-axiomatic Approach to Third Order Inpainting 6.12 Inpainting of Gibbs/Markov Random Fields 7

xi 301 302 303 307

Image Segmentation 309 7.1 Synthetic Images: Monoids of Occlusive Preimages 309 7.1.1 Introduction and Motivation 309 7.1.2 Monoids of Occlusive Preimages 310 7.1.3 Mimimal and Prime (or Atomic) Generators 315 7.2 Edges and Active Contours 318 7.2.1 Pixelwise Characterization of Edges: David Marr's Edges . . . .318 7.2.2 Edge-Regulated Data Models for Image Gray Values 320 7.2.3 Geometry-Regulated Prior Models for Edges 322 7.2.4 Active Contours: Combining Both Prior and Data Models . . . . 325 7.2.5 Curve Evolutions via Gradient Descent 327 7.2.6 F-Convergence Approximation of Active Contours 329 7.2.7 Region-Based Active Contours Driven by Gradients 331 7.2.8 Region-Based Active Contours Driven by Stochastic Features 332 7.3 Geman and Geman's Intensity-Edge Mixture Model 338 7.3.1 Topological Pixel Domains, Graphs, and Cliques 338 7.3.2 Edges as Hidden Markov Random Fields 339 7.3.3 Intensities as Edge-Regulated Markov Random Fields 342 7.3.4 Gibbs' Fields for Joint Bayesian Estimation of u and F 343 7.4 The Mumford-Shah Free-Boundary Segmentation Model 344 7.4.1 The Mumford-Shah Segmentation Model 344 7.4.2 Asymptotic M.-S. Model I: Sobolev Smoothing 345 7.4.3 Asymptotic M.-S. Model II: Piecewise Constant 347 7.4.4 Asymptotic M.-S. Model III: Geodesic Active Contours 351 7.4.5 Nonuniqueness of M.-S. Segmentation: A 1-D Example 355 7.4.6 Existence of M.-S. Segmentation 355 7.4.7 How to Segment Sierpinski Islands 359 7.4.8 Hidden Symmetries of M.-S. Segmentation 362 7.4.9 Computational Method I: F-Convergence Approximation . . . .364 7.4.10 Computational Method II: Level-Set Method 366 7.5 Multichannel Logical Segmentation 369

Bibliography

373

Index

393

This page intentionally left blank

List of Figures 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

1.9 1.10 1.11

1.12 1.13 1.14 1.15 1.16

An ideal image: noiseless, complete, and in good contrast A low-contrast version of the ideal image Degradation by additive Gaussian noise (Chapter 4) Salt-and-pepper noise with 10% spatial density (Chapter 4) Degradation by out-of-focus blur: the digital camera focuses on a fingertip about 1 inch away while the target scene is about 1 foot away (Chapter 5) Degradation by motion blur due to horizontal hand jittering during a single exposure (Chapter 5) 150 8-by-8 packets are randomly lost during transmission (Chapter 6). The goal of error concealment (or more generally, inpainting) is to develop models and algorithms that can automatically fill in the blanks [24, 67, 166]. Such cartoonish segmentation seems trivial to human vision but still remains the most fundamental and challenging problem in image processing and low-level computer vision (Chapter 7). (The background segment Wo is not shown explicitly here.) A binary set A and a structure element S (for dilation and erosion). ... Dilation D s ( A ) (left) and erosion E s ( A ) (right): dilation closes up small holes or gaps, while erosion opens them up A digital image and its discrete Fourier transform. The example reveals a couple of salient features of Fourier image analysis: (1) most highamplitude coefficients concentrate on the low-frequency band; (2) dominant directional information in the original image is easily recognizable in the Fourier domain; and (3) the coefficients, however, decay slowly for Heaviside-type directional edges (i.e., a jump line with distinct constant values along its two shoulders) An example of a (mother) wavelet by Daubechies' design [96]. Localization and oscillation are characteristic to all wavelets A noisy 1-D signal and its optimal restoration according to (1.7) A trefoil evolves under mean-curvature motion (1.11) or (1.12) Different approaches are intrinsically connected: an example Organization of the book

XIII

7 7 8 9 10

10 11

11 13 14

15 16 19 21 22 24

xiv

List of Figures

2.1 2.2 2.3 2.4

2.5 2.6 2.7

2.8

2.9

2.10

2.11 2.12 2.13

2.14 2.15 2.16 2.17

A planar curve: tangent t, normal n, and curvature K Second order local geometry: surface normal N and two perpendicular principle curvature circles A surface patch z — h(u, v] — cos(u)cos(u) (top), its mean curvature field H (lower left), and Gauss curvature field K (lower right) A set E (left) and a covering A (right). Covering elements like the huge box with dashed borders (right) fail to faithfully capture small-scale details. The definition of Hausdorff measures forces covering scales to tend to zero The Hausdorff dimension dimH(E) of a set E is the critical d, above which E appears too "thin," while below it appears too "fat." Three 1-D images with TV[f] = TV[g] = TV[h] = 2 An example of L 1 -lower semicontinuity. The sequence (un) of 1-D images on [0, 1 ] converges to u = O i n L ' since ||u,n+1 — u||L1 2-n. Notice that TV(u) = 0 while TV(u n ) = 2, which is consistent with the property of lower semicontinuity: TV(u) lim inf n TV(u n ). In particular, strict inequality is indeed realizable Geometric meaning of TV: the co-area formula. For smooth images, TV[w] is to sum up the lengths of all the level curves, weighted by the Lebesgue element dk. Plotted here is a discrete approximation: TV[«] ~ £„ length(« = X,,) AA Gibbs' CE: higher temperature T corresponds to smaller /3 and more uniform distribution; on the other hand, when T — 0, the system exclusively remains at the ground state (leading to superfluids or superconductance in physics) Gibbs' entropies for two imaginary 4-state systems (with K set to 1). Entropy generally measures the degrees of freedom of a system. A lower entropy therefore means that the target system is more restrained. (This observation led Shannon [272] to define negative entropies as information metrics, since less randomness implies more information.) An example of linear anisotropic diffusion by (2.38) with diagonal diffusivity matrix D = diag(D v , D v ) and Dy : D v = 10 : 1. The image thus diffuses much faster along the vertical y-direction The median of a 5-component segment from an imaginary 1-D signal (*„). An example of using the median filter (with 7 x 7 centered square window) to denoise an image with severe salt-and-pepper noise. Notice the outstanding feature of median filtering: the edges in the restored image are not blurry Two examples of Marr's wavelets ("Mexican hats") as in (2.50) A generic approach to designing the window template w(x) via symmetrically constructing the profile of its square w2(x) MRA as (Hilbert) space decompositions: the finer resolution space V-± is decomposed to the detail (or wavelet) space W\ and coarser resolution space V\. The same process applies to V\ and all other V/s [290] A pair of compactly supported scaling function 4>(x) and mother wavelet iAU) by Daubechies' design [96]

32 38 40

45 45 47

49

54

59

60 68 70

71 75 79 84 87

List of Figures 2.18 3.1

3.2 3.3

3.4

3.5 3.6

3.7

4.1 4.2

4.3 4.4

4.5

Fast wavelet transform via filter banks: the two-channel analysis (or decomposition) and synthesis (or reconstruction) banks

xv

88

Images as distributions or generalized functions. Test functions then model various biological or digital sensors, such as retinal photoreceptors of human vision or coupled charge devices in CCD cameras 94 Three Haar mother wavelets in two dimensions via tensor products. . . . 103 Besov norms in B®(LpYs measure the strength of signals in the spacescale plane: Lp for intrascale variations, Lq for inter- or cross-scale variations (in terms of dX = dh/h), while h~a for comparison with Holder continuity 108 Ingredients of Markov random fields by examples: a neighborhood Na (left), two doubleton cliques C e C (middle), and locality of conditional inference p(ua \ «n\ a ) = p(ua \ uNu) (right) 120 An example of branches LX,AA.'S with A A, = h and A. = nh. The branch L^,h contains two leaflets 132 The(xn, bn, gn) -representation of a (compactly supported) piecewise smooth signal u. On smooth regions where the signal u varies slowly, gn is often small and only a few bits suffice to code them 135 Poisson representation of a piecewise smooth signal u. Shown here is only a single representative interval. The signal u is represented by the two boundary values «+ and u ~+, and its second order derivative / = u" (corresponding to the source distribution in electromagnetism). Reconstructing u on the interval then amounts to solving the Poisson equation (3.62). The advantage of such representation is that / is often small (just like wavelet coefficients) for smooth signals and demands fewer bits compared with u 136 A sample 2-D Brownian path W(t) with 0 < / < 4, and the reference circle with radius 2 = V4, the standard deviation of W (4) A sample random signal generated by the SDE dX = 2Xdt + XdW with X(Q) = 1. The smooth dashed curve denotes the mean curve x(t) = E X ( t ) = e2'. The random signal is clearly not stationary since both the means and variances of X ( t ) evolve Hard thresholding T\(t) and soft shrinkage 5^(0 The shrinkage operator a* = Sff(a0) achieves minimum shrinkage (Theorem 4.4) among all the estimators a that satisfy the uniform shinkage condition (4.25) Two singleton error functions e\(t\ao) with A. = JJL = 1 and ao = 1,2. Notice that the optimal estimators (i.e., the valleys) a = 0, 1 are precisely the shrinkages Sa(ao) with a = fji/'k = 1

152

155 161

163

165

xvi

List of Figures

4.6

4.7 4.8

4.9 4.10 4.11

4.12 4.13

5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9

6.1 6.2

6.3

An example of shrinkage-based image denoising, with Gaussian white noise level a =0.1, and four levels of wavelet decompositions based on one of Daubechies' orthonormal wavelets [96]. The three bottom panels demonstrate the effect of the threshold parameter A on the shrinkage operator Sx '• a too large X causes overshrinkage and oversmoothing, while a too small one leads to the opposite 171 Same example as in Figure 4.6 with increased Gaussian noise level a =0.2 172 Denoising effects of moving means vs. moving medians. Both employ the same symmetric (moving) window of four neighbors on either side. As expected, the median filter preserves sharp edges better than the mean filter. (The dashed curve denotes the ideal clean signal.) 175 An example of Rudin, Osher, and Fatemi's TV denoising 179 A target pixel O and its neighbors 181 An example of applying Perona and Malik's anisotropic diffusion to a noisy image. The nonlinear diffusivity coefficient D(\ VM|) for this example is the one in (4.93) and (4.94). Notice the remarkable feature of edge preservation, as contrast to linear heat diffusions that invariably blur sharp edges while removing noises 194 Two examples of salt-and-pepper noises and the denoising performance of median filtering. As consistent with theoretical analysis, median filtering works very efficiently for low spatial densities but poorly for high ones. 202 An example of applying the TV inpainting technique (4.103) (also see Chan and Shen [67] or Chapter 6) to the denoising of high-density saltand-pepper noises 203 An example of motion blur (due to camera jittering) Motion blur of the image of a point source An example of out-of-focus blur Geometric optics of out-of-focus imaging An example of Wiener filtering for both denoising and deblurring Deblurring an out-of-focus image (with known PSF) Deblurring an image blurred by horizontal hand jittering Another example of restoring an image blurred by motion A computational example for the double-BV blind deblurring model (by Chan and Wong [77]). Left: the blurry image; Right: the deblurred image The goal of inpainting is to reconstruct the ideal image u on the entire domain based on the incomplete and often degraded data u° available outside the missing (or inpainting) domain D The effect of local and global pattern recognition on image inpainting: the black color seems to be a reasonable inpainting solution in the left panel, while for the right chessboard pattern, the white color becomes more plausible (Chan and Shen [67]) The effect of aspect ratios on image inpainting [67]

209 210 212 213 219 225 226 226 243

246

257 258

List of Figures

6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17 6.18 6.19 6.20 6.21 7.1

7.2 7.3 7.4 7.5 7.6

xvii

Harmonic inpainting of a smooth image (u = r = y/x 2 + j2) and an ideal step edge 263 Can the TV inpainting model explain Kanizsa's Entangled Man? 272 The model for Kanizsa's Entangled Man 273 Inpainting a noisy edge 277 Inpainting occluded bars 277 Inpainting a noisy face 278 Inpainting for text removal 278 Digital harmonic zoom-in and TV zoom-in (Chan and Shen [67]) 279 Edge decoding by TV inpainting (Chan and Shen [67]; also see the recent pattern-theoretic approach by Guo, Zhu, and Wu [148]) 280 The aspect ratios of inpainting domains influence the performance of inpainting schemes (Chan and Kang [59]) 281 Inpainting based on the F-convergence approximation (6.42) and its associated elliptic system (6.47) 284 Text erasing by inpainting based on the Mumford-Shah image model. . .284 Two examples of elastica inpainting, as compared with TV inpainting. In the case of large aspect ratios [68], the TV inpainting model fails to comply to the connectivity principle 286 Inpainting based on the Mumford-Shah-Euler image model can satisfactorily restore a smooth edge as expected 289 An example of inpainting a noisy and incomplete BV image with missing wavelet components using the model (6.64) (Chan, Shen, and Zhou [74]). 294 Another example of inpainting a noisy and incomplete BV image with missing wavelet components using the model (6.64) (Chan, Shen, and Zhou [74]) 294 TV fails to realize the connectivity principle in inpainting problems with large scales (or aspect ratios) 301 An example of CDD inpainting for scratch removal (Chan and Shen [68]) 302 Role of occlusion in visual perception. Left panel: a gray vase against a bright background or two opposite human faces against a gray background? This is a typical illusion problem caused by the lack of occlusion cues; Right panel (a trefoil): the role of occlusion in knot theory [253]. An example of two preimages aY and b$ and their occlusion ay H bs. Complex image patterns in the real world often originate from simple objects via occlusion Left: a diskette preimage; Right: a planar preimage On an edge pixel x (with y ( x ) — 1), p — \ Vu\(x) is more likely to be large according to (7.17) An example of geodesic active contouring: the given image and three different stages of the contour evolution by (7.32) An example of Chan and Vese's Active Contour model (7.52) without gradients (i.e., region based and driven by mean fields): three different evolution stages for a sectional brain image (from Chan and Vese [75]).

.310 312 316 322 329 . 337

xviii 7.7

7.8 7.9

7.10

7.11 7.12

7.13

7.14 7.15 7.16 7.17

List of Figures An example of texture segmentation by combining Chan and Vese's Active Contour model (7.52) and Gabor filters (7.53) (shown at increasing times; see [65, 260]) A pixel lattice £1 = [a,b, ...} and an edge lattice E = [a, @ , . . . } Six binary edgel configurations on a maximal edge clique and their empirical potentials V^(F)'s assigned by Geman and Geman [130] (where A is a constant) The total curvatures (in the interior of the pixel squares) of closed tours around the six edge configurations qualitatively match Geman and Geman's empirical choice of potentials in Figure 7.9 Visualization of the example: the set S consisting of the circles at all scale levels has the weird property H* (S) < oo but ~Hl (S) = oo, since S = Q. Left panel: the complement G of the classical Sierpinski gasket (Mandelbrot [206]) consists of all the shaded open triangles (excluding their sides). Right panel: Sierpinski Islands G(p) is the open set obtained from G by contracting each triangle of G with some suitable rate controlled by some p e (0, 3/2) (see the text). And the target image up discussed in this section is the indicator function of G(p) An example of F-convergence-based approximation to the M.-S. segmentation model (7.97): the optimal image estimation u = u £ ( x ) and its associated edge "canyon" function z — z£(x). Computation has been based on the alternating minimization scheme (7.98) An example of M.-S. segmentation using the level-set algorithm by Chan and Vese [75, 76] A synthetic example of an object in two different channels. Notice that the lower left corner of A\ and the upper corner of A 2 are missing Different logical combinations for the sample image: the union, the intersection, and the differentiation Region-based logical model on a medical image. In the first channel A \, the noisy image has a "brain tumor," while channel AI does not. The goal is to spot the tumor that is in channel A\, but not in AT, i.e., the differentiation A\ D —>A2. In the right column, we observe that the tumor has been successfully captured

338 340

341

342 357

360

366 369 369 370

372

Preface No time in human history has ever witnessed such explosive influence and impact of image processing on modern society, sciences, and technologies. From nanotechnologies, astronomy, medicine, vision psychology, remote sensoring, security screening, and the entertainment industry to the digital communication technologies, images have helped mankind to see objects in various environments and scales, to sense and communicate distinct spatial or temporal patterns of the physical world, as well as to make optimal decisions and take right actions. Image processing and understanding are therefore turning into a critical component in contemporary sciences and technologies with many important applications. As a branch of signal processing, image processing has traditionally been built upon the machinery of Fourier and spectral analysis. In the past few decades, there have emerged numerous novel competing methods and tools for successful image processing. They include, for example, stochastic approaches based upon Gibbs/Markov random fields and Bayesian inference theory, variational methods incorporating various geometric regularities, linear or nonlinear partial differential equations, as well as applied harmonic analysis centered around wavelets. These diversified approaches are apparently distinct but in fact intrinsically connected. Each method excels from certain interesting angles or levels of approximations but is also inevitably subject to its limitations and applicabilities. On the other hand, at some deeper levels, they share common grounds and roots, from which more efficient hybrid tools or methods can be developed. This highlights the necessity of integrating this diversity of approaches. The present book takes a concerted step towards this integration goal by synergistically covering all the aforementioned modern image processing approaches. We strive to reveal the few key common threads connecting all these major methodologies in contemporary image processing and analysis, as well as to highlight some emergent integration efforts that have been proven very successful and enlightening. However, we emphasize that we have made no attempt to be comprehensive in covering each subarea. In addition to the efforts of organizing the vast contemporary literature into a coherent logical structure, the present book also provides some in-depth analysis for those relatively newer areas. Since very few books have attempted this integrative approach, we hope ours will fill a need in the field. Let u denote an observed image function, and T an image processor, which can be either deterministic or stochastic, as well as linear or nonlinear. Then a typical image

XIX

xx

Preface

processing problem can be expressed by the flow chart

where F represents important image or visual features of interest. In the present book, we explore all three key aspects of image processing and analysis. • Modeling: What are the suitable mathematical models for u and Tl What are the fundamental principles governing the constructions of such models? What are the key features that have to be properly respected and incorporated? • Model Analysis: Are the two models for u and T compatible? Is T stable and robust to noises or general perturbations? Does F = T[u] exist, and if so, is it unique? What are the fine properties or structures of the solutions? In many applications, image processors are often formulated as inverse problem solvers, and as a result, issues like stability, existence, and uniqueness become very important. • Computation and Simulation: How can the models be efficiently computed or simulated? Which numerical regularization techniques should be introduced to ensure stability and convergence? And how should the targeted entities be properly represented? This view governs the structure and organization of the entire book. The first chapter briefly summarizes the emerging novel field of imaging science, as well as outlines the main tasks and topics of the book. In the next two chapters, we introduce and analyze several universal modern ways for image modeling and representation (for u), which include wavelets, random fields, level sets, etc. Based on this foundation, we then in the subsequent four chapters develop and analyze four specific and significant processing models (for T) including image denoising, image deblurring, inpainting or image interpolation, and image segmentation. Embedded within various image processing models are their computational algorithms, numerical examples, or typical applications. As the whole spectra of image processing spread so vastly, in this book we can only focus on several most representative problems which emerge frequently from applications. In terms of computer vision and artificial intelligence, these are often loosely categorized as low-level vision problems. We do not intend to cover high-level vision problems which often involve pattern learning, identification, and representation. We are enormously grateful to Linda Thiel, Alexa Epstein, Kathleen LeBlanc, Michelle Montgomery, David Riegelhaupt, and Sara Murphy of the SI AM Publisher for their constant encouragement and care throughout the project. It has been such a wonderful experience of planning, communication, and envisaging. We also owe profound gratitude to the following colleagues whose published works and personal discussions have greatly influenced and shaped the contents and structures of the current book (in alphabetical order): Antonin Chambolle, Ron Coifman, Ingrid Daubechies, Rachid Deriche, Ron DeVore, David Donoho, Stu Geman, Brad Lucier, Jitendra Malik, Yves Meyer, Jean-Michel Morel, David Mumford, Stan Osher, Pietro Perona, Guillermo Sapiro, Jayant Shah, James Sethian, Harry Shum, Steve Smale, Gilbert Strang, Curt Vogel, Yingnian Wu, Alan Yuille, and Song-Chun Zhu, and the list further expands. The book would be impossible without the generous support and help of many friends: Doug Arnold, Andrea Bertozzi, Carme Calderon, Charles Chui, Babette Dalton, Bjorn

Preface

xxi

Enquist, Bob Gulliver, Xiaoming Huo, Kate Houser, Yoon-Mo Jung, Dan Kersten, Paul Schrater, Mitch Luskin, Riccardo March, Willard Miller, Peter Olver, Hans Othmer, Fadil Santosa, Zuowei Shen, Chi-Wang Shu, Luminita Vese, Kevin Vixie, Tony Yezzi, Hong-Kai Zhao, and Ding-Xuan Zhou. Tony Chan would like to personally thank his colleagues for numerous inspirations and interesting discussions: Emmanuel Candes, Raymond Chan, Li-Tien Cheng, Ron Fedkiw, Mark Green, Michael Ng, Stefano Soatto, Xue-Cheng Tai, Richard Tsai, and Wing Wong. Tony Chan is also deeply grateful to the following students and postdoctoral fellows for the privilege of working with them in this new and exciting field: Peter Blomgren, Jamylle Carter, Selim Esedoglu, Sung-Ha Kang, Mark Moelich, Pep Mulct, Fred Park, Berta Sandberg, Jackie Shen, Bing Song, David Strong, Justin Wan, Yalin Wang, Luminita Vese, Chiu-Kwong Wong, Andy Yip, Hao-Min Zhou, and Wei Zhu, Jackie Shen wishes to personally thank Professors Gilbert Strang, Tony Chan, Stan Osher, David Mumford, and Stu Geman for their profound influence on his scholastic growth in the field, as well as numerous personal friends for warming and shining up each ordinary day during the project: Tianxi Cai, Shanhui Fan, Chuan He, Huiqiang Jiang, Ming Li, Tian-Jun Li, William Li, Chun Liu, Hailiang Liu, Huazhang (Andy) Luo, Conan Leung, Mila Nikolova, Jiaping Wang, Chao Xu, Jianling Yuan, Wen Zhang, and Qiang Zhu. We are very grateful to Jean-Michel Morel and Kevin Vixie specifically for their insightful and constructive comments on an early version of the manuscript for this book, which have significantly improved the quality. During this book project, the authors are enormously grateful for the support from the National Science Foundations (NSF-USA), the Office of Naval Research (ONR-USA), as well as the National Institute of Health (NIH-USA). In particular, Tony Chan would like to thank Wen Master (at ONR) for the continuous support to this novel area in applied and computational mathematics. We also acknowledge the tremendous benefits from the participation of numerous imaging sciences related workshops at the Institute of Pure and Applied Mathematics (IPAM) at UCLA, the Institute of Mathematics and Its Applications (IMA) at the University of Minnesota, the Institute of Mathematical Sciences (IMS) at the National University of Singapore, the Mathematical Sciences Research Institute (MSRI) at Berkeley, as well as the Center of Mathematical Sciences (CMS) at Zhejiang University, China. Finally, this book project, like all the others in our life, is an intellectual product under numerous constraints, including our busy working schedules and many other scholastic duties. Its contents and structures as presented herein are therefore only optimal subject to such inevitable conditions. All errata and suggestions for improvements will be received gratefully by the authors. This book is absolutely impossible without the pioneering works of numerous insightful mathematicians and computer scientists and engineers. It would be our great pleasure to see that the book can faithfully reflect many major aspects of contemporary image analysis and processing. But unintentional biases are inevitable due to the limited views and experiences of the authors, and we are happy to hear any criticisms from our dear readers.

This page intentionally left blank

Chapter 1

Introduction

The last few decades have witnessed a new era of imaging sciences. From satellite imaging and X-ray imaging to modern computer aided CT, MRI, and PET imaging in medical sciences, mankind's naked vision is no longer restricted by the conventional notions of space, time, scales, and visibility. Consequently, just as Mother Nature has achieved for biological vision, there is tremendous need for developing the virtual brain structures for all these imaging areas, i.e., models and algorithms for making these mechanical or artificial vision systems more efficient, accurate, and stable, as well as artificial intelligence that can successfully process, communicate, and interpret the resultant images. The current book attempts to offer a well-organized view of this rapidly evolving field. This introductory chapter first gives an overview of this exciting emerging field, and then introduces the scopes and structures of the current book. A detailed guidance on how to efficiently read and use the book is also speculated for different groups of potential readers.

1.1

Dawning of the Era of Imaging Sciences

Imaging sciences consist of three relatively independent components: image acquisition, image processing, and image interpretation. This book mainly focuses on the second subject, though we shall first give a quick overview of all three.

1.1.1

Image Acquisition

Image acquisition (orformation) studies the physical mechanisms by which particular types of imaging devices generate image observations, and it also investigates the associated mathematical models and algorithms employed by the computers integrated into such imaging devices. We now briefly describe several important image acquisition problems. 1

2

Chapter 1. Introduction

Human Vision The human vision system is Mother Nature's biological device for carrying out intelligent visual perception. The imaging formation system mainly consists of optical peripherals including the cornea, pupil, iris, and lens, as well as two different types of photoreceptors of the retinas, cones and rods, which are responsible for capturing colors and gray value intensities separately [168, 227]. The left and right retinas make up the frontal end of a complex vision system that transforms light or photons into electrochemical signals and prepares for further image processing and interpretation in the inner visual cortices through optical neuron fibers. The perceptible light spectra for average naked vision are about between 400 and 700nm (nanometer, or 10~9 meter) in terms of wavelengths. On the other hand, human vision achieves astounding sensitivity to a large dynamic range of light intensities due to its delicate adaptive mechanism described by the celebrated Weber's law in vision psychology, psychophysics, and retinal physiology [121, 168, 275, 315]. Night Vision and Thermal Imaging Wavelengths longer than 700nm belong to the infrared (IR) range. Mother Nature forces human beings to sleep and replenish during the night by making IR light imperceptible to normal naked eyes. However, assisted by proper signal amplification and enhancement devices (e.g., night vision goggles), human beings are capable of perceiving IR images ranging from just above 700nm to a few microns (10~6 meter). Furthermore, the majority of light waves in the IR spectra (ranging from 3 to 30 microns) are not reflected by objects as in the visible spectra but emitted by objects due to their intrinsic thermal and statistical activities at atomic levels. For such weak photons, thermal imaging devices can still collect, enhance, and process them, producing the socalled thermograms that map the detailed distribution of temperatures of a target scene. Assisted with proper signal processing hardware and software, one can then visualize the otherwise imperceptible images. Radar Imaging and SAR Radar is a device that employs RAdio waves to make Detection And Ranging of targets. The signals are in the microwave range between 1cm and 1m wavelengths and polarized. Radar antennas emit these polarized radio wave pulses and receive their echoes reflected from objects. The waiting time between sending and receiving, as well as the strength of the echoes or backscalters, reveal important information such as ranges, local configuration, and material differentiation of the targets. Applied to imaging, radar becomes a powerful tool for remote sensoring. To acquire ground images, a radar device is usually attached to a flight vehicle such as a space shuttle and continuously sends, receives, and processes signalling microwave pulses. The radar image is then constructed along its flight path. The resolution of the acquired image is inversely proportional to the size of the radar; that is, larger antennas result in higher resolution images. This hardware restriction on image resolution can, however, be effectively loosened using clever models, algorithms, and software that can extract extra spatial information

1.1. Dawning of the Era of Imaging Sciences

3

from the correlation among the echoes. This amounts to increasing the effective size of the real antenna being employed, and thus the resultant imaging technique is called synthetic aperture radar (SAR) imaging. Aperture refers to the opening area of an antenna that collects echoed pulses. Ultrasound Imaging Like ultraviolet light, ultrasound refers to acoustic waves whose frequencies are higher than the audible range of normal human ears, i.e., about between 20Hz and 20KHz. Medical ultrasound is often above 1MHz and below 20MHz. The different responses of the various tissues and structures to ultrasound pulses form the physiological foundation for ultrasound imaging. By measuring the time intervals between emitting and receiving sound pulses, as well as their strengths, the range and local topological information of target physiological entities can be inferred and displayed using computers. Furthermore, using Doppler shifts, ultrasound imaging can be employed to study organs in real-time motion (e.g., heart valves and blood vessels). Recall that the Doppler shift refers to the phenomenon that the reflected echoes from a moving object have higher or lower frequencies than the incident waves depending on whether the motion is toward or away from the wave source. Computed (Axial) Tomography (CAT or CT) Ever since the legendary discovery of X-rays by the German physicist Wilhelm Rontgen in 1895, who was then awarded the first Nobel Prize of physics in 1901, X-rays have become the most useful probing signals in medical imaging and diagnosis. Unlike the visible spectra, IR range for night vision or thermal imaging and the microwaves for radar, X-rays have wavelengths in the nanometer (10~ 9 ) or picometer (10~ 12 ) scales, which are energetic enough to penetrate most materials including biological tissues. This makes X-rays ideal for noninvasive medical probing and imaging. Different tissues and organs in the human body absorb X-rays at different rates. Denser structures like the bones absorb more X-rays than softer tissues. Thus conventional X-ray imaging is very similar to ordinary photographing, with a film directly exposed to the penetrated X-rays, and the different intensity patterns on the film differentiate different organs and structures. Computed (axial) tomography, often referred to as CAT or CT scanning, improves the above conventional X-ray imaging method by more delicate hardware design as well as computer software for reconstructing the final images. CT can be called rotational and localized X-ray imaging technology since it consists of a localized X-ray source and a detector which are put opposite to each other and rotate around the human body's longitudinal central axis. A two-dimensional (2-D) slice image is then synthesized from the collected X-ray signals along all the directions. Since the mid 1970s, CT scanning has witnessed broad applications in medical imaging and diagnosis, thanks to its fast speed and patient-friendly nature, as well as its unique capability in imaging a wide variety of soft tissues, hard bones, blood vessels, and so on.

4

Chapter 1. Introduction

Magnetic Resonance Imaging (MRI) The 2003 Nobel Prize in Physiology or Medicine was awarded to Paul Lauterbur and Peter Mansfield for their contributions on magnetic resonance imaging (MRI). The MRI imaging technology has been built upon the remarkable physical theory of nuclear magnetic resonance (NMR), discovered and developed by two American physicists, Felix Bloch and Edward Purcell, who won the Nobel Prize in 1952. A nucleon such as a proton or neutron possesses spin, which comes at a multiple of 1/2 with plus or minus differentiation. Paired opposite spins can effectively cancel each other out, leading to no evident observables. In an external magnetic field, unpaired or net spins tend to align with the orientation of the field, and statistical mechanics describes the distribution of two opposite alignments. The two alignments in different energy levels can transit to each other by either absorbing or releasing a photon energy proportional to the external magnetic field. For hydrogen atoms in clinical applications, for example, the frequencies of such photons are in the radio frequency (RF) range, i.e., MHz. Simply speaking, MRI imaging works as follows. First the atomic spins (e.g., hydrogen nuclei) of the human body are aligned by applying a strong external magnetic field. Then high frequency RF pulses are emitted into a slice plane perpendicular to the external field. By absorbing the photons in the RF waves, the aligned spins can transit from the lower energy alignment to the higher (or excited) one. After the RF waves are turned off, the excited nucleus ensemble starts to gradually resume to its original distribution, a process known as relaxation. It is the relaxation time that reveals or differentiates different tissues or body structures. MRI imaging hardware and software measure the distribution of relaxation times and construct the associated spatial images. MRI technology has many new important developments, including functional MRI (fMRI) for brain mapping [36] and diffusion tensor imaging (DTI) for studying neural fibers [314, 319]. Computer Graphics and Synthesized Images All the above image acquisition processes concern how to design hardware and combine it with powerful software to construct targeted images in the real three-dimensional (3-D) world. Computer graphics, on the other hand, is a purely so/rvvare-based image acquisition approach, or rather, an art that creates images of imaginary 3-D scenes. As a result, computer graphics has become increasingly important in the movie industry (e.g., popular movies like Shrek, Finding Nemo, Titanic, The Matrix, and The Lord of the Rings, to name only a few). Computer graphics can be formulated as follows. With an imaginary scene in mind, a designer first creates the 3-D coordinates G for shapes and geometric information of imaginary objects, as well as their optical properties such as the degree of transparency and surface reflectance R. The designer then designs a light source / for illuminance, which could be a point source or scattered sources such as in a normal daytime environment, and also picks a viewpoint 9 of an imaginary observer. The final imaginary image u is then synthesized based on some proper optical relation F:

1.1. Dawning of the Era of Imaging Sciences

5

If motion is involved, such as in movie production, time t has to be included as well:

which then becomes a dynamic process of image sequence generation. Human dreaming can be considered as Mother Nature's subconscious art of graphics, not by computers but by human visual neurons and cortices. In addition to all the aforementioned optical, electromagnetic, or nuclear imaging areas, there also exist other imaging techniques such as underwater or geophysical acoustic imaging (see, e.g., Papanicolaou [244]).

1.1.2

Image Processing

Any image processor can be abstractly formulated by an input-output system:

The input UQ denotes an already acquired image, which could be degraded due to either poor imaging conditions or problems during storage and communication. Mathematically speaking, an image processor T could be any linear or nonlinear operator that operates on the inputs and produces targeted features or patterns F. The design and development of image processors are often driven by specific applications or tasks. For example, in situations where observed images are degraded, one wishes to restore or enhance them to obtain high-quality images, based upon which other important visual decisions can be made more robustly and safely. Tumor detection from lowly contrasted CT images is such an example, while detecting celestial details such as binary stars from a blurry Hubble telescope image is another one. The main challenge for designing and computing T is that most image processing problems are ill-posed inverse problems. That is, there often exists a. forward problem that generates an image observation UQ from the targeted feature variable F if F is known or given. Most image processing problems are, however, to recover or detect F from UQ, and are therefore inverse problems. Ill-posedness often results from the nonuniqueness or instability of direct inversion. The output F = T[UQ] can be generally any collection of features that are visually meaningful. For example, F could denote the location of corners, the layout of object boundaries, the disjoint regions associated with distinct objects, or the relative order or depth of different objects in a scene. When the input UQ is a degraded image, the output F could also simply be its enhanced ideal version, often denoted by u, The input UQ can contain more than one image, such as in surveillance cameras, video or movie processing, and continuous medical monitoring of a single patient. For video or movie processing in particular, the time variable can lead to new features such as velocity distribution (or optical flows), dynamic evolution of the shapes and areas of interesting targets (or target tracking), and image registration. This book mainly focuses on image processing, and more specifically, on four classes of processing tasks which frequently arise from numerous applications: denoising, deblurring, inpainting or image interpolation, and segmentation. We shall further elaborate on them in the next major section.

Chapter 1. Introduction

6

By far it seems that image processing can be separated apart from image acquisition. This is, however, only partially true and is dependent on the real tasks at hand. A doctor who tries to decipher the image patterns of a patient certainly does not need to know nuclear physics or the physics details of MRI. But knowledge of an image acquisition process can often help develop more powerful image processing models or algorithms. One also must be aware that many image acquisition devices have their built-in image processors in order to construct high-quality images during the image formation stage.

1.1.3

Image Interpretation and Visual Intelligence

The ultimate purpose of imaging sciences is to be able to "see," monitor, and interpret the targeted portion of the world being imaged, whether it is the surface of Mars, the newly discovered tenth planet of the solar system, Sedna, or an MRI image slice of the brain. Therefore image interpretation is of fundamental importance for imaging sciences and is intrinsically connected to vision-based artificial intelligence. Classical image processing is often loosely identified with low-level vision, while image interpretation generally belongs to high-level vision. An ideal artificial vision system shall be able to perform all the major functions of the human vision system, including depth perception, motion estimation, and object differentiation and recognition. A combination of image processing and image interpretation is often necessary for these tasks. Image interpretation is partially the inverse problem of image acquisition. The latter studies how to form 2-D images from 3-D structures and patterns, which are arranged or positioned with relative orders, depths, opaqueness or transparency, etc. Image interpretation, on the other hand, attempts to reconstruct or interpret the 3-D world from 2-D images. For example, a well-trained radiologist can identify some abnormal image patterns in an MRI image with tumor tissues in a real 3-D body. Mathematically speaking, the foundation of image interpretation is a rapidly growing field called pattern theory, which mainly covers pattern modeling and representation, learning theory, and pattern detection and identification. This developing field in mathematics has been pioneered by, for example, S. Geman and D. Geman [130], Grenander [143], Mumford [224], Poggio and Smale [252], and Vapnik [307].

1.2

Image Processing by Examples

This section explains some typical image processing tasks through tangible examples. Some will be explored in later chapters with great details, while others can be found in more classical image processing books [140, 194]. Throughout this section, the working example is the ideal image displayed in Figure 1.1.

1.2.1

Image Contrast Enhancement

Figure 1.2 shows an image with very low contrast, which makes it difficult for human vision to clearly perceive important visual features. Assume that u is a gray-scale image in the range of 0 (black) and 1 (white). Let M denote the maximum value of u and m the minimum value. Then low contrast often means

1.2. Image Processing by Examples

7

Figure 1.1. An ideal image: noiseless, complete, and in good contrast. that M and m are not separated far enough. It may appear that a simple linear transform like can renormalize the range and suffice to enhance the image. But imagine the following scenario. Let w be a normal image with natural contrast and w(x) e [0, 1] for any pixel x. Suppose it is distorted by some nonlinear process to a low-contrast image u by

Figure 1.2. A low-contrast version of the ideal image.

Chapter 1. Introduction

8

Then for u, 0.4 < m < M < 0.5, and indeed the contrast is quite low. Furthermore, due to the nonlinearity of the distortion, the simple linear transform u —> i> proposed above cannot work faithfully to approximate the original image w. In real scenarios, there are infinitely many transforms like (1.1) that could lead to lowcontrast images. Further complexities may arise when the transforms are pixel dependent. It thus starts to become an interesting and nontrivial problem: without any specific knowledge on low-contrast transforms like (1.1), how could one design a scheme that can optimally restore the contrast, and perhaps even more importantly, what are the proper criteria for optimality? A very important solution to this problem in classical image processing is the histogram equalization technique, which can be read in many image processing books (e.g., [ 140,194]).

1.2.2

Image Denoising

Figure 1.3 shows a noisy image UQ\

where n(x) denotes additive Gaussian white noise and u the ideal "clean" image. To extract M from its noisy version MO is the denoising problem. The above additive Gaussian white noise model is popular for testing the performance of denoising models and algorithms. Another common noise model is the salt-and-pepper noise, which is nonlinear:

where u(x) e [0, 1] stands for the ideal clean image, c(x) a random field of independent 0-1 (with probabilities 50%:50%) binary random variables, and s(x) a random field of independent 0-1 binary random variables with p — Pmb(s(x) — 1) fixed for all .v e £i. c

Figure 1.3. Degradation by additive Gaussian noise (Chapter 4).

1.2. Image Processing by Examples

9

Figure 1.4. Salt-and-pepper noise with 10% spatial density (Chapter 4). and s are assumed to be independent. A sample of such noises can be simulated as follows. Take any two imaginary coins, one biased for s so that the chance of heads (s = 1) is p and the other unbiased for c. At each pixel x, first flip the s-coin. If s is a tail (s = 0), leave u(x) as it is and move to an unvisited new pixel. Otherwise (i.e., 5 = 1), further flip the c-coin, and change u(x) to the outcome of c. Repeat this game until each pixel has been visited exactly once. Then the resultant distorted image gives a sample of MO- Figure 1.4 shows an image degraded by the salt-and-pepper noise with p = 10%. Chapter 4 will further elaborate on noises and various approaches to suppress them.

1.2.3

Image Deblurring

Figure 1.5 shows a real blurry image out of focus. Blur is a common problem in satellite imaging, remote sensoring, and cosmological observations. For instance, before the Hubble telescope was launched into space in 1990, astronomic image observations had always been made by telescopes directly mounted on the surface of the Earth. As a result, images qualities were often inevitably blurred by atmospheric turbulence. Image blur occurs not only in ground telescope observations but also in other scenarios including fast motion or camera jittering (see Figure 1.6). Chapter 5 further explains in details different types of blurs, and models and algorithms to correct them, which is the problem of image deblurring.

1.2.4

Image Inpainting

Figure 1.7 shows a simulated image with 150 8-by-8 blocks randomly missing. The phenomenon often occurs in image communication where the packets are uncontrollably lost, for example, during wireless transmission. To restore the missing image information is in essence an interpolation problem, and is professionally called the error concealment problem by the communication community [166, 266].

10

Chapter 1. Introduction

Figure 1.5. Degradation by out-of-focus blur: the digital camera focuses on a fingertip about 1 inch away while the target scene is about 1 foot away (Chapter 5).

Figure 1.6. Degradation by motion blur due to horizontal hand jittering during a single exposure (Chapter 5).

In Chapter 6, readers see that error concealment is just one particular example of general image interpolation problems, which arise from numerous important applications including restoration of ancient paintings, movie production, image compression and construction, zooming and superresolution, and visual interpolation. The word "inpainting" is the artistic synonym for image interpolation and has been circulating among museum restoration artists for a long time [112, 312]. It was first introduced into digital image processing in the work of Bertalmio et al. [24] and later further popularized in many other works [23, 67, 61, 116].

1.2. Image Processing by Examples

11

Figure 1.7. 150 8-£>j-8 packets are randomly lost during transmission (Chapter 6). The goal of error concealment (or more generally, inpainting) is to develop models and algorithms that can automatically fill in the blanks [24, 67, 166].

1.2.5

Image Segmentation

Image segmentation crucially bridges low-level and high-level computer visions. Figure 1.8 shows an ideal image and its visually meaningful partition (or segmentation) of the entire image domain into different regions and boundaries. Almost trivial to human vision, image segmentation still remains one of the most challenging and most studied problems in image processing, image understanding, and artificial intelligence. A general image segmentation can be formulated as follows. Given on a 2-D domain £1 an image observation UQ, which could be degraded by noise or blur, find a visually

Figure 1.8. Such cartoonish segmentation seems trivial to human vision but still remains the most fundamental and challenging problem in image processing and low-level computer vision (Chapter 1). (The background segment QQ is not shown explicitly here.)

12

Chapter 1. Introduction

meaningful partitioning of the image domain,

for some image-dependent N, such that each component £2n (n > 1) visually corresponds to an "object" (besides the "background" patch S70)- Beneath this intuitive description of the segmentation problem lie the unavoidable complications on how to properly characterize the notions of (1) "visually meaningful," (2) image-dependent N, and even (3) "objects." These questions have been partially answered, e.g., in the celebrated works of Geman and Geman [130] and Mumford and Shah [226].

1.3

An Overview of Methodologies in Image Processing

In this section, we give an overview of some major approaches to image processing. It is unfair to say which one is necessarily superior to all the others. The efficiency and advantages of a particular methodology often depend on the concrete tasks at hand, as well as the classes and data structures of the images provided.

1.3.1

Morphological Approach

Since images capture objects, image processing naturally studies operations of objects. 2-D objects can be modelled as sets or domains of the 2-D plane R2 in a continuum setting or of the 2-D lattice Z2 in the discrete or digital setting. Equivalently, an object A can be identified with its binary characteristic function:

A morphological transform T is a map among objects,

which is often local in the sense that whether or not a pixel x belongs to B can be completely determined by the local behavior of A in the vicinity of x. Two most fundamental morphological transforms are the dilation operator D and erosion operator E, both depending on a structure element 5, or equivalent!y, a local neighborhood template. For example in the discrete setting, one could take

which is a 3-by-3 square template. The associated dilation and erosion operators are defined by [140]

Here 0 denotes the empty set and y + S = {y + s \ s e S}. It is evident that

1.3. An Overview of Methodologies in Image Processing

13

as long as (0, 0) e S, from which come the names "erosion" and "dilation." It is evident that both operators are monotonic:

Dilation and erosion provide convenient ways to define the "edge" or boundary of objects. For instance the dilative edge of an object can be defined by

while the erosive edge can be defined by

Edges are very crucial in visual inference and perception [163, 210]. Figures 1.9 and 1.10 show the effect of dilation and erosion on a binary object. In applications, two other morphological transforms generated from dilation and erosion play even more important roles. The first is the closing operator C$ — E$ o DS, i.e., with dilation followed by erosion. When applied to shapes, C$ can close small holes or gaps. And the second is the opening operator O$ = DS o E$, i.e., with dilation following after erosion. The net effect is that O$ can often erase narrow connectors and open up large holes or caves. The above binary morphological transforms can be naturally extended to general gray-level images via their level sets. Take the dilation operator DS for example. Let u be a general gray-level image defined on the lattice Z2. For each gray-level A., define the cumulative A-level set of u by

Then one defines another gray-level image v = D$(u], a dilated version of u, by

Figure 1.9. A binary set A and a structure element S (for dilation and erosion).

Chapter 1. Introduction

14

Figure 1.10. Dilation DS(A) (left) and erosion E S ( A ) (right): dilation closes up small holes or gaps, while erosion opens them up. By the monotonicity condition (1.2), FA (v) defined in this way indeed satisfies the necessary condition for cumulative level sets:

Therefore v is uniquely determined by the reconstruction formula

1.3.2

Fourier and Spectral Analysis

Fourier or spectral analysis has been one of the most powerful and favored tools in classical signal and image processing. If an image u is considered as a continuous function on the canonical rectangular domain £2 = (0, I) 2 , with periodic extension, the information of u can be completely encoded into its Fourier coefficients:

The completeness is in the sense of L 2 . Furthermore, in the digital setting when the image domain is actually a finite lattice Q — (0 : N — 1) x (0 : Af — 1) and the image u = (uj) = («;,./,) with j = (ji, 72) e Q a square matrix of data, one resorts to the discrete Fourier transform (DFT):

DFT is actually an orthonormal transform after a multiplicative scalar renormalization. Furthermore, DFT allows fast implementation known as Ihefast Fourier transform (FFT), which greatly facilitates numerous computational tasks in signal and image processing [237] (see the example in Figure 1.11). As in the spectral analysis of linear differential equations, Fourier transforms and their variations (e.g., cosine or sine transforms, and Hilbert transforms [237]) have been

1.3. An Overview of Methodologies in Image Processing

15

Figure 1.11. A digital image and its discrete Fourier transform. The example reveals a couple of salient features of Fourier image analysis: (1) most high-amplitude coefficients concentrate on the low-frequency band; (2) dominant directional information in the original image is easily recognizable in the Fourier domain; and (3) the coefficients, however, decay slowly for Heaviside-type directional edges (i.e., a jump line with distinct constant values along its two shoulders). universally applicable in many image processing tasks, including linear filtering and filter design, shift-invariant linear blurs, as well as classical image compression schemes such as the old JPEG protocol. In image analysis and processing, the Fourier approach has, however, been greatly challenged since the 1980s by another even more powerful tool—wavelet analysis.

1.3.3 Wavelet and Space-Scale Analysis Fourier transform mixes long-range spatial information, which makes it less ideal for handling localized visual features like edges. One extreme example is Dirac's delta image:

Then its Fourier coefficients are

That is, all the Fourier coefficients respond indiscriminately despite such a simple image. The harmonics are therefore inefficient in representing and coding localized image information [122]. On the other hand, cognitive and anatomic evidences in vision research have shown that human vision neurons are cleverly organized to be able to resolve localized features more efficiently [122, 152]. Wavelets in a certain sense ideally embody the idea of locality by resorting to localized bases, which are organized according to different scales or resolutions [96, 203]. Simply put, the wavelet representation of a given image u is defined by where i/ra's are wavelets, with the index set A resolving scale-adapted spatial locations.

Chapter 1. Introduction

16

Figure 1.12. An example of a (mother) wavelet by Daubechies' design [96]. Localization and oscillation are characteristic to all wavelets. The pioneering works of Nobel laureates Hubel and Wiesel [ 152] revealed the remarkable physiological fact that the simple and complex cells of human vision system behave like differentiators. Likewise, each wavelet \(fa possesses the similar differentiation property

implying that it remains dormant to featureless constant images. More generally, wavelets often satisfy the so-called vanishing-moment condition:

for some positive integer m. Coupled with spatial localization (i.e., rapid decay at infinities or being compactly supported [96, 204, 290]), this annihilation property immediately implies that the wavelet coefficients of a generic piecewise smooth image are mostly negligible except for those along important visual cues such as jumps or edges. Thus wavelets are efficient tools for adaptive image representation and data compression. Figure 1.12 displays a typical example of Daubechies' wavelets.

1.3.4

Stochastic Modeling

For images with notable stochastic nature, statistical methods become more suitable than deterministic approaches. There are two different sources leading to the statistical nature of a given image observation UQ. (a) UQ is the composition of an ideal (and often deterministic) image u with some random effect X: where the function F can be either deterministic or stochastic. For example, F(u, X) = u + X, with X = n denoting the Gaussian white noise, represents a noisy image. Or, when X denotes a spatial Poisson process and F is defined by

MO becomes a stochastic image with salt-and-pepper noise.

1.3. An Overview of Methodologies in Image Processing

17

(b) Treat individual images as the typical samples of some random fields, as in Geman and Geman's celebrated work [130] on Gibbs and Markov random fields. Most natural images such as the pictures of trees, cloudy skies, sandy beaches, or other natural landscapes are often more properly handled by the stochastic framework. For images with statistical nature, stochastic methods are the most ideal and relevant tools for modeling image processing. Especially important are statistical pattern theory [143, 224], learning theory [252], Bayesian inference theory for signal and parameter estimations [176], and stochastic algorithms such as Monte-Carlo simulation, simulated annealing, and the EM algorithm [33]. Combined with filtering techniques, a stochastic approach can become even more powerful in transformed feature spaces than directly in pixel domains. The Bayesian framework is worth extra highlighting among all stochastic methods due to its fundamental role. Let Q denote some observed data and F some hidden features or patterns embedded within that have generated Q. Bayesian inference of F is to maximize a posteriori probability (MAP):

The feature distribution probability p ( F ) is called the prior model since it specifies the a priori bias among the targeted patterns and is independent of data observation. The conditional probability p(Q \ F) is called the (generative) data model since it describes how the observation Q is distributed once F is specified. As far as MAP estimation is concerned, the denominator p(Q) is simply a probability normalization constant and plays no essential role. Without any prior knowledge,

alone would also give an optimal feature estimation, and is called the maximum likelihood (ML) estimator. In image processing, however, due to the huge degrees of freedom of images as high-dimensional data sets, prior knowledge becomes crucial for effective signal or feature estimation. The Bayesian estimator combines both the prior and data knowledge

Such a Bayesian principle frequently emerges in numerous image processing models in later chapters.

1.3.5

Variational Methods

Variational approach could be formally considered as the deterministic reflection of the Bayesian framework in the mirror of Gibbs' formula in statistical mechanics [82, 131]:

where ft = \/(kT) denotes the reciprocal of temperature T multiplied by the Boltzmann constant k, and Z = Z^ denotes the partition function for probability normalization. The

18

Chapter 1. Introduction

formula directly expresses the likelihood p(F) of a feature configuration F in terms of its "energy" E[F]. Therefore under any given temperature, the Bayesian MAP approach in (1.5) amounts to the minimization of the posterior energy:

where an additive constant (or ground energy level) independent of F has been dropped. When F and Q belong to certain functional spaces, such as Sobolev or bounded variation (BV) spaces, the minimization of this posterior energy naturally leads to a variational model (see, e.g., the monograph by Aubert and Kornprobst [15]). We must mention that here F denotes a general feature or pattern variable, which could further contain multiple components, say F = ( F \ , F~L, ..., FN). According to the telescoping formula of conditional probabilities, one has

In the important case when the feature components have a natural Markov-chain structure, i.e., one has

In terms of energy formulation, it amounts to

In later chapters, readers will witness this universal structure in several celebrated models including Geman and Geman's mixture image model [130], Mumford and Shah's free boundary segmentation model [226], and Rudin, Osher, and Fatemi's total variation-based image restoration model [258]. As an example, consider the following additive noise model:

Assume that (a) n is a homogeneous field of Gaussian white noise with zero mean; and (b) V« is a homogeneous random field of isotropic Gaussian white vectors with zero means. By the variational formulation, estimating F = u from Q — UQ can be achieved by

where the two weights are inversely proportional to the variances. The preceding stochastic language is a helpful metaphor or rationale for such variational formulations [73, 223, 276]. Figure 1.13 shows an example of applying the above variational model to denoise a smooth signal polluted with noise.

1.3. An Overview of Methodologies in Image Processing

19

Figure 1.13. A noisy \-D signal and its optimal restoration according to (1.7). The second theoretical foundation for the variational approach is the Tikhonov regularization technique for solving inverse problems [298]. Many image processing tasks carry the nature of ill-posed inverse problems, and the Tikhonov regularization technique is a powerful tool to make processing tasks well-posed. In essence, Tikhonov regularization is equivalent to the above Bayesian variational formulation:

where now the prior model E[F] is understood as the regularization term and the data model E[Q | F] as the fitting or fidelity term.

1.3.6

Partial Differential Equations (PDEs)

The successful applications of partial differential equations (PDEs) in image processing can be credited to two main factors. First, many variational problems or their regularized approximations can often be effectively computed from their Euler-Lagrange equations. Second, PDEs as in classical mathematical physics are powerful tools to describe, model, and simulate many dynamic as well as equilibrium phenomena, including diffusion, advection or transport, reaction, and so on [15, 117, 146, 154, 172, 254, 262, 305, 317]. By calculus of variation, for example, the denoising model (1.7) amounts to solving the following elliptic boundary value problem:

20

Chapter 1. Introduction

or dynamically via gradient descent marching, the diffusion-reaction equation

with some suitable initial condition u(x, t = 0). On the other hand, PDE modeling in image processing does not have to always follow a variational model. This is well known in physics, with famous examples of Navier-Stokes equations in fluid dynamics, Schrodinger equations in quantum mechanics, and Maxwell equations in electromagnetics. One most remarkable example of PDEs in image processing is the mean curvature motion (MCM) equation [40, 119, 241]

under suitable boundary as well as initial conditions. The MCM equation shares the double features of diffusion and advection. First, note that by naively crossing out the two | V0|'s, one arrives at the more familiar diffusion equation 0, = A0. On the other hand, define N = V0/| V0| to be the normal direction of level curves. Then the MCM equation (1.10) can be rewritten as

which is in the form of an advection equation with characteristics (i.e., particle motion) given by the curvature motion [117]

where at each pixel a, K(a) is precisely the (signed) scalar curvature of the level curve {* | (*) = 0( fl )}. More interestingly, let s denote the (Euclidean) arc length parameter of a level set and T = xs its unit tangent vector. Then K N = Ts = xss, and the above particle motion driven by curvature simply becomes

which by itself is a system of two heat equations coupled by the arc length s. Therefore, it is conceivable that under the MCM equation (1.10), all the level curves tend to become more and more regular. As a result, the MCM equation can be employed for tasks such as image denoising. Figure 1.14 displays the evolution of a trefoil shape under the mean curvature motion. Another important example of PDE processing model is Perona and Malik's anisotropic diffusion for image denoising and enhancement [251, 317]:

1.3. An Overview of Methodologies in Image Processing

21

Figure 1.14. A trefoil evolves under mean-curvature motion (1.11) or (1.12).

where the diffusivity coefficient D depends on the gradient information of u. D often inversely responds to p = V«|, for example, D(p) = 1/(1 + p2). The Perona-Malik model improves the ordinary constant-coefficient diffusion equation in its adaptivity to edge information: faster diffusion in interior regions where image changes mildly and slower diffusion in the vicinities of object boundaries. Thus the model can maintain sharpness of edges while diminishing the effect of noises, which is desirable for image analysis and processing. Well-posedness of the Perona-Malik model is, however, not so trivial as its intuitive behavior might have suggested (e.g., see the work by Catte et al. [51]). PDEs involving important geometric features are often called geometric PDEs. As one shall see in later chapters, geometric PDEs can usually be constructed by either optimizing some global geometric quantities in the variational setting (e.g., lengths, areas, and total squared curvatures) or by geometric invariance under certain transform groups. High order PDE models also frequently emerge in modern image processing [24, 26, 68, 61, 184, 262].

1.3.7

Different Approaches Are Intrinsically Interconnected

The current book attempts to present most modern image processing approaches and reveal their qualitative or quantitative connections. Some of such efforts have already been made in the literature (see, e.g., [270, 274, 283, 330]). To illustrate how different approaches are intrinsically connected, Figure 1.15 gives a very specific example on a classical image denoising model, which has just been mentioned in the preceding sections. In the stochastic setting (lower left panel in the figure), the noise model is given by uQ = u + n,

with n denoting Gaussian white noise of zero mean.

22

Chapter 1. Introduction

Figure 1.15. Different approaches are intrinsically connected: an example.

Given a typical (in the information-theoretic sense [93]) single observation w 0 , one attempts to remove noise n and extract the ideal image u. By the Bayesian MAP estimation theory, one attempts to maximize the posterior likelihood

Since the generative data model is straightforward in the Gaussian white noise case, the performance of such a MAP estimator is therefore mostly determined by the image prior model p ( u ) . Using Gibbs' ensemble models for images as random fields [130], the prior p ( u ) is then completely determined by its "energy" E[u]. In the case when (i) a standard Cartesian graph (or topological) structure is endowed to the underlying pixel lattice Q and (ii) the energy E[u] is built upon dipolar (or doubleton cliques; see Chapter 3) quadratic potentials, the MAP estimator becomes the variational model in the upper left panel of Figure 1.15 (after being normalized by some multiplicative constant):

1.3. An Overview of Methodologies in Image Processing

23

where the integrals are understood as discrete summations if the pixel domain £2 is a Cartesian lattice. In the continuum limit, it is naturally required that the compatible images u belong to the Sobolev class //' (Q,). Applying the first variation u -> u + 8u to the energy functional E[u \ UQ], one then obtains the generalized differential (or Euler-Lagrange equation)

in the distributional sense. Here the boundary term is understood as an element in the Hilbert space of L 2 (3£2, %'), i.e., all boundary functions that are square integrable with respect to the one-dimensional (1 -D) Hausdorff measure of 9^2. Therefore, by either gradient descent time evolution or directly solving the equilibrium equation, one arrives at the two equations in the upper right panel of Figure 1.15. Notice that the boundary term gives rise to the Neumann boundary conditions for both equations. Finally, from the equilibrium equation

one has for the optimally denoised estimator u,

This solution can be understood as a new way to decompose the given image UQ into two components: u and if. u belongs to the Sobolev space //' and is therefore a smooth or regular component, w, on the other hand, is oscillatory since it is the distributional Laplacian of u: w — — Aw j\ and generally belongs only to L 2 (£2). w is large along edges where u experiences noticeable jumps. In combination, locality (since differentiation is local), oscillations, and strong responses to edges qualify the w-component as a generalized wavelet projection of UQ that encodes detailed features. By analogy, the smoother w-component behaves more like a coarse-scale projection in the multiresolution setting of wavelet theory [96, 204, 290]. (The two are, however, not exactly orthogonal.) This brings us to the lower right panel of Figure 1.15. Naive as it may have sounded, this approach indeed led Burt and Adelson [35] to introduce the famous Gauss/Laplace pyramid algorithm for image coding and compression, which has been acknowledged as one of the few earliest wavelets and multiresolution ideas [96, 204, 290]. In Burt and Adelson's pioneering work [35], one should make the following adjustment in (1.13):

where Ga — G\(x/a}/o2 with G\(x} — G\(x\, jc2) denoting the 2-D canonical Gaussian. A similar observation was also made earlier by Gabor in 1960 in the context of image deblurring, as pointed out by Guichard, Moisan, and Morel [146]. In Chapter 3, readers shall also be able to see that for Besov images, there indeed exist precise variational formulations that can lead to Donoho and Johnstone's thresholding-based

24

Chapter 1. Introduction

denoising method in wavelet domains [106, 109]. Such formulations are due to the works of DeVore, Lucier, and their collaborators [55, 101]. Notable efforts have also been made by Steidl et al. [286] to understand the intrinsic connections among wavelets, PDEs, and variational models.

1.4

Organization of the Book

The current book can be divided into three major parts (see Figure 1.16): (A) Chapter 2: general mathematical, physical, and statistical background for modern image analysis and processing; (B) Chapter 3: several generic ways to model and represent images; and (C) Chapters 4,5,6, and 7: models and computation of four most common image processing tasks: denoising, deblurring, inpainting or interpolation, and segmentation. Chapter 2 introduces several general topics significant for modern image analysis and processing, which can be very helpful for understanding the materials of later chapters. It also provides independent reference sources for most contemporary works in the field. These topics are: (a) differential geometry of curves and surfaces in two and three dimensions; (b) the space of functions with bounded variations (BV); (c) elements of statistical mechanics and their implications for image analysis; (d) an introduction to the general framework of Bayesian estimation theory; (e) a compact theory of filtering and diffusion; and (f) elements of wavelet theory.

Figure 1.16. Organization of the book.

1.4. Organization of the Book

25

On one hand, these topics reveal the broad scopes of modern image analysis and processing, and on the other, they make the current book more self-contained. Readers can find numerous other monographs for these six individual topics, but few have been exclusively and systematically devoted to image analysis and processing. Chapter 3 introduces several ways to model and represent images, which in [273] the second author has characterized as the fundamental problem of image processing. Proper mathematical models for images are crucial for designing efficient models and algorithms to process them. Some image models are more general and universal than others. It is, however, often difficult to say which one offers the best model or representation, since optimality is inseparable from the concrete tasks at hand. This is quite similar to the particle-wave duality in quantum mechanics. Particle theory and wave theory are just two different approaches to modeling the same phenomenon of microscopic particles. Particle theory is more convenient in explaining behavior like the mass of photons according to Einstein's equation E = me2, while wave theory becomes superior in other occasions such as to explain the diffraction phenomenon of electrons through slim slits. To conclude, each model or representation has its own strengths and limitations. In Chapter 3, we first discuss several deterministic image models in the context of real or functional analysis, including distributions or generalized functions, L^-functions, Sobolev functions, functions with BV, as well as Besov functions for which the multiscale nature of images can be conveniently revealed using wavelet bases. We then turn to stochastic models of images, as originated from Geman and Geman's celebrated work on Gibbs' image models [130]. It is also discussed how to learn image field characteristics by combining filtering techniques with the maximum entropy principle, as first proposed in the remarkable work of Zhu and Mumford [328] and Zhu, Wu, and Mumford [329]. This chapter concludes with two other geometry-based image models: the level-set representation and Mumford and Shah's free boundary image model. These two models play fundamental roles in designing and understanding numerous variational or PDE models in modern image processing. Based on these preparations, the remaining four chapters study individually four significant image processing tasks and their models: noise and denoising (Chapter 4), blur and deblurring (Chapter 5), inpainting or image interpolation (Chapter 6), and edge extraction and image segmentation (Chapter 7). All these four tasks are in essence inverse problems. We thus usually start with explaining the corresponding "forward" problems, such as the precise meaning or models of noises, the physical mechanisms leading to various blurs and their mathematical models, as well as how to design generative or synthetic image models leading to multiple segments in general images. We then explore various approaches to their inverse problems, i.e., suppressing embedded noises (denoising), reducing blurry effects (deblurring), completing lost image information (inpainting), and partitioning into different objects (segmentation). To sum up, the current book has been organized in the logical order of

and the central line of philosophy or conclusion is

26

Chapter 1. Introduction

These characteristics well echo the Bayesian inference theory introduced in the preceding section (see (1.5)):

That is, proper prior knowledge on the target features or patterns F can introduce necessary bias among a large pool of candidates and lead to efficient schemes for estimation and decision making.

1.5

How to Read the Book

A monograph like the present book has multiple purposes: (a) to demonstrate the intrinsic logical structure behind this nontraditional novel field of image analysis and processing in contemporary applied mathematics; (b) to offer a well-organized presentation of modern image processing for interested students or scientific researchers in many related disciplines; (c) to summarize and highlight the most significant contributions by mathematicians during the past few decades, especially for our colleagues in engineering, computer sciences, astronomy, medical imaging and diagnosis, and brain and cognitive sciences; (d) to explain the vast range of mathematical as well as computational techniques and challenges in contemporary image processing, especially for general applied mathematicians working on computation, modeling, and applied analysis; (e) to reveal to the mathematicians of other areas many fascinating applications of the knowledge of topology, geometry, real, functional, and harmonic analysis, invariant theory, probability theory, partial differential equations, calculus of variations, and so on, as well as the genuine impact of mathematical power on modern societies and humanity; (f) to highlight and emphasize the interconnections among different mathematical tools and methodologies in contemporary image processing, with more weight on the variational PDE approach, an area in which the authors have made many contributions; and (g) as a whole, to serve as an effective communication channel bridging different communities related to image processing, including computer and human vision, computer graphics, signal and data analysis, statistical inference and learning, and numerous aforementioned application fields in medicine, astronomy, and scientific probing and sensoring. With such a broad range of readers in mind, the authors thus by no means expect a uniformly optimal approach to reading the present book. Each individual reader should come up with his or her own comfortable usage of the book, subject to important factors such as one's background training in imaging sciences or mathematics, one's primary focuses and motivations, as well as potential time pressure, and so on. Below we try to envision several classes of readers and offer our unnecessarily optimal personal suggestions. These categories are by no means mutually exclusive.

1.5. How to Read the Book

27

Scientific Workers with Specific Tasks in Image Processing Some readers are being driven by certain specific image processing tasks at hand. Examples include an astronomy data analyst who is bothered by blurs in telescope images, a medical radiologist who attempts to remove noises in some CT images, or some security developer who engages in designing automatic tracking software for surveillance cameras. We then recommend such a reader to jump right into the image processing chapters (i.e., Chapters 4 to 7; see Figure 1.16) that are most pertinent to the problem at hand and skip over Chapters 2 and 3 in the beginning. That is, for example, the previously imagined security developer could start directly with Chapter 7 on image segmentation (since segmentation is the most essential ingredient for object tracking). After getting the general ideas and an integrated picture of those relevant chapters on image processing, the reader can then selectively consult related topics in Chapter 2 and some necessary materials in Chapter 3. Depending on his or her mathematical background, the reader could even consult some external monographs on specific topics there, for example, differential geometry [103], elements of statistical mechanics [82, 131], or partial differential equations [117, 291]. (We have made great efforts in making the current book more self-contained, but real situations often vary among different individuals.) With these investments, the reader can then frequent those image processing chapters that are directly related to the tasks at hand, gain more detailed understanding of the modeling ideas and techniques, and eventually start to develop his or her own novel versions of models or algorithms. Graduate Student Readers As far as the current book is concerned, there are mainly two "dual" classes of graduate student readers: (A) one already with some direct experience or knowledge on image processing but less mathematical training and (B) the other with more mathematical training but lacking previous exposure to imaging sciences. For example, an average graduate student in electrical or medical engineering or computer sciences falls into the first category, while a typical mathematics graduate student belongs to the second. For graduate students of class (A) who are working on some image processing project towards theses, we recommend to first read Chapter 2 on several important mathematical elements or techniques for modern image processing. Afterwards the reader can proceed directly to the image processing chapters (i.e., Chapters 4 to 7) that are intimately connected to his or her current research project. With adequate understanding of the main lines of those chapters, he or she can then further invest time and efforts in Chapter 3 in order to gain more fundamental and systematic understanding of image analysis and modeling, as well as their application in developing good image processing models. For mathematical graduate students of class (B), we recommend to follow the existent logical structure of the book. The reader can start with Chapter 2 and skip whatever he or she is familiar with and then proceed to Chapter 3 to understand how knowledge in real, functional, harmonic, and stochastic analysis can become very useful in describing various classes of images. Afterwards comes the freedom of order in reading the last four chapters on image processing since Chapters from 4 to 7 are relatively independent. The reader is then suggested to start with whatever topic he or she is more comfortable with. For example, if

28

Chapter 1. Introduction

one is familiar with function interpolation in numerical or approximation theory, the reader can first read Chapter 6 on image inpainting. However, please keep in mind that Chapter 4 on denoising contains many important techniques described earlier and is always a good starting spot if no preference of a specific reading order arises. General Mathematicians Interested in Image Processing There are a growing number of mathematicians in pure, applied, or computational mathematics who have no previous direct involvement in image processing research but wish to broaden their own horizon of view by understanding the main ideas of image processing. We then recommend to follow the logical structure of the book (Figure 1.16). They could first selectively read topics that are most recognizably connected to their own field of interest, based upon which gradually expand interest into other topics. A geometer, for example, can first read the geometric topics (Sections 1 and 2) in Chapter 2 and perhaps first skip the details of other sections in order to directly proceed to the similar topics in Chapter 3 with geometric flavors. He or she can then peruse some selected topics in Chapters from 4 to 7 to see how geometric ingredients are built into various image processing models and what the main resultant challenges are in modeling, analysis, and computation. A computational mathematician, for another instance, may be mainly thirsty for finding challenging models to compute or novel algorithms to analyze. He or she can then quickly go over both Chapters 2 and 3 and focus more on Chapters from 4 to 7 on image processing models and algorithms. To him or her, Chapter 2 or 3 can serve as either intellectual enlightenment on new areas or good reference materials on image processing. There are indeed numerous opportunities in image processing for further computational development. Other General Scientists Interested in Image Processing The group could include, for example, psychologists working on human visual perception, computer graphics designers, or scientists who have been mainly working on acoustic signal processing and just recently become curious as to how visual signals (or images) are being processed. These scientists have already been exposed to signal or image processing concepts previously but never yet had a good chance to adequately acquaint themselves with major image processing ideas and methods. They have no specific image processing tasks at hand and are willing to invest some of their time and efforts because of the kinship between their own fields and image processing. Some mathematical tools and topics presented in Chapter 2 or 3 may in the beginning sound intimidating, but these readers should always keep in mind that tools always follow the footprints of tasks. Without good understanding on the real tasks or missions, difficulty could multiply in comprehending the meaning and languages of their tools. Therefore, such readers can start from Chapters 2 and 3 to get basic tones of the book but should never be discouraged by some of the mathematical details, notations, or formulations that are temporarily beyond their comprehension. They can thus in the beginning spend more time focusing on Chapters 4 to 7 on four very specific but important image processing tasks. The tasks themselves are often relatively much easier to grasp.

1.5. How to Read the Book

29

When coming to the approaches to carrying out those tasks, the readers then become highly motivated in terms of the main challenges and key ideas. It is from this point on that our readers can start to use Chapters 2 and 3 and gradually understand in depth the real motivations, advantages or disadvantages, and the key ingredients of the mathematical entities. Depending on their mathematical background, these readers are also encouraged to read some external introductory materials on specific mathematical topics in Chapters 2 and 3. Finally, the authors wish that our readers can genuinely enjoy the numerous unique views and state-of-the-art methods presented in the book. Together we can further advance this exciting novel field which is benefiting profoundly not only contemporary sciences and technologies but also the society, civilization, and the very nature of humanity. Seeing is not passive watching or reception but active believing, inferencing, and decision making based upon conscious or subconscious computations in the human brain.

This page intentionally left blank

Chapter 2

Some Modern Image Analysis Tools

In this chapter, we help readers review several major tools in modern mathematical image and vision analysis, whose fundamental roles will manifest in later chapters. These topics reveal the geometric, stochastic, real analysis, and harmonic analysis aspects of this emerging field. In addition to providing useful background mathematical knowledge for later chapters, the present chapter can also serve as a valuable reference source for researchers in the field.

2.1

Geometry of Curves and Surfaces

Curves and surfaces are basic geometric elements in image and vision analysis, and computer graphics. For example, they define and reveal the information of an automobile (in automatic traffic control), a star or planet (in astronomic imaging), a human body (in video surveillance), and an internal organ (in medical imaging and 3-D reconstruction). In this section, we review some basic theories on curves and surfaces in two or three dimensions. Readers are also referred to the classical geometry book by do Carmo [103] or some other monographs on geometric image analysis and computation, e.g., Romeny [256] and Kimmel [173]. 2.1.1

Geometry of Curves

Local Geometry A parametric curve is often denoted by x(t), with t varying on an interval andx = (x\,X2) or x — (XI,XI,XT,) depending on the dimension. Differentiation with respect to t will be denoted by an overhead dot. We shall assume that x(t) is smooth enough, and that the parametrization is regular, meaning that x ( t ) is nonzero for all t. The local infinitesimal displacement is given by dx = xdt, leading to the Euclidean parameterization ds = \dx\ = \x\dt. Differentiation with respect to the arc length df/ds will be denoted by /' in this section for any function / defined on the curve. Then t = t(s) = x'(s) is the unit tangent vector. In two dimensions, rotation of the tangent by 90 degrees (see Figure 2.1) gives the normal vector n. More generally in 31

32

Chapter2. Some Modern Image Analysis Tools

Figure 2.1. A planar curve: tangent t, normal n, and curvature K. differential geometry, the normal is a second order feature, as manifest in the curvature formula In three dimensions, the displacement of n further introduces the third degree of local geometry, i.e., the binormal vector b and the torsion scalar T. Concisely, these combine into the well-known frame formula of Frenet [103]:

The skew-symmetry of the connection matrix results directly from the orthonormality of the moving frame (t,n,b). The torsion and binormal often appear in computer graphics and designs, while in 2-D image processing, curvature plays a more explicit role. In vision psychology, it has been observed that human subjects can easily observe and process convex or concave features of shapes and boundaries [163, 170]. That is to say, our biological vision networks are indeed able to evaluate and process the curvature information. The importance of curvature in image and vision analysis is perhaps rooted in the most fundamental laws of physics. From Newton's law in classical mechanics, Hamiltonian dynamics, and the Second Law in thermodynamics to Schrodinger equations in quantum mechanics, physics seems to be perfectly complete working with up to second order information. Take the Second Law for example; first order derivatives (of free energies) often define some key quantities such as temperature, pressure, and chemical potential, while the second order ones well explain their statistical fluctuations and introduce all the rest of the key quantities such as heat capacity and thermal compressibility. Though direct connection has not yet been established, we do believe that curvature-based vision cognition must result from the physics or biophysics of the vision neural networks (see, e.g., the work by Hubel andWiesel[152]). Under the arc length parametrization, assuming K(S) > 0 and applying the exterior product, we have

2.1. Geometry of Curves and Surfaces

33

Back to the general parametrization x(t), it gives

As an application, consider the graph curve of a function y = f ( x ) on the x-y plane. Take t = x as the natural parametrization. Then

Therefore, the curvature is given by

with positivity corresponding to convexity. For some applications in computer graphics or medical imaging (e.g., the action potential and its shock front on a normal heart), one is interested in the behavior of curves on curved 2-D surfaces instead of planes. One could generally study parametric curves x (t) on a smooth Riemannian manifold M [85, 103]. Then the tangent x belongs to the tangent space TXM, which is endowed with some local Riemannian inner product {• , •}. The arc length parameter is derived from

On a general Riemannaian manifold, high order derivatives are taken based on its LeviCivita connection D. For a local vector field X defined on a neighborhood of x e M, D X is a type (1,1) tensor in TXM, meaning that for each v e TXM, DVX e TXM. Thus the LeviCivita connection generalizes the ordinary (Euclidean) notion of directional derivative: DVX is the rate of change of X along the direction v. In particular, one could define the Levi-Civita derivative x"LC(s) = Dx>(S)x'(s). Like straight lines in Euclidean spaces, a curve with x"LC(s) ~ 0 for all s is interpreted as always heading straight forward on M, or simply, a geodesic. In image and vision analysis, most 2-D surfaces are naturally embedded in R3, and their Levi-Civita connections become more explicit. Let M denote a 2-D surface embedded in R3, and let x (t) be a parametric curve on M. Then x (t) has a natural representation under the canonical coordinates of R3: x = (x, y, z). Taking the ordinary first derivative in R3 gives the tangent x(t) & T X M . Under the arc length parameter s, the second order derivative in R3: x"(s) = (x"(s), y"(s), z " ( s ) ) often "sticks out" from the tangent plane TXM. To obtain the intrinsic or Levi-Civita second order derivative x"LC, one applies an extra step of orthogonal projection: where N is the normal of TXM in R3. In particular, a "straight line" with x"LC = 0 on M can still bear a hidden degree of bending, namely, along the direction of N. The big circles on a sphere are household examples. Furthermore, suppose that the surface M is defined as the zero level set of some function 0:

34

Chapter 2. Some Modern Image Analysis Tools

Let x(s) be a parametric curve on M with arc length s. Then (x' , V0) = 0, where V denotes the gradient. Taking differentiation again leads to

where D2 denotes the Hessian operator and the second term for the Hessian bilinear form. Notice that the surface normal is given by N = V0/| V0|. Therefore,

where the last step follows from (2.2). One clearly sees how the geometry of the surface and that of the curve interact. Variational Geometry In many situations, a collection of curves have to be dealt with simultaneously. In image segmentation for example, one attempts to select the "best" candidate among all possible boundary curves [75, 226, 299]; also for image interpolation or disocclusion in computer vision [67, 115, 116, 234], missing or incomplete boundary curves are usually fitted by some optimal candidates. In such situations, it is more ideal to treat a curve as a "point," out of an ensemble of candidate curves. In order to weigh the quality of each individual curve y, the first step is to establish some suitable measure E [ y ] . For curves generated by random walks, E can be specified by probability measures. For Brownian paths for example, the Wiener measure becomes natural and even unique under suitable assumptions. In this section, we shall deal with smooth curves, and regularity functionals are then the ideal candidates for E [ y ] . Let x(t) be a parametric curve with t e [a, b\. Then the first natural geometric measure is the length L:

At first sight, it seems necessary that x be C', or at least with L' integrable x. But the length L can be interpreted as the total variation of x(t), and it therefore suffices that x(t) belongs to the space of bounded variation (B V). Under variational processes, global geometry always points to local geometry. For length, consider a small perturbation of a given curve: x(t) —>• x(t) + Sx(t), a < t < b. By Taylor expansion, up to the linear order of Sx, the change of length 8L is given by

2.1. Geometry of Curves and Surfaces

35

Integration by parts further leads to

Thus away from the two ends, the rate of change of the tangent, —i, characterizes the sensitivity of length growth. Under the arc length parameter s,—t' = —Kn exactly defines the curvature information. A general second order (Euclidean) measure for a smooth curve y is given by

where /(/c) is some suitable function of curvature. In two dimensions for instance, /(/c) = K leads to

where 0 is the signed angle of t with respect to a fixed reference direction (such as the .x-axis) and changes continuously. For a closed smooth curve, E\ has to be an integer multiple of 2n, and is thus called the total circulation. In many situations, E\ is not a good measure since it does not penalize local oscillations as long as they cancel each other. If /(/c) = |K |, then the total absolute curvature

is exactly the total variation of the heading direction 9, It is not as memoryless as the circulation measure E\, in that local oscillations do not cancel each other. The total curvature measure £2, though capable of memorizing local ripples, still does not effectively discriminate local behavior and global trends. For instance, E2\y] equals 2n for any simple and closed convex curve y, including one that contains kinks or corners. In order to more effectively penalize local singularities, one weighs the arc length element ds by some w(|/c|), with w being nonnegative and increasing with respect to |/c|:

In the case when W(\K\) = \K\P~] for some p > 1, one obtains

A special but very useful case is Euler's elastica measure Ee\:

36

Chapter 2. Some Modern Image Analysis Tools

for some positive constant weights a and b. It first appeared in Euler's work in 1744 on modeling torsion-free elastic rods and has been reintroduced into computer vision by Mumford [222] for measuring the quality of interpolating curves in disocclusion. The ratio a/b indicates the relative importance of the total length versus total squared curvature. Variational calculus leads to the equilibrium equation of the energy Ee\:

which is nonlinear and could be solved explicitly using elliptic functions, as Mumford developed in [222].

2.1.2

Geometry of Surfaces in Three Dimensions

Local Geometry In abstract differential geometry, surfaces are 2-D Riemannian manifolds on which inner products are defined pointwise and smoothly. Their local differential structures are identical to R2. Klein bottle is a well-known example which, however, cannot be faithfully visualized infl 3 [103]. In computer graphics and 3-D image reconstructions, surfaces are always naturally embedded in their 3-D environment /?3, which makes their theory and computation relatively simpler. There are two alternative ways to define a surface embedded in /?3, parametric or implicit. A parametric surface is defined as a smooth map x(u, v) from a parameter domain ft c R2 to /?3: with the nondegeneracy condition xu x xv ^ 0. As a familiar example in image analysis, a color image could be considered as a parametric surface in the RGB color space, with the two Cartesian coordinates in the image plane as natural parameters. An implicit surface, on the other hand, is defined as a zero level set: (#) = 0 for some suitable function 0, often with the nondegeneracy condition V0 ^ 0 on the surface. By the implicit function theorem, an implicit surface could be locally parameterized. On the other hand, a parametric surface could also be made implicit locally using the distance function. Let £ c R3 denote the image of a parametric surface x(u, v). Define the distance function d^ by Then £ is exactly the zero level set of d-z. The distance function has many applications in image and vision analysis, especially in the level-set method of Osher and Sethian [241]. Consider a parametric surface x(u, v) : £1 —» R 3 . The tangent plane Tx is spanned by xu andxv. Define

Then the inner product on Tx is completely determined by the 2-by-2 structure matrix M — (E, F\ F, G), which contains all the first order information of the surface. A parametric

2.1. Geometry of Curves and Surfaces

37

curve on the surface is in the form of x(u(t), v ( t ) ) . Computing its arc length leads to the so-called first fundamental form /(-, •) of the surface:

Generally, I(u, v) is precisely the inner product for all tangent vectors u, v e TX*L. To study the second order geometric information of a surface E, one has to understand how to take derivatives for tangent vector fields. In general Riemannian geometry as already mentioned, this is realized by the Levi-Civita connection, or the covariant differential operator D. For a given tangent vector field X on Z and a direction v e TxE, DVX e 7^ E generalizes the ordinary notion of directional derivative of X along v. Due to the product rule of the covariant derivative,

and the linearity of Dv on v, it suffices to understand the following four basic actions:

If the surface is flat, these are simply the ordinary second order derivatives xuu, xuv, and X vv

For a general embedded surface, xuu may stick out from it. But creatures living on the surface cannot sense what is away. Therefore, xuu needs to be projected back onto the surface: where N denotes the unit normal in the direction of xu x xv. (If the codimension of the embedding is not 1, N has to be replaced by the entire normal space.) This discussion certainly applies to the other three as well. The quantities (x.. , N) completely measure the deviation of the covariant derivatives from the ordinary ones. Therefore they encode the local information of how far the surface deviates from being a plane, i.e., the bending and curving of the surface. Conventionally, they are denoted by

Since (x. , N) = 0, one has

In combination, they define the second order geometric information and lead to the second fundamental form //, which is defined as follows. By abusing the covariant derivative symbol D, for any v e Tx E, let us denote by Dv N the ordinary directional derivative of N along v in R3. Since along v (i.e., along a curve segment on E which is tangent to v at x), (N , N) = 1, one must have (DVN , N) = 0, implying that DVN e T^E! Thus D.N could be considered as a linear transformation in TXY. Indeed, —D.N is called the Weingarten linear map [103]. The second fundamental (bilinear) form in Tx E is then defined as

38

Chapter 2. Some Modern Image Analysis Tools

Figure 2.2. Second order local geometry: surface normal N and two perpendicular principle curvature circles. Noticing that DXu N = Nu and DXiN — Nv, one has

In particular, dx = xudu + xvdv and

The second form is thus symmetric as well. In order to understand how the second fundamental form encodes second order geometric information, consider a unit vector u e Tx S. The two orthonormal vectors u and N span a unique plane though jc, which locally intersects Z along a planar curve y. Applying Frenet's frame formula (2.1) to y at x, one has

In the current context, one has

Therefore, II(u, u) is exactly the curvature of y, which is denoted by KU and gives the explicit geometric meaning of the second fundamental form. Furthermore, one can introduce the Rayleigh quotient for the symmetric second form as in linear algebra [138, 289]. For any u e TXZ [103], define

Generically, its maximum K\ and minimum KI are called the principal curvatures (see Figure 2.2). Gaussian curvature is defined as K = K\K2, while mean curvature H = (K\ + K2)/2. From the perception point of view, it is unclear to which degree of precision the human vision system can detect these two curvatures. Psychological evidences show that human vision is at least good at detecting the sign of K, which corresponds to the visual differentiability between ellipticity (i.e., convexity) and hyperbolicity (i.e., saddle structures).

2.1. Geometry of Curves and Surfaces

39

Let M and L denote the two structure matrices of the first and second forms under the natural parametric basis (xu, xv):

Then from the computational point of view, the two principal curvatures are exactly the two eigenvalues of the generalized eigenvalue problem (L, M): Lu = KMu. Therefore, by linear algebra [139, 289],

or more explicitly,

The ubiquitous denominator EG — F2 = det M is the squared area of the parallelogram spanned by xu and xv. In particular, if xu and xv just happen to be orthonormal at x, then

The graph of a gray-scale image h(u, v) : Q -> R could be considered as a parametric surface in (u, v, h) e R3. Then it is easy to establish that

which is the determinant of the Hessian of h, normalized by the graph area element yi + hi + h\. Similarly for the mean curvature,

Figure 2.3 shows a surface patch and its associated Gauss and mean curvatures. Now consider an implicit surface £ defined by a potential function 0:

Working near a fixed target point :c, one may take any parametrization x (u, v) locally near x on E. In addition, assume that at this target point (only!), (xu, xv) is an orthonormal system of Tx E. Differentiating (f)(x(u, v)) — 0 gives the first order information:

Taking differentiation again to the first equation with respect to u gives

40

Chapter 2. Some Modern Image Analysis Tools

Figure 2.3. A surface patch z = h(u, v) = cos(w) cos(t') (top), its mean curvature field H (lower left), and Gauss curvature field K (lower right). where D20 denotes the Hessian matrix in /?3 and its associated bilinear form. Noticing that the surface normal /V = V(/>/| V|, we obtain the formula for e in the second form:

Similarly, one has

Since (xu, xv) is an orthonormal frame at the target point x, (2.5) leads to

These two beautiful formulae are independent of the parametrization and completely selfcontained. The mean curvature H has been frequently studied in dynamical processing of shapes and surfaces, mainly due to its variational meaning (see the coming subsection). In the literature, another more explicit formula for H is also commonly used based on the trace formula The first term is exactly the 3-D Laplacian

Therefore, —2H equals

and the mean curvature is thus the divergence of the normal field N (up to a constant)!

2.1. Geometry of Curves and Surfaces

41

Variational Geometry One common global feature of a given surface £ is its area A. Under the parametric form x(u, v) : £2 —»• H, the area element is given by da = \xu x xv\dudv. Thus

As for curve length, under a small perturbation,

the first order variation of A must be in the form of

where H is a vector field on E that takes values in R3. Its exact form can be worked out as follows. Let X = xu x xv. Then the normal is N = X/\X\. Notice that 8\X\ = (N , 8X) and

From the vector calculus of inner and exterior products, one has

Applying integration by parts and noticing that 8x(u, v) has been assumed compactly supported, one has

Compared with (2.10), it gives

Notice that all the four vectors on the right are in the tangent plane Tx E. Therefore the two exterior products are both in the normal direction N. Assume

for some scalar H. Then by (2.11),

42

Chapter 2. Some Modern Image Analysis Tools

Notice that the denominator is exactly EG — F2 from the first fundamental form. Applying the two rules in vector calculus

one establishes that

where e, f , and g are the coefficients of the second form. Referring to the mean curvature formula (2.4), we end up with Therefore, we have shown that

which is the global (or variational) meaning for the local mean curvature H. It also indicates that the most area-efficient displacements are along the normals, which well coincides with the daily life intuition. Consider, for example, the (top) surface area of a glass of still water. Smooth and slow rotation (i.e., 8x's are on the tangent planes) brings no noticeable change of surface area, but any wave ripples in the vertical direction do. A surface with mean curvature H = 0 is called a minimal surface, which has been a fascinating subject in geometry, geometric measure theory, and nonlinear PDEs ever since J. L. Lagrange [233]. Minimal surfaces are also very useful in computer graphics and geometric designs [92]. In image analysis, one often considers the graph surface Z/, of a given gray-scale image h(u,v) on an image domain £2. Then the area is given by

where the convenient notation |X| fl denotes ^/a2 + \X\2 for any vector X . In this situation, we are interested in the sensitivity of A(h) with respect to small changes of h: h —> h+8h. Noticing that 8\X\\ = ( X / \ X \ \ , 8X}, we have for the first order variation of A(h)

where the perturbation 8h has been assumed to be compactly supported. It was first derived by Lagrange in 1762 (see more detailed discussion in [233] for example). Notice that for image graphs, x = (u,v,h(u, v}}, 8x = (0, 0, 8h), and (N , 8x} = 8h/\Vh\\. Comparing (2.14) with the most general form (2.13), one concludes that the surface mean curvature is given by

2.1. Geometry of Curves and Surfaces

43

For smooth surfaces, we could consider high order global features as well. For example, the total curvature

the total absolute curvature

and the total squared curvature

More generally, let f ( K , H) = g(K\, ^2) be a proper function of the two principal curvatures. Then one could attempt to minimize the cost function

Many interesting results exist regarding these energies. For example, suppose the Gaussian curvature K does not change signs. Then from the Gauss map theorem [103]

one has

which is exactly the total parametric area of the normal map N : £2 -> S2. Suppose N is injective. Then this latter area could never be bigger than 4jr, the total area of S2, which implies that EI < 4rt in this situation. Sometimes one may also become interested in energies that are not intrinsically geometric, meaning that they do depend upon specific representations. For example, for a graph surface representation x(u, v) = (u, v, h(u, v)), one could attempt to minimize

By the total curvature formula (2.6),

which is not intrinsic due to the h-factor but is indeed a close approximation to the total squared curvature £3 when Vh is small. Finally, we would like to mention that besides curves and surfaces, high-dimensional differential or Riemannian geometry could also be very useful in image processing (see, e.g., [174]).

44

2.1.3

Chapter 2. Some Modern Image Analysis Tools

Hausdorff Measures and Dimensions

In image and vision analysis, one has to deal with curves and surfaces that are less regular, for example, fractal landscapes. Then differential geometry often demands too much, while the notions of Hausdorff measures and dimensions can become more powerful. Let E be a given set in R2 or /?3, or more generally any subset in a metric space with distance function. Suppose it is known a priori that £ is a d-dimensional object, with d = 1.618 say. Then a natural question is how to measure its J-dimensional "volume." Hausdorff measure Jid provides one good answer. A covering A of E is a (countable) collection of subsets A's whose union includes E as a subset. The scale ||^4|| of A is defined as

where diam(A) stands for the diameter of A. Furthermore, one defines

Then the J-dimensional Hausdorff measure 1-Ld(E) is defined as

where A's must be coverings of E. Notice that the e limit is well defined since the infimum as a function of € is nonincreasing. That is to say, the € limit process is in fact a supremum one. Thus the definition is a minimax statement. It is easy to see that all the ingredients in the definition are natural. (a) When d is an integer and A is an ordinary ^/-dimensional cube or ball, diam( A)d correctly measures the volume up to a multiplicative constant. (b) Taking infimum is to diminish the effect caused by inefficient overcovering or repetition. (c) Finally, by demanding 6 to vanish, a Hausdorff measure can potentially resolve details at all scales (see Figure 2.4). From the measure-theoretic point of view, T-Ld defined above is only an "outer" measure. Any outer measure, once restricted on the class of measurable sets, becomes a true measure. It is easy to see that for any d,t > 0,

from which one has, for any 6 > 0, Therefore, as long as T-Ld(E) is finite for some d, H.S(E} = 0 for all s > d. Therefore, there exists at most one d, so that l-Ld(E) is nonzero and finite, which leads to the notion of Hausdorff dimension dim//:

If dim//(E) > 0, then for any t e [0, dim//(E)), n'(E) = +00 (see Figure 2.5).

2.2. Functions with Bounded Variations

45

Figure 2.4. A set E (left) and a covering A (right). Covering elements like the huge box with dashed borders (right) fail to faithfully capture small-scale details. The definition of Hausdorjf measures forces covering scales to tend to zero.

Figure 2.5. The Hausdorff dimension dim//(£) of a set E is the critical d, above which E appears too "thin," while below it appears too "fat. "

2.2

Functions with Bounded Variations

Functions with bounded variations (B Vs) are ideal deterministic image models which allow the existence of jumps or edges, and are, however, still mathematically tractable. Theory of BV functions and its extensions [7, 9, 118, 195, 137] provide the necessary mathematical foundations for celebrated image processing models such as Rudin, Osher, and Fatemi's total variation denoising scheme [257, 258] and Mumford and Shah's segmentation model [226].

46

2.2.1

Chapter 2. Some Modern Image Analysis Tools

Total Variation as a Radon Measure

Let £2 c R2 denote a bounded open domain and u = u(x, y) belong to L 1 (£2). Motivated by both analysis and image processing, £2 is often assumed to be a Lipschitz domain. If u is smooth, then its total variation (TV) TV[w] is defined as

A function (or an image) u with TV[w] < oo is said to have B V. The notation B V(ft) denotes all functions in L 1 (ft) with BV. How smooth should a function be for the definition (2.18) to make sense? At first sight, it appears that at least u should belong to the Sobolev space WIA (ft), i.e., functions with integrable first order partial derivatives. But the power of TV and BV in image analysis exactly arises from the relaxation of such constraints. That is, the space BV(ft) is in fact much larger than W 1 - 1 (ft) to include many interesting image functions. To get acquainted with such relaxation, let us first look at a related but simpler example. For any function f(x, y) e L 1 (ft), we could define the total integral / by

The condition L 1 seems natural but not exclusively necessary for /[/] to be well defined. For example, assume that (0, 0) € ft, and denote by 8(x, y) the Dirac delta function. Then one can still define More generally, let \JL be a nonnegative measure on all Borel sets of ft with p.(K) < oo for all compact subsets K c ft. Such measures are examples of Radon measure [126, 118]. Then one could still define, as for Dirac's delta,

where /z could be generally a signed Radon measure. Signed Radon measures generalize the space of L1'oc(ft). Furthermore, in a certain proper sense (i.e., linear functionals) [126], it could be shown that signed Radon measures are indeed the most general class of mathematical entities for /[•] to be well defined. The same story applies to the TV functional. In some sense, TV[w] is the integrated version of /[/] by taking / = VM. Consider the case in Rl when /(*) = 8(x). Then u' = 8 would imply that u = H(x)—the Heaviside 0-1 function, provided that u(—00) = 0. Since H ( x ) contains a jump at the origin, it stays outside the locally Sobolev class H^,1. Inspired by the discussion on /, one could, however, still define

Generally, let ^t(jc) be any finite and signed Radon measure in /?', and let u = u(x) denote its cumulative distribution function (c.d.f. in probability theory). One could then define

2.2. Functions with Bounded Variations

47

Figure 2.6. Three \-D images with TV[/] = TV[g] = TV[/z] = 2. Here |/z| denotes the TV measure of /JL. That is, if /z = /JL+ — ^T denotes the Jordan decomposition into two mutually singular and nonnegative Radon measures [126], then \fji\ = ^L+ + IJL~. Figure 2.6 shows the TVs of three 1-D signals. With these insights in mind, we now formalize the definitions of TV[w] and BV(£2). Let U c Q c R2 denote any open set of £2. For any u e Ll (£2), define

where B2 denotes the open unit disk in /?2, and the admissible class C\(U, B2) denotes

If / | D*u | < oo for any open subset U with compact closure contained within £2, we say that locally u has BV and denote by BVioc(£2) the collection of all such functions in Ll(Q). For any u e BVioc(£2), by the standard practice in measure theory, one could construct an outer measure by defining, for any subset E c Q,

where U are open subsets. Restricted on measurable sets, it becomes a Radon measure on Q, which shall be denoted by / \Du\, and defines the TV of u (on £2) by

The notation BV(£2) denotes the collections of all functions with finite TVs in L 1 (fi). Assume that u e W M (n). Then for any g e C'(^, fi2),

48

Chapter 2. Some Modern Image Analysis Tools

by the Gauss-Green divergence theorem (or integration by parts). Since B2 is closed under radial reflection g -» — g, one has

On the other hand, by Lusin's continuity characterization, Uryson's extension lemma, and the mollification technique [126], for any specified precision e and -) = AH(x — 1/2), a product of an amplitude constant A > 0 and the shifted Heaviside function. It represents a binary black-white image of two adjacent vertical bands. Then

where D[ means 1-D TV. It is not a coincidence in this case that the TV measure is the product of the size d — c and strength A of the target edge. Example. Let £2 = /?2, and let aB2 denote the disk centered at the origin with radius a > 0. Let XaB2 = XaB2(x, y) be its indicator function. Then by definition,

where 3 (A) stands for the (topological) boundary of a set A, K 1 the 1-D Hausdorff measure, and n the outer normal. On the other hand, it is easy to construct a g e C\ (R2, B2), so that

2.2. Functions with Bounded Variations

49

Figure 2.7. An example of L 1 -lower semicontinuity. The sequence (un) of 1D images on [0, 1] converges to u = 0 in L 1 smce ||«n+i — u\\i\ < 2~n. Notice that TV(«) = OwhileTV&n) = 2, which is consistent with the property of lower semicontinuity: TV(w) < liming TV(« n ). In particular, strict inequality is indeed realizable. This result holds for any smooth domain E; namely, the TV of its indicator function XE is exactly its perimeter. In fact, this approach leads to a more general definition of perimeters for nonsmooth domains [137, 195].

2.2.2

Basic Properties of BV Functions

In this section, we discuss several major properties of BV functions that are important in image analysis. They are L 1 -lower semicontinuity, completeness (i.e., Banach), trace and divergence theorem, smoothness and regularity, and relative compactness in L'. Theorem 2.1 ( L1-Lower Semicontinuity). Suppose un (x, v) -> u(x, y) in L' (£2). Then

In particular, if(un)n is a bounded sequence in BV(£2), then u belongs to BV(fi) as well. The proof almost immediately follows from the definition of TV. For any g e Clc(Q,B2),

Taking supremum confirms Theorem 2.1 (see Figure 2.7). Among the many interesting results that Theorem 2.1 makes contribution to, perhaps the most fundamental one is that BV(£2) is a Banach space. First, it is easy to see that BV(£2) is a normed linear space under the BV norm:

which is stronger than the L 1 norm. To establish the completeness, let (un)n be a Cauchy sequence under || • ||BV- Then (un)n is also Cauchy in L 1 , and let M e L 1 be its limit. By

50

Chapter 2. Some Modern Image Analysis Tools

semicontinuity, u e BV(ft). Furthermore, applying semicontinuity again to (u,n — un}n with ra fixed, we have

which converges to 0 as m —> oo since (un)n is a Cauchy sequence. Thus un —>• u in BV(ft) under || • ||BV. The behavior of a BV function u 6 BV(ft) near the boundary 9ft is reflected in the notion of trace fu = T(u) e L ' ( 9 f t , Hl). In what follows, we shall assume that ft is a Lipschitz domain. Roughly speaking, the definability of /„ = T (u) indicates that a B V function u cannot oscillate too much near the boundary 9ft. Instead, it converges to some function (in a certain suitable sense) / along 9£2. To see clearly why this is the case, let K\ c K^ c • • • be a sequence of compact sets so that ft = [_J^° Kn. Then by the basic continuity property of a finite measure,

Therefore, roughness indeed gradually fades away on the "narrow rings" £2 \ Kn, and the function tends to converge. More explicitly, as far as a local segment of rectifiable boundary is concerned, one could "zoom in" and simply assume that ft = R x (0, oo). Take K€ = R x (e, oo), e > 0, and define f £ ( x ) = u(x, s) in the L 1 sense (by Fubini's theorem). For any 8 < e,

Therefore ( f s ) £ is a Cauchy sequence in L' (/?), and the trace fu = T(u) along the boundary is defined as its L 1 limit. The concept of trace leads to the Gauss-Green divergence theorem for BV images. Suppose that u € BV(£2), U d £2 (i.e., U is compact in ft) and is Lipschitz, and g e C,S(ft,IR 2 ). Then

where n denotes the outer normal of dU. As a Radon measure, for any compact set K c ft,

In image analysis and processing, of particular interest is the case when K is a compact piece of Lipschitz curve y. y is said to be singular for a given image u if

Notice that under the 2-D Lebesgue measure, y is always only a null set.

2.2. Functions with Bounded Variations

51

Let / + and /~ denote the corresponding traces of u along the two sides of y. Then one can show that the integration under the Hausdorff measure of the jumps along the curve! In image and vision analysis, intensity "jumps," usually referred to as "edges," have been considered as crucial visual cues. Equation (2.21) therefore partially explains why the TV measure and the B V space could have become so convenient and powerful in image analysis and processing. Another key question is how close W 1 ' 1 ^) is to BV(fi). By adaptive mollification based on the partition of unity [126], one could establish the following result. Theorem 2.2 (Mollification of BV). For any u e BV(£2), one can find a sequence of approximations (un)n such that 1. un e C°°(Q)forn = 1 , 2 , . . . ; 2. un -H» u in L' (Q) as n —> oo; 3. JR Dun | —> fft | Du \, and their traces gn = g for n — 1 , 2 , . . . . Theorem 2.2 is the best that one could wish for, since for a general BV image u it is infeasible to ask for a sequence of smooth functions (un)n so that

This could be easily explained by (2.21). Suppose u e BV(£2) and that a compact regular curve y c £2 is singular, i.e., f \Du\ = a > 0. Then for any smooth function «„,

That is, convergence under the TV seminorm is impossible for such an image. However, mollification provided by Theorem 2.2 has already been quite powerful. It allows one to conveniently transplant many results from the classical Sobolev space Wl'l(tt). Let u € BV(£2) and / = T(u) be its trace. A functional L[u, f] on BV(£2) is said to be L 1 -lower semicontinuous if, for any sequence (un}n e BV(£2) with un —> u in L 1 (£2) and /„ = / as n —> oo, one has

Corollary 2.3 (SoboIev-BV Transition). Suppose L[u, f] is an L[-lower semicontinuous functional in BV(£2) and that E = E(t) is a continuous function oft e [0, oo). If then the same inequality must hold for all u e BV(£2).

52

Chapter 2. Some Modern Image Analysis Tools

The proof follows readily from the semicontinuity property and Theorem 2.2 on mollification approximation. As an application, define in R" (with n > 1) p = ^-j-, and

Then the celebrated Sobolev inequality in W I A ( R " ) states that L < E for all compactly supported u e WIA(R") for some suitable constant c which depends only on n. By the preceding corollary, the same inequality must hold for any compactly supported u €. BV(£2) as well. The second most important application of the mollification theorem (Theorem 2.2) is the weak compactness property, which has become a fundamental tool for establishing the existence of solutions to various BV-based image processing models. Theorem 2.4 (Weak Compactness of BV). Let (un)n be a bounded sequence in BV(£2) where Q is a Lipschitz domain. There must exist a subsequence which converges in L 1 (Q). The proof is again immediate from the mollification theorem (Theorem 2.2) and the weak compactness of W L 1 (£2) in L'(£2). By Theorem 2.2, we could approximate un by wn € W{A(Q) such that

for each n. Thus (wn)n must be bounded in WL[(Q). By the weak compactness property of W L { ( £ 2 ) (or the Rellich theorem [3, 137, 193]), (u;,,),, contains a subsequence indexed by n(k), k = 1, 2 , . . . , which converges in L 1 (£2). Then the subsequence of (u,,)n with the same indices n(k) must also converge in L 1 (fi).

2.2.3

The Co-Area Formula

This beautiful formula was initially due to Fleming and Rishel [125], and De Giorgi gave the complete proof [134, 137]. It reveals the geometric essence of the TV measure, which has been fundamental for numerous applications in image analysis and processing. For smooth functions, it essentially results from a special change of variables, at least locally. Let u(x, >') denote a smooth image on a domain Q. Assume that Vw(jc 0 , >'o) / 0. Say uy(xQ, yo) 7^ 0, without loss of generality. Then by the implicit function theorem, locally near (JCQ, >'o), we could explicitly solve v:

For each given X, one could reparameterize the 1-level curve (x, y(x, 1)) using arc length 5, so that Since the reference points with 5 = 0 may also depend on A, one has

2.2. Functions with Bounded Variations

53

which is well defined at least locally. In summary, this change of variables is governed by

In the (s, A) system, taking partial derivatives with respect to s and A in the first equation gives The first identity, combined with the second one in (2.22), gives Vu = ±| Vu\(—ys,xs). Therefore, the second identity becomes

and with the Jacobian, we are able to conclude that

Assuming that this change of variables is global on £2, one has

where the second integral is over an appropriate domain in the (s, A) plane. Let y^ denote the A-level curve. Then

which is the celebrated co-area formula for regular functions. More generally, for any function

\ , 0)2) decays sufficiently fast as \co\ = J u>\ + w\ —>• oo.

66

Chapter 2. Some Modern Image Analysis Tools

Here K stands for the 2-D Fourier transform

In signal processing, they are called the lowpass conditions, since in combination when applied to images, such PSFs almost preserve low frequency components (condition (a)) while suppressing high frequency ones (condition (b)). A radially symmetric PSF K is said to be isotropic. That is, there exists a single variable function k(r2) so that K(x, y) = k(x2 + y 2 )- More generally, K is said to be orientation selective or polarized, if

for some positive definite matrix A = [a, b; b, c]. A nonnegative PSF ^(jc, y) > 0 for (x, y) e R2 could be naturally treated as a probability density function due to condition (a), which is also equivalent to the Markov condition for Markov transitions, as to be made clear below. Given a PSF K(x, y } , the associated spreading or blurring transform for images is given by the convolution

In the Fourier domain, it is simply v = K • u. In particular, u(0, 0) = w(0, 0), which is exactly the conservation law if u (x, y) is understood as a photon density function for optical signals:

The spreading behavior could be easily visualized in the ideal case when u = 1 K (x, y) and K = l r (jc, y)/(nr2), the indicator functions of disks BR and Br centered at (0, 0) with radii R and r. Then it can be easily shown that v = K * u is supported over BR+r and vanishes outside. That is, the original image information on BR spreads out to an expanded region BR+r. When K is nonnegative, point spreading carries the beautiful stochastic interpretation by Markov transitions or random walks. Assume that u(x, y) denotes the spatial density of photons. Imagine (which may not be physically plausible!) that each photon is moving randomly in the (jc, y)-visual field and is subject to the one-step transition law

which describes the percentage of photons in the infinitesimal box (p, p + dp) x(q,q+dq) to be observed within (x, x + dx) x (y, y + dy) after one step. The Markov condition for random walks is then automatically satisfied:

2.5. Linear and Nonlinear Filtering and Diffusion

67

which simply says that no photons are created or annihilated, and each photon either stays at the original place or transits to a new location after one step. Suppose the initial photon distribution is given by u(x, y), and denote the new distribution after one step by v(x, v). By the transition law, out o f u ( p , q}dpdq photons, P(x, y\p, q)dxdy percentage would be transferred to (x, x + dx) x ( y , y + dy}. Therefore,

which is exactly the spreading transform of u when P is given by (2.35). Notice that so far the actual time lapse 8t of the single transition step has been ignored. In the continuous and evolutionary setting, the transition should be time dependent:

For Brownian motion [164,228] for example, one has P,+!i = Pto Ps, when Pt is understood as density transform operators as in (2.36). Furthermore, the dynamic transition law for Brownian motion is well known to be

the 2-D isotropic Gaussian with variance 2t (one t from each dimension). More generally in image and vision analysis, point spreading via random walks could be realized by the associated infinitesimal generators of parabolic PDEs, which we will now discuss.

2.5.2

Linear Filtering and Diffusion

Let us start with the digital setting of which (w/y) is a given digital image defined on the canonical Cartesian lattice (/, j ) e Z2. In the linear filtering theory, one has

where e > 0 is a constant weight, and /z's are filter coefficients. The notation (k, /) ~ (/, 7) refers to all pixels (k, /) that are "connected" to (/, j). Connectivity is usually established by specifying neighboring windows, or edges in graph theory [87]. To respect the conservation law, the filter h must be highpass in the sense of

Consider, for example, the 5-pixel stencil

and the associated Laplacian filter

68

Chapter 2. Some Modern Image Analysis Tools

Figure 2.11. An example of linear anisotropic diffusion by (2.38) with diagonal diffusivity matrix D = diag(D x , D v ) and Dy : Dx — 10 : 1. The image thus diffuses much faster along the vertical y-direction. Then one has

Introducing spatial step AJC = Ay and time step Ar, we choose e = DAr/(A.v) 2 in such a way that D is a fixed constant of order 1. Then

Letting Ar and Ax = A>' pass towards 0, one ends up with the standard heat equation

provided that Ujj = u(i AJT, j Ay, t) and u\j = u(i AJC, j Ay, t + Ar). Thus the iteration of this classical linear filter is tantamount to the diffusion process. More generally, by allowing D = D(;t, y) to vary spatially, or even to be a 2-by-2 positive definite matrix D = [a(x, y), b(x, y); b(x, y), c(x, y)], one obtains the PDE version of spatially varying linear stationary filters:

These equations are accompanied by suitable boundary conditions as well as the initial condition u(x, y, 0) = UQ(X, y)—the given image. Figure 2.11 shows an example. Physically, diffusion is a typical self-driven evolutionary process, and the driving "force" is spatial inhomogeneity. It is also manifest in the flux formula

As long as D = D(x, y) = [a, b\ b, d] is positive definite, such flux must be gradient descent, since

2.5. Linear and Nonlinear Filtering and Diffusion

69

Therefore, diffusion discourages the development of large gradients and diminishes fast local oscillations caused by random noise. Furthermore, during a linear diffusion process, oscillations with different spatial wave numbers decay in different rates. For simplicity, first consider the case when £2 = R 2 and D — [a,b\ b, c] is a positive definite constant matrix. Denote by L the spatial linear operator V • D V. A monochromatic oscillation in the direction of 9 with spatial frequency k is given by0O) = ei(k-*\ with A; = (kcosO, k sin6») and x = (x,y). Then

with A = kDkT > 0. Thus 0 is L-invariant since L0 is still in the "direction" of , or loosely speaking, 0 is the eigenvector of L associated with the eigenvalue X. Invariance allows one to restrict the diffusion process ut = Lu within the 0-space and to try a solution in the form of u(x, t) = g(t)(/)(x). Then

which shows that the decaying rate for 0 = el (k ' x) is X — kDk7. Since D is a fixed positive definite matrix, one has A = O(k2}, confirming the earlier assertion that faster oscillations die out faster during the diffusion process. The above analysis beautifully carries on to linear diffusion processes on general bounded image domains with proper boundary conditions. The starting point is to first solve the invariant equations (2.39): L0 = A0 for suitable A's to characterize different oscillatory modes on the domain. It is an eigenvalue problem since L, together with boundary conditions, behaves very much like a nonnegative matrix. The squared frequencies are generally quantized to and their associated eigenfunctions depict the associated spatial patterns. Then the decay law (2.40) still holds. Following this line, one observes the striking feature of diffusion: asr —>• oo, all oscillatory components with nonzero A's must vanish. Under the Neumann adiabatic boundary condition for instance, only the direct-current component <j>\ = 1 associated with X\ = 0 remains in the end. Furthermore, by the conservation law (guaranteed by the adiabatic condition), one must have

It is of course unwanted in image processing since all image features are wiped out. Therefore in both practice and theory, it is a crucial problem to decide when to turn off the diffusion process before all major image features and visual patterns are gone. This time Ts is called the optimal stopping time. The answer clearly depends on the task at hand. For the removal of white noise with variance a2 for example, it could be shown that Ts = O(o~2). For spatially inhomogeneous noises, one could imagine that a uniform stopping time Ts cannot be sufficient to fully characterize the spatial inhomogeneity, unless the diffusivity coefficient D has already been adapted to the noise.

70

Chapter 2. Some Modern Image Analysis Tools

Figure 2.12. The median of a 5-component segment from an imaginary- l-D signal (xn).

2.5.3

Nonlinear Filtering and Diffusion

Compared with most acoustic or sound signals, images differ in that they are discontinuous functions. The discontinuities in 2-D images are commonly associated with the boundaries of objects in the 3-D world, and therefore intrinsic and visually important. Effective filtering and diffusion processes must be able to distinguish such singularities from noisy oscillations. Linear diffusions fail on this aspect due to the lack of adaptivity to given images. Thus nonlinearity is naturally demanded in image analysis and processing. The most well known and simplest nonlinear filter is perhaps the median filter:

Here as in the preceding subsections, interpixel connectivity ~ is often specified by a local window, for example, the 5-by-5 Cartesian neighborhood. Figure 2.12 shows a 1 -D example with a window of size 5. Median filtering is a special example of the more general class of order-statistics filtering. Recall that for 2n + 1 observations (JC_,M . . . , XQ, ..., xn) (a window at index 0), the order statistics x(-n) < x(-n+\) < • < x(n) are established by sorting the data in increasing order. Then a general lowpass order statistic filter H is specified by

for some suitable set of filter coefficients /z's. Notice that nonlinearity is caused by the association of the fixed filter coefficients with the order statistics, instead of with the original raw data. (The choice of XQ is purely for the sake of illustration, and in principle one could use the same expression to update the value at any fixed index m between —n and n.) For median filtering, hk = 8k, the discrete Dirac delta sequence. Due to its statistical nature (i.e., inside a window with enough data samples), median filter easily applies to noise removal. But why is it capable of preserving edges? As an example consider a simple l-D "cliff" or Heaviside-type signal: Assume that XQ — 1, where the cliff is located. Applying the median filter with a symmetric window of length 5 (i.e., n — 2), one obtains the output sequence near the cliff,

2.5. Linear and Nonlinear Filtering and Diffusion

71

Figure 2.13. An example of using the median filter (with 7 x 7 centered square window) to denoise an image with severe salt-and-pepper noise. Notice the outstanding feature of median filtering: the edges in the restored image are not blurry. Although the values at m = 1 and 2 seem to be updated unwisely for such a piecewise smooth signal, the steep cliff at m = 0 is indeed well preserved. Figure 2.13 shows an example of image denoising via median filtering. Now consider a general ideal 2-D Heaviside cliff defined by

where a,b,c are constants and H = H(s) is the Heaviside function: H(s) = 1, s > 0, and H(s) = 0 otherwise. The cliff spreads out along the straight line ax + by + c — 0. For any pixel (jt+, v + ) so that ax+ + by+ + c > 0, any median filter with a radially symmetric window (i.e., a disk) would lead to

simply because more pixels inside the disk window carry the value 1 than 0. Similarly, for any pixel (jc_, j_) with ax- + by- + c < 0,

Therefore, the output from median filtering preserves the cliff in this highly idealized case. Nonlinear filtering has become a very important research topic in digital signal and image processing [63]. We now turn to the class of nonlinear digital filters which are associated with nonlinear diffusions. In order to macroscopically preserve edges, the randomly walking particles must not go across edge "cliffs." Thus in the density transition formula

the transition probabilities Pf- w 's have to be adaptive. For example, considering u as an energy function, as motivated by Gibbs' distribution law in the preceding section, one may define

72

Chapter 2. Some Modern Image Analysis Tools

Here the controlling parameter e plays the role of temperature, and the partition function Z£ normalizes the total outgoing transition probability to 1:

The power parameter a is also tunable. By this law, transitions between two neighboring pixels with huge difference in u are highly unlikely. To derive the macroscopic diffusion PDEs, assume that HJJ = u(ih, jh} with A.X = Ay = h and that the lowpass condition holds:

(P satisfying both (2.44) and (2.45) is said to be doubly stochastic.) Then the above density transition formula becomes

Suppose the neighboring connectivity is specified by 5-pixel stencils as in the linear case before, and the transition probabilities bear the nonlinear form of

for some suitable continuous function PE(u, V«). To be precise, we now adopt the notation of Pf^ j+,_ for such Pfj kl. Then

Here the two difference-differential approximations are both accurate to the leading orders on h. We now choose P£ so that

By identifying w, 7 = u(ih, jh, t) and «,y = u(ih, jh, t + e), we therefore obtain the nonlinear diffusion equation in the limit e, h —>• 0:

Notice that to make the lowpass condition (2.45) valid, it suffices to require PE < 1/4, or £\\D\\oc/h2 < 1/4. In computational PDEs, this is well known to be the CFL condition named after Courant, Friedrichs, and Lewy [287]. In [251], Perona and Malik first proposed to properly choose the dependence of D on Vw in order to preserve edges in images.

2.6. Wavelets and Multiresolution Analysis

73

Similarly in the preceding example (2.43) inspired by Gibbs' distribution law, take e = h a / f i for some fixed parameter ft of order 1. Then to the leading order of h as h —» 0,

Define

with At denoting the time lapse to be associated with a single-step transition of the process. Let At, h -> 0 in such a way so that D\ and D^ both converge. Then we obtain the nonlinear anisotropic diffusion

with D = diag(D), D2) being a diagonal matrix. Furthermore, as h -» 0, the partition function Ze converges to

Thus if we let /i, Af -»• 0 so that h2/ At ->• 3, then

The artificial choice of 3 leads to the simple yet tight estimation: D\ < I and D2 < 1.

2.6

Wavelets and Multiresolution Analysis

In this section, we intend to give a concise introduction to the classical wavelet theory [96, 204, 215, 290]. The most recent developments on geometric and nonlinear wavelets can be read in, e.g., [41, 107, 108, 155, 248].

2.6.1

Quest for New Image Analysis Tools

Unlike many other signals, images result from the random positioning of individual objects in the 3-D world. The most striking features of these objects are manifest in the large dynamic range of spatial scales and their easily recognizable (to human vision) differences in terms of surface geometry, reflectance, and various patterns of surface textures. Being able to easily decorrelate and identify individual image objects and their patterns is thus demanded by image and vision analysis. Wavelets and their associated techniques have so far met such requirements most satisfactorily. An individual wavelet, literally speaking, is a localized small wave. It acts as a virtual neuron that fires strongly whenever localized visual features are presented to it. Due to localization, it can respond strongly only when its window captures the target feature within.

74

Chapter 2. Some Modern Image Analysis Tools

Wavelet analysis studies how to design, organize, and analyze such wavelets and achieve efficient computational schemes. A significant portion of its mission is to scientifically model human and machine vision and to effectively perform various image processing tasks. To better understand and appreciate the characteristics and advantages of wavelets, let us first briefly return to the analysis of conventional signals that are generated by typical physical systems. Near their equilibrium states, such systems could usually be well approximated by linear systems. An input / and output u are then connected by a linear system A: For any system such as a string of a guitar, the vocal channel of a human being, or the surface of a drum, the spatial part of A could often be well modelled by a second order elliptic operator where x and t denote the spatial and temporal variables, and D and b specify the physical nature of the target system. In practice, suitable boundary conditions have to be specified as well. If the nature of the system is independent of time t: D = D(x} and b — b(x), its intrinsic states or eigenmodes are L-invariant and specified by the eigenvalue problem

where the eigenvalue X reflects the energy level of the associated state. Generically, the system possesses a collection of eigenmodes,

The general theory of second order linear elliptic operators [117, 132] claims that X,, —>• oo as n —>• oo, and (0n);, provides a natural linear basis for studying signal generation and evolution under the linear system A, which is typically connected to 3, — L if A is purely diffusive or dn — L if A is purely convective. For instance, consider the small-amplitude vibration of a string with length a and two ends fastened. Its eigenmodes are given by

Then it is well known that

so that (4>n)n is an orthonormal basis of harmonic waves for L 2 (0, a). As is clear from the expressions, these harmonic modes neatly combine both the motion law (of strings) and the global geometry (i.e., size a), and therefore cannot be local. The nonlocality can be further made explicit by considering the signal diffusion problem

2.6. Wavelets and Multiresolution Analysis

75

with the initial shape given by Dirac's delta u(x,0) = f = 8(x - a/2), an idealized local impulse at the middle point. The signal diffusion can be expressed by u(x, t) = . Notice that the initial impulse equally contaminates all the odd eigenmodes since

Spatial localization is therefore not intrinsic for such systems. One has to turn to wavelets for localized signal analysis or processing.

2.6.2

Early Edge Theory and Marr's Wavelets

Edge is the most common and significant visual feature in images, which defines and differentiates objects and leaves crucial cues for 3-D spatial orders [234]. Consequently, it is one of the most fundamental problems in image and vision analysis to properly define and extract edges from 2-D images, as initiated by Marr [210], and Marr and Hildreth [211] and further advanced by many others [44, 130, 165, 226]. Marr defined edges of an image u to be the locations where the zerocrossings of the Laplacian of the image u mollified to the scale of a. More specifically, let x = (x1, X2) and r = x |, and let

so that Ga is a radially symmetric probability density with variance 2„)„ are "hard cutters,"

i.e., nonsmooth indicator functions associated with a partition of the whole space, one could easily show that indeed redundancy is completely gotten rid of since (^n.k/V^n)n.k is an orthonormal basis of L2(R). But such cheap orthonormality comes at the heavy cost of approximation precision: for a generic smooth image u ( x ) , a simple nonconstant linear function for example, the approximation error decays at a sluggish rate of O(\k\~l). Therefore, much "softer" or smoother windows are preferred for securing high order accuracy in image analysis and synthesis. The key questions are, Using smooth windows, can one indeed design an orthonormal basis that is very similar to (^n.k)n,k^ and In which aspects can V^./t's be further polished? The answer to the feasibility part is affirmative. Among many of the existing designs, we now give a brief introduction to the remarkable approach of Malvar [205] and Wilson [321], which leads to what are called Malvar-Wilson wavelets and their variants [91,215]. First, choose a smooth window template w(x) which satisfies (see Figure 2.15): 1. w is supported on [—TT, 3n] and symmetric around x = TT: w(2jr — x) = w(x); 2. restricted on [—TT, n], w(x) satisfies the nonaliasing condition w2(x) + w2(—x) = 1. Generate a sequence of windows by 2n -translations: wn(x) = w(x — 2n7t), n — 0, ±1, ±2, . . . . Then the symmetry and nonaliasing condition lead to

Notice that the infinite sum consists only of at most two terms since supp wn and supp wm never overlap for \n — m\ > 1. If the old construction (2.54) were to be followed, one would have to utilize e'5 x , k e Z, to be associated with each window wn(x) since /„ = /o = 4;r, or equivalently,

2.6. Wavelets and Multiresolution Analysis

79

Figure 2.15. A generic approach to designing the window template w(x) via symmetrically constructing the profile of its square w2(x). But two adjacent windows wn and wn+\ do clearly overlap and inevitably introduce redundancy, which can, however, be cleverly removed by an alternated "downsampling" of factor 2. That is, associated with any even window W2n(x), one applies only half from (2.56),

while for any odd window W2n+\, one applies the other half,

Therefore, we have derived the set of so-called Malvar-Wilson wavelets:

where CQ — 1 and Q = V?-, k > 1 for L2 normalization. The beautiful fact is that the Malvar-Wilson wavelets (Vv*)«,* consist into an orthonormal basis in L2(R)\ The interorthogonality between two adjacent wavelets

results from the symmetry (or evenness) of the template w(x) and the oddness of the product of a cosine with a sine, while the intraorthogonality between two wavelets within a common

80

Chapter 2. Some Modern Image Analysis Tools

window

is guaranteed by the nonaliasing condition of the template as well as the fact that both

are orthonormal bases for L2[0, 2n}. The very last statement follows from the fact that these cosines are exactly the eigenmodes of the Sturm-Liouville system

while the sines are the eigenmodes of

This alternated cosine-sine association was first suggested by Nobel laureate Kenneth Wilson in the study of renormalization group theory [321], was later systematically developed in the remarkable paper by Daubechies, Jaffard, and Journe [97], and was also independently rediscovered in digital signal processing by Malvar [205]. The Malvar-Wilson wavelets could be extended into a broad class of orthonormal space-frequency wavelets, as done by Coifman and Meyer [91]. For example, one could replace the translation invariance of the windows

by dilation invariance

Such windows have variable window lengths, and accordingly the associated cosine and sine waves should be scaled proportionally. Much more generally, one could completely abandon a single-window template and the translation- or scaling-invariant property and individually design each window wn. As long as one 1. maintains the nonaliasing condition for intraorthogonality and 2. well coordinates the symmetry behavior of any pair of adjacent windows and their associated harmonic waves in the overlapping region for interorthogonality, an orthonormal basis of space-frequency wavelets is still easily achievable by this design.

2.6.5

The Framework of Multiresolution Analysis (MRA)

From the designs of Gabor to Malvar-Wilson, wavelets still live in the mighty shadow of Fourier frequency analysis. Historically, the situation did not make an energetic turn until Meyer and Mallat introduced the independent and general framework of multiresolution analysis (MRA). Though mainly originated from methodologies in signal and image

2.6. Wavelets and Multiresolution Analysis

81

processing, MRA does not look completely strange to scientists working on multiscale phenomena such as multiscale turbulence flows and multigrid methods. The contributions of the MRA framework to the development of wavelets have been extraordinary for a number of crucial reasons. First, woven into the core principles of MRA are a few ingenious novel ideas in modern signal and image processing. Second, MRA provides a systematic approach for wavelet construction and analysis. More remarkably, it allows fast and efficient digital implementations for image analysis and synthesis by wavelets. Finally and above all, it lifts wavelet theory out of the shadow of Fourier frequency analysis and puts it on the solid and self-contained foundation of multiscale analysis. MRA is an elegant mathematical framework for representing and analyzing images and general signals in multiple scales. An orthonormal MRA of L2(R) is an ordered chain of closed subpspaces: which satisfies the following three conditions: 1. [Completeness] lim^oo Vj = L2(R), andlim^-oo Vj = {0}. 2. [Dyadic Similarity] u(x) e Vj u(2x) e Vj+\. 3. [Translation Seed] There exists a function 0 e VQ, so that (0(jc — k))k is an orthonormal basis of VQ. Thus if such an MRA indeed exists, the entire L2(R) world can be coded by a single "seed" function 0. In fact, define

By dyadic similarity, (0/,fc)fc is also an orthonormal basis of Vj. By the completeness postulation, any image u(x) e L2(R) can be well approximated to any desired precision by its projection HJ = PjU onto Vj:

In this sense, in the early days of wavelet theory, 0 was often called the "father" wavelet. It is now usually called the scaling function, or the shape function by some computational scientists. As a necessary condition, the scaling function must satisfy the two-scale similarity condition K,

ft,

with the coefficients (/z/Ofc given by

In practice, (2.59) provides the equation for designing the scaling function 0 (x) from appropriate coefficients (/z^)/t and is often called the two-scale relation, the refinement equation, or the dilation equation.

82

Chapter 2. Some Modern Image Analysis Tools

Another characteristic condition on the scaling function is that fK ^ 0, assuming that 0 e L2(R) (~}Ll(R). To avoid unnecessary technicality, assume for simplicity that is compactly supported. (Generally, a moderate decay rate such as = O ( \ x \ ~ l ~ f ) at x = ±00, combined with tools such as the Lebesgue dominated convergence theorem, will still be able to accomplish the following proof.) Suppose otherwise fR

-f oo. Since all Oy) y >o can be easily shown to be uniformly compactly supported, one must have, as j -> oo,

which contradicts (2.60). This seemingly trivial condition has profound influences. First, it allows one to normalize the scaling function by fR 0 = 1. Then the dilation equation (2.59) implies that

which is called the lowpass condition in digital signal processing when h = (hk)k is considered as a digital filter. Second, applying Fourier transform to the dilation equation (2.59) leads to where H(a)) = £^ hke~lk(a is the impulse response of the filter h. In combination with 0(0) = 1, iteration of (2.61) gives the explicit expression of the scaling function in the frequency domain:

Thus the scaling function is completely determined by the lowpass filter h. Consequently, it can be expected that the orthonormality condition on ((x — k)k) are orthonormal, one has

Still denote the inner product in L2(0, 2n} by {• , •}. The last identity is then equivalent to

which immediately implies that

By the general identity (2.64), we have thus derived a necessary condition on the lowpass filter h = (/i*)* in an orthonormal MRA:

The lowpass condition //(O) = 1 then immediately implies that H(n) = 0, and the digital filter indeed suppresses high frequencies. In practice, identity (2.65) provides a good starting point of designing orthonormal MRA [96]. Thus far we have only discussed the scaling function 0 and its associated lowpass filter h = (h^k- The real beauty of MRA of course cannot terminate here since wavelets have not yet emerged from the analysis. In MRA, wavelets correspond to the details that are lost during the transfers from fine scales to coarser ones. Due to the dyadic scale similarity in MRA, it is sufficient to focus on a representative pair of adjacent scales: V0 c V\. V\ has finer details, and the orthogonal projection of an image u\ e V\ onto UQ E VQ: UQ = PQU\ is to wipe out some details of MI that are impossible to detect in VQ. Let W0 denote the range of (/ — PQ) \ v , or equivalently, the space of "details." Then

Generally, let Wj denote the space of details at scale A = 2 ;, i.e., V}+i = Vj © Wj. Then it is easy to see that the dyadic scale similarity of MRA is faithfully inherited:

84

Chapter 2. Some Modern Image Analysis Tools

Figure 2.16. MRA as (Hilbert) space decompositions: the finer resolution space ¥2 is decomposed to the detail (or wavelet) space W\ and coarser resolution space V\. The same process applies to V\ and all other Vj 's [290]. The completeness of MRA is now read as

where J is any arbitrary reference scale level, and the two infinite sums always denote the smallest closed subspaces of L2(R) that contain all the relevant wavelet spaces. The formula shows that a general image u(x) e L2(R) is the accumulated effect of its details at all scales. Figure 2.16 gives a low-dimensional visualization to the relations among different subspaces. Introduce the dilation notation for a given subspace U c L2(R),

Then (2.67) could also be written as

which clearly shows that WQ completely determines the representation and analysis of images in L2(R). To further reveal the structure of WQ, consider a general element 77(jc) e WQ. By (2.66), one must have

or, in the frequency domain,

The set of coefficients g = (gk)k have to reflect the nature of "details," i.e., WQ _L V0, which amounts to (77 , 0(;c + fc)) = 0, k e Z, or in the frequency domain,

2.6. Wavelets and Multiresolution Analysis

85

After 27r-periodization, it leads to

Having co replaced by 2co, and applying the two-scale relations (2.61) and (2.69), we obtain

where A is as given in (2.63). Since A(w) = 1 due to the orthonormality of (0(jc — k ) ) k , we have Define A.(<w) = G(a>)elLa}/H(co + TT) for a fixed odd integer L. Then the last equation becomes implying that X(o)) must be a n-periodic function, or C(&>) = X(a>/2) ITC-periodic. Define a template Then G(a>) — Go(tt>)C(2<w), and

where ^ = (Go(f )0(f )) v , defined by the template GO, is called the mother wavelet. Assume the Fourier coefficients of C(&>) are (Q)*- Then it has been established that a general detail element rj e WQ must be in the form of

Therefore, the detail space WQ is spanned by all the integer translates of the mother wavelet VK*). We further show that (ijs(x — k))^ is in fact an orthonormal basis of WQ. As for the scaling function, define

Then it is adequate to show that B(a>) = 1. Similar to the scaling function, the two-scale relation (2.69) easily leads to

Here we have applied A(w) = 1 and (2.65). By scale similarity,

86

Chapter 2. Some Modern Image Analysis Tools

is an orthonormal basis for L2(R), in which sense that ^ has been fondly called the mother wavelet. It satisfies the so-called wavelet equation

for a set of suitable coefficients g = (gk) as discussed above. Let G(a>) denote the Fourier transform of g. Then the necessary condition (2.70), combined with //(O) = 1 andH(n) = 0,leads to implying that fR -^ = 0. These conditions on G and V are called the highpass conditions of wavelets, since they both filter out low frequency components near co % 0. In summary, an orthonormal MRA carries four pieces of crucial data: the scaling function 0 and its lowpass filter h connected by the dilation equation

and the mother wavelet \js and its highpass filter g connected to the scaling function by the wavelet equation A.

For orthonormality, the two filters must satisfy

The actual design of orthonormal MRA starts from these last two equations on the lowpass filter h and highpass filter g. In most applications in signal and image processing both h and g contain only a finite number of nonzero coefficients; then it essentially becomes a polynomial (or Laurent polynomial) design problem [95, 96, 290]. Such filters are said to have finite impulse responses (FIR) in digital signal processing [237]. Figure 2.17 shows a pair of scaling function and mother wavelet.

2.6.6

Fast Image Analysis and Synthesis via Filter Banks

We now further explain how the MRA framework leads to efficient computational algorithms for multiscale image analysis and synthesis. By the completeness postulation, any image u(x) G L2(R) could be approximated to any precision by its projections HJ = PjU onto increased resolutions V/:

The efficiency of MRA very much results from the intrinsic connections of such expansions at different scales, as nourished by the dyadic similarity postulation in the design of MRA. Thus in situations where multiple users are interested in different tasks and scales, there is a universal data storage and retrieval structure which could efficiently meet most of the

2.6. Wavelets and Multiresolution Analysis

87

Figure 2.17. A pair of compactly supported scaling function (f)(x) and mother wavelet iff(x) by Daubechies' design [96]. demands. In particular, it becomes unnecessary to work out the expensive integrals for all j and k in (2.76). Such efficiency results from the two fundamental equations (2.72) and (2.73). By the dyadic similarity postulation, one has

Suppose that Uj+\ = Uj + Wj e Vy 0 Wj. Then

Similarly for the wavelet coefficients,

Define the transposes of the two (real-valued) filters by spatial reversal:

and as in digital signal processing also define the downsampling operator in I2:

by dumping all the odd components. Then the above results can be concisely expressed by

88

Chapter 2. Some Modern Image Analysis Tools

Figure 2.18. Fast wavelet transform via filter banks: the two-channel analysis (or decomposition) and synthesis (or reconstruction) banks. where the asterisk symbols denote discrete convolutions of sequences. It is exactly here that the mathematical framework of MRA neatly meets modern digital signal processing (DSP). In DSP, these two identities lead to a two-channel analysis filter bank (h, g), which takes cj+[ as the input and (cj, d j ) as the output. The iteration of this simple DSP unit leads to fast wavelet decomposition: starting from any fine scale X. = 2~J,

Therefore, once the representation at a fine scale is available, all the coarser-scale decompositions can be obtained by simply cascading the two-channel analysis bank, which spares the costly computation of all the inner product integrals. The analysis bank (in Figure 2.18) is accompanied by its companion synthesis or reconstruction bank. In the orthonormal MRA framework, the synthesis bank is the "transpose" of the analysis bank, very much like the well-known linear algebra fact that inverting an orthogonal matrix amounts to transposing it. The transpose of the analysis filter pair (h, g) is simply (h, g), the original filter pair in the dilation and wavelet equations. The transpose of the downsampling operator (4 2), as also well known in multigrid methods [34], is the upsampling operator (t 2):

Then the synthesis formula is given by (see Figure 2.18)

In combination, MRA goes completely digital, and in between the analysis bank and the synthesis bank in Figure 2.18, i.e., on the wavelet domain, lies the vast freedom for image coding, transmission, processing, or analysis. The filter bank structure also immediately suggests the way for designing the socalled biorthogonal MRA in which (VO',/0./,* is onty a Riesz basis for L 2 (7?) and (0/,*)t for Vj. Digitally, it means that the synthesis filter pair (/i y , gs) is no longer required to be the transpose of the analysis pair (ha, ga). Instead, a biorthogonal MRA imposes only the condition of perfect or lossless reconstruction:

2.6. Wavelets and Multi resolution Analysis

89

Of course, like computing the inverse of a general affine matrix A in linear algebra, one has to make sure that filter banks are well conditioned, so that both the analysis and synthesis steps are stable. The notion of biorthogonality greatly expands the family of MRA-based wavelets. And even more remarkably, it allows the scaling functions, wavelets, and their filters to be symmetric or to have even phases [290], a crucial technical virtue in signal and especially image processing.

This page intentionally left blank

Chapter 3

Image Modeling and Representation

The goal of image modeling or representation is to find proper ways to mathematically describe and analyze images. It is therefore the most fundamental step in image processing. One, however, must realize that there is no absolutely the best representation, since optimality inevitably depends upon specific processing tasks, as in the case of different representations of natural numbers: the decimal system is more convenient than the dyadic one in daily life, but the latter is more natural for digital or quantum computers. The current chapter introduces five general and useful approaches to image representation, based upon which many successful image processors are to be developed in later chapters.

3.1

Modeling and Representation: What, Why, and How

Digital images are most commonly presented as a matrix of scalars for gray-scale images or vectors for color images, since they are often captured by charge coupled device (CCD) arrays (as in digital cameras) or displayed by liquid crystal arrays (as for laptop or pocket computers). Pixel-matrix representation is, however, by no means the most efficient from the information-theoretic point of view. In what follows, we shall call such direct matrix representation u = (M/J) or its analog idealization u = u(x, y) with (x, y) e £2 = (a, b) x (c, d), the physical image. A representation of a given class U of physical images refers to a transform T under which any image u from the class is transformed to a new data type or structure w = Tu. Let W denote the range space of 7~, which is also called the transform space. Then we have, for short, A representation is said to be linear if 1. both the image class U and its transform space W are linear spaces, or convex cones of linear spaces, so that for example, au \ + bu^ belongs to U as long as u \ and u^ do and a, b e R or R+; and 2. the transform is linear.

91

92

Chapter3. Image Modeling and Representation

For instance, consider the class of all two-pixel images: U = {u = ( s , t ) : s.t e R+}, which is a convex cone in R2 (i.e., the first quadrant). Let the transform space be w = (a, d) € W = R2, and let T be the average-difference transform of Haar [96, 290]:

It is the key linear transform for the construction of Haar's wavelets [96, 290]. Mathematically, even the human vision system could be considered as a (biological) representation operator T/,. The photoreceptors (i.e., cones and rods) of the two retinas are indeed laid out very much in a 2-D manner. Let U denote the class of images u = (HI, ur) projected onto the retinas of both left and right eyes. Let W denote the class of electrochemical signals encoded by all the major visual cortices (e.g., VI and MT [227, 168]). Then the human vision system (from the photoreceptors, ganglion cells, and lateral geniculate nucleus (LGN) to the primary visual cortex V1 and so on) is a biological transform Th realized by this complex neural and cellular network. Numerous evidences show that it is nonlinear [150, 152]. How does one judge whether a particular type of image representation is good for image and vision analysis? Images as a special class of signals carry rich material, geometric, and positional information about the 3-D world. A good representation should be able to highlight such information and efficiently catch the associated key visual features. This is the most general guideline for efficient image modeling or representation. A representation T is said to be lossless if any u can be perfectly reconstructed from its representation w = Tu. That is, there exists another (reconstruction) transform 'R, from W to U so that A necessary condition for a representation T being lossless is being injective; namely, two distinct images u and v should be transformed into different Tu and TV, since Tu = TV would immediately imply

Haar's representation in (3.1) is apparently lossless since the reconstruction 72. is easily realized by More generally in terms of linear algebra, the lossless nature of Haar's representation is guaranteed by the invertibility of the associated matrices. In signal and image analysis, the two transforms T and 72. are usually called the analysis and synthesis transforms separately. Analysis is to analyze a given image signal and extract its key features and information, while synthesis is the effort of reconstructing the original signal from the output of the analysis step. If one has a complete code book or dictionary at hand, being injective is already sufficient to be lossless. For instance, if 10 complicated images have been a priori stored and labelled by M I , 1*2, . . . , U\Q, the simple injective representation

3.2. Deterministic Image Models

93

is already sufficient to reconstruct the images as long as the code book is available. Imagining that the number k = 3 is received, one can then simply turn to the code book to retrieve the image w 3 . Unlike any national language whose vocabulary is relatively stable and finite, it is practically impossible to create a dictionary which contains all the images of the world. As a result, the preceding dictionary-based representation idea can never work effectively in the real world. Instead, most good representations focus on catching the intrinsic visual features of images. By ignoring visually unimportant image information, a good representation is inevitably lossy instead of being lossless. However, the associated reconstruction processes often lead to reasonably good synthesis since no significant visual information has been ignored. This is indeed the case in the wavelet-based JPEG2000.

3.2

Deterministic Image Models

We first discuss several deterministic image models that prevail in the literature. They are all mathematical models that approximate real images to certain levels of faithfulness.

3.2.1

Images as Distributions (Generalized Functions)

Treating images as distributions or generalized functions is the broadest approach for deterministic modeling [225] and has profound merits in image understanding as explained below. As before, let £2 denote an open and bounded 2-D Lipschitz image domain. The set of test functions are defined as

Each test function 0 e D(£2) could be considered as a linear sensor for capturing image signals. An image u on £2, being treated as a distribution, is a linear functional on D(£2):

where the inner product symbol has been used formally, as is customary in distributional theory [3, 193]. Thus the image space is the distribution space D'(£2). Though such an image u may not look like any familiar function at all, it definitely outputs a single response (u , 1) norm, the resultant Sobolev image spaces are often denoted by Wn-p(£l) [3, 126, 193],

3.2.4

BV Images

From general distributions to Lp images, and further to Sobolev images, each time when more regularities are added, the spaces become more specific and better manageable. However, there is a tradeoff between the regularity of image spaces and the fidelity of image modeling. An optimal solution is to find the balance, namely, to construct an image space that is both mathematically tractable and practically faithful in representing the key features of general images. This is the space of BV images, as already introduced in the preceding chapter, whose role in image analysis and processing has been fundamental and influential ever since the well-known work of Rudin and Osher [257], and Rudin, Osher, and Fatemi [258]. From the distributional point of view, an L' image u belongs to B V(£2) if and only if its distributional gradient VM satisfies

where ||0||oc = sup ;cen (0f(jc) + ^(x))1/2. If so, this supremum is called the TV of u and denoted by TV[w] or \Du\(£2). BV(£2) is a Banach space under the natural norm

3.3. Wavelets and Multiscale Representation

99

The remarkable features of BV images are characterized as follows. 1. All W1'1 images (including //' as a subspace since Q, is bounded) are BV images, but the former do not allow edges (i.e., intensity jumps along curves). From the very beginning of artificial intelligence research, edges have been recognized as crucial visual cues in perception and image understanding [44, 183, 210, 226]. 2. While Lp images also permit edges, they keep no specific track of local oscillatory irregularities. But BV images do. Thus the notion of B V images indeed achieves a reasonably good balance between penalizing irregularities (often due to noise) and respecting intrinsic image features such as edges. The concept of BV images is also compatible to the occlusion-generative imaging process of Lee, Mumford, and Huang [191]. Most objects in the 3-D world are not transparent, causing the generic phenomenon of occlusion. On a 2-D image domain £2, occlusion is modelled as follows. Suppose there are N Lipschitz subdomains of £2, ordered from foreground to background: whose boundaries all have finite 1-D Hausdorff measures. Assume that each Qn has a constant gray level of un, n — 1 : N. Then the synthetic image u generated from such a scene is provided that the background intensity is 0. Then the TV is

implying that such a synthetic image must belong to BV(£2). If the w w 's are allowed to be smooth functions instead of constants, the resultant image u is piecewise smooth and still belongs to BV(fi).

3.3

Wavelets and Multiscale Representation

In this section, following the basic knowledge prepared in the preceding chapter [86, 96, 290], we further explore the role of wavelets as efficient tools for multiscale image modeling and representation.

3.3.1

Construction of 2-D Wavelets

We first introduce two different approaches for constructing 2-D wavelets from 1-D ones. Approach I: Three Mother Wavelets at Equal Scales

As in the preceding chapter, let ty (t) and 0(0 denote the mother wavelet and scaling function for an orthonormal multiresolution analysis in R 1 :

100

Chapters. Image Modeling and Representation

so that (0(r — k) | k e Z) is an orthonormal basis for Vb, and (\J/(t — k) \ k e Z) for the wavelet space WQ, which orthogonally complements V0 in the more "detailed" space V\. Then is an orthonormal wavelet basis for L2(R). To avoid complication, it shall be assumed in this section that the image domain Q = R 2 . In practice, many ways exist for extending image data on a bounded square domain to R 2 , such as zero-padding, symmetric reflection, periodic extension, extrapolation, etc (see, e.g., Strang and Nguyen [290]). There also exist rigorous ways of constructing wavelets on bounded domains (see, e.g., the celebrated work of Cohen et al. [89]). Recall that the tensor product of two linear spaces X and Y is defined by

If both X and Y are Banach spaces of infinite dimensions, the above definition should be followed by a further completion process so that X Y is Banach as well. If X and Y are both Hilbert spaces, and (xn \ n e Z) and (ym \ m e Z) are their orthonormal bases separately, then under the natural tensor-product extension of inner products, the tensor products of the bases vectors

become an orthonormal basis for X Y. The construction of 2-D MRA starts from the following relation. Theorem 3.5 (Separability). Treating R2 = R © R (direct sum of vector spaces), one has, under Lebesgue measures,

It claims that for any /(jc) = f(x\, x{) e L 2 (R 2 ), and e > 0, there exists a function fs(x) in the form of

so that || / — / e ||z. 2 (R 2 ) < £• There are many different approaches leading to this fact, and perhaps none could be more convenient than the Fourier transform. Proof. Notice that the test functions Z)(R 2 ) (i.e., smooth functions with compact supports) are dense in L 2 (R 2 ). For any test function define

where a. = (a,b) denotes the four types. Then the following theorem is straightforward.

102

Chapters. Image Modeling and Representation

Theorem 3.6 (2-D Wavelet Decomposition). 1. At any scale j, (^ k \ k e Z 2 ) is an orthonormal basis for V(j). 2. Let W(j) denote the orthogonal complement ofV(j) in V(j+\). Then

is an orthonormal basis for W (y) . 3. Thus V(i) = V(0) © W(Q) has two sets of orthonormal bases: (i/rj\' k \ k e Z2) and the combination of othornormal bases from V(o) and W(Q) just established above. They are connected by the two-scale relations

where ha 's are four digital filters, lowpassfor a — (0, 0) and highpassfor the other three. Computationally, the four 2-D digital filters ha's are also the tensor products of the two 1-D digital filters defining the scaling function and wavelet in 1-D MRA. Therefore, the filter bank realizing the above change of bases is the tensor product as well. More precisely, suppose u € V(j+\) has the representation

After passing through the 1-D analysis filter bank (i.e., lowpass and highpass channels with downsampling; see Section 2.6) along the k\-direction for any fixed &2> (c(j+\\ I ^ e ^ 2 ) splits into two intermediate fields of coefficients:

where the uniform scale index splits into (j, j + 1) since the ^-direction is still in the scale level of j + 1. Next, after passing through the same 1-D filter bank along the ^-direction for each fixed k\, the two intermediate fields further split into the desired four fields of coefficients of u in the bases {^"J'kK \ k e Z2, a. — (0|1, 0|1)}:

Example. As an example, consider the 2-D Haar wavelets:

3.3. Wavelets and Multiscale Representation

103

Figure 3.2. Three Haar mother wavelets in two dimensions via tensor products. Figure 3.2 helps visualize these mother wavelets. What looks very appealing is that the three wavelets imitate the three differential operators (under the finite difference scheme)

up to some multiplicative constants. Finally, let us go to the spatial frequency domain to understand the very nature of such 2-D MRA. Either in the case of Shannon wavelets [96, 204, 216] or in terms of the asymptotic behavior of Daubechies's family of wavelets [280, 281], MRA clearly parallels Littlewood and Paley's theory on dyadic partitioning of the Fourier domain [127]. Thus in the 1-D case, each resolution space Vj can be approximately thought of representing all signals bandlimited to \a)\ < 2'jr and its orthogonal complement Wj in Vj+i all signals bandlimited to the high frequency band \a>\ e jr(2 y , 2 7+1 ]. In fact, it is the primary change of notion in wavelet theory to think of these bands as scale resolutions instead of individual frequencies [96, 204, 215]. This view naturally extends to the 2-D MRA just constructed above. The mechanism of tensor product implies that each resolution V(j} approximately represents all 2-D images bandlimited to the dyadic cube

and its wavelet complement W(7) denotes the high frequency Cartesian ring inside Qi+\ but outside QJ. Each wavelet space W (j) further orthogonally splits into three subspaces represented by the three mother wavelets

corresponding to the further partitioning of the Cartesian ring into three bands:

104

Chapters. Image Modeling and Representation

Approach II: Single Mother Wavelet and Anisotropic Scale Mixing One characteristic feature of the above construction is that each space V{J) exactly corresponds to a single resolution level j in both directions of jci and x2. There exists, however, another way of tensor-product-based construction that allows anisotropic scale mixing. In addition to the translation vector k = (k\, k2), define j = ( j \ , j2) e 1r to be the anisotropic resolution vector. According to the general construction formula (3.7), the tensor products consist into an orthonormal basis for L 2 (R 2 ). To further investigate the intrinsic structures of this family of 2-D wavelets, for each anisotropic resolution vector j = ( j \ , j2), define the diagonal matrix

Then from linear algebra,

Treating the continuous pixel variable x e R2 as a column vector (x \, x2)T and k = (k\,k2)T as well, we then have for each j e Z2,

which looks more like the 1 -D formula. For each anisotropic resolution vector j = ( j \ . j2), define its associated wavelet subspace Following the discussion in the preceding subsection, in the Fourier domain, Wj corresponds to the four rectangles of high frequencies:

If \j\ ~ J2\ ^> 1» Qj represents a rectangle with large aspect ratio, resolving image features that are relatively highly oscillatory in one direction while smoother in the other. Thus, one could imagine that this second construction is particularly efficient in representing images structures such as thin silky cloud, hair or fur, tall grass, or straws. In the physical domain, such structures must also bear large aspect ratios according to Heisenberg's uncertainty principle [288J.

3.3.2

Wavelet Responses to Typical Image Features

In this section, we discuss wavelet responses to typical image features: smoothness and edges. To illustrate the main ideas, we shall restrict the analysis only to the 1-D case.

3.3. Wavelets and Multiscale Representation

105

Response to Smooth Images A mother wavelet VOO € L 2 (E) is said to be r-regular if it belongs to C r (M). For orthonormal MRA, there is a beautiful theorem connecting the smoothness of the mother wavelet to the number of its vanishing moments, belonging to Meyer [216], Battle [18], and Daubechies [96]. An improvement on Daubechies' result was also made by Cai and Shen [38]. Theorem 3.7 (Vanishing Moments [96,38]). Suppose the mother wavelet iff (x) is r-regular. In addition, (i) T/T has up to rth order moments, and (ii) 1^(^)1 < C(\ + \x\Y, x e R, for some fixed constant C. Then all the moments vanish up to the rth order:

See Daubechies [96] and Cai and Shen [38] for a proof. Notice that conditions (i) and (ii) are trivially satisfied for compactly supported wavelets, which are the most common in applications [95]. In particular, we have the following results. Theorem 3.8 (Decaying Rate and Smoothness of Wavelets). (1) A compactly supported mother wavelet \fs cannot be C°°. (2) More generally, a mother wavelet \jf that decays at certain exponential rate (i.e., O(e~a\x\}for some a > 0 at x = ooj cannot be C°°. (1) results from (2), while the proof of (2) is straightforward from the preceding theorem; otherwise, all moments of \js would vanish, or equivalently, all the Taylor coefficients of \}r (Fourier transform) vanish at the origin, which is impossible as ty is a nonzero function and analytic on a horizontal stripe containing the real axis (due to the exponential decay). In what follows, we shall assume that the mother wavelet \fr is r-regular and compactly supported on some finite interval [a, b}. As for Daubechies' wavelets [96, 95], r could be any desired integer order. Suppose u(x) is any given Cm 1-D image with 0 < m < r + 1. We expect the wavelet coefficients to be small in small scales (i.e., large j) when the signal behaves smoothly relative to the tiny wavelets. By definition, the wavelet coefficient GJ^ is

where XQ = Xj^ denotes 2~ ; fc. By Taylor's expansion at JCQ,

For convenience assume that Cm = Cm(u) = ||M (m) ||oo < °o. Then the residual Rm can be easily bounded:

106

Chapters. Image Modeling and Representation

On the other hand, since r > ra — 1, the inner product of ^00 and Pm-\(2~jy) must vanish. Therefore,

Define the mth order absolute moment Mm of iff by

Then we have established the estimation for the wavelet coefficients

where as in computational PDEs, hj denotes the scale 2~-/ at level j. Let Dj denote the orthogonal projection from L2 to the wavelet space Wj, and Wj(x) = D j u ( x ) . Due to the limited overlapping among V0'.t's at each fixed resolution j, the estimation in (3.20) immediately implies that

Notice that the extra factor hi bounding cy^'s in (3.20) has been cancelled out by the normalization factor h~- in iffj.k(x) = h^l/ if/(h~](x — */.*))• In both formulae (3.20) and (3.21), #(•) can be improved to o(-) after some small modification of the argument [38]. Besides this classical or pointwise approach, estimations such as (3.21) can be established for Sobolev images in even more convenient ways, which we now explain. In terms of its Fourier transform u(co) [193], the condition u e //"'(R) is equivalent to Assume that u has a power law decay

for some optimal a. Then (3.22) implies that (in one dimension)

Suppose the regularity r of the mother wavelet if/ is much larger than m. Then in the Fourier domain, the jth level wavelet component Wj of u is approximately equivalent to (with a. — m — 1/2 denoted by £ > 0)

where in reality the indicator function should be replaced by a smooth function (depending on the construction of the associated MRA [96, 216]). Therefore, for some constant C,

Thus from the approximation theory point of view, for smooth images, high resolution wavelet coefficients and components can all be set to 0.

3.3. Wavelets and Multiscale Representation

107

Response to Image Edges Consider a normalized ideal edge of a 1-D Heaviside image

where xe is the edge or jump pixel. As in the preceding section, write Xj^ = 2~^k and hj = 2~j, and assume that suppi/f = [ a , b ] . Then the wavelet coefficient is

If the edge pixel xe j,k+i = 1, there are at most N = \b — a\ + 1 number of nonzero wavelet coefficients at each resolution level j. Notice that N is independent of j. Thus on the space-scale plane (i.e., x versus — Iog2 h plane with x = 2~jk and h = 2~j), the influence domain of such an ideal edge bears a parasol shape [96, 290]. More generally, if the 1-D image signal u is piecewise smooth and has many edge points, then each edge point carries its own parasol influence domain, and they overlap at coarse resolutions. Within each parasol, the wavelet coefficients are bounded by

while all the rest are still subject to the estimation developed in the preceding section.

3.3.3

Besov Images and Sparse Wavelet Representation

As a multiscale tool, wavelets are particularly powerful for studying a class of images known as Besov images, whose multiscale nature is intrinsic from their definitions. A Besov class on M. is often denoted by B"(LP), with three indices a, p, and q: 1. a is called the regularity index (or power), measuring the degree of smoothness. 2. Lp = L P (R) denotes the intrascale metric controlling the finite differences at each scale. We shall call p the intrascale index. 3. q or Lq (dh/h, R + ) denotes the interscale metric (on the scale space h e R+ with the logarithmic measure dh/h) that controls the overall regularity across all scales. We thus call q the interscale index. In what follows, we shall start out with the classical definition of Besov spaces and then state the equivalent one under wavelet representations. For simplicity, we focus on the 1-D case when the image domain £2 = R.

108

Chapter 3. Image Modeling and Representation

Besov Images and Multiscale Characterization It shall be assumed throughout this section that the regularity index a. e [0, 1) and u(x) e LP(R). For any scale h > 0, define the /^-modulus of continuity of u to be

It is uniformly bounded by 2\\u \\LP due to the triangular formula of Lp norms. An image u is said to belong to the Besov class B"(LP) if

If so, the homogeneous Besov (semi-)norm of B*(LP) is defined to be

and the true norm of B^(LP) is defined by

As in the classical Lp theory, if either p or q is oo, the corresponding norm is understood in the sense of essential supremum [126, 193]. Figure 3.3 helps visualize the meaning of the three parameters a, p, and q in the definition of Besov norms.

Figure 3.3. Besov norms in B" (Lp) 's measure the strength of signals in the spacescale plane: Lp for intrascale variations, Lq for inter- or cross-scale variations (in terms ofdX = dh/h), while h~a for comparison with Holder continuity.

3.3. Wavelets and Multiscale Representation

109

In the definition, h = oo causes no problem in the integral of (3.27) as long as a and q are both positive. The sensitive end is the ultraviolet limit (borrowing Mumford and Gidas's insightful terminology [225]) as the "wavelength" h —>• 0. Suppose for a given image u e B"(LP), its /?-modulus of continuity has a strict power law decay at the ultraviolet end: Then the decay has to be faster than a: ft > a; otherwise, the ultraviolet end will cause blowup in the integral (3.27). It is out of this consideration that a. is called the regularity index. It is easy to see that /^(L00) is precisely the Holder space Ha, which is usually defined by the single formula

Thus Besov spaces naturally generalize Holder spaces. If the scale space h e R+ is endowed with the measure dfjLa^(h) = h~^~aqdh, then the Besov condition simply becomes cop(u, •) e L 9 (E + , d/v

Our partners will collect data and use cookies for ad personalization and measurement. Learn how we and our ad partner Google, collect and use data. Agree & close