Calibrating image roughness by

estimating Lipschitz exponents, with

applications to image restoration

Alfred S. Carasso

Mathematical and Computational Sciences Division

Andr ́as E. Vlad ́ar

Precision Engineering Division

National Institute of Standards and Technology,

Gaithersburg, MD 20899

Phone 301-975-2705. FAX 301-975-3553.

Email alfred.carasso@nist.gov, andras.vladar@nist.gov

1

NISTIR 7438

August, 2007
Abstract

Most images f(x, y) are not smoothly differentiable functions of x and

y, but display edges, localized singularities, and other significant fine scale

roughness or texture. Correct characterization and calibration of image roughness is vital in many image processing tasks. The L

1 Lipschitz exponent α,

where 0 < α ≤ 1, measures fine scale image roughness provided the image is

relatively noise free. This paper describes a recently developed mathematical

technique for estimating α. The method is based on successively blurring the

image by convolution with increasingly narrower Gaussians, using commonly

available Fast Fourier Transform algorithms.

Instructive examples are used to illustrate the quantitative changes in

α that occur when an image is either degraded or restored. Of particular

interest are the documented changes in α that accompany APEX blind deconvolution of real images from the Hubble space telescope, from MRI and

PET brain scans, and from state of the art nanoscale Scanning Electron Microscopy. (APEX is the actual name of a blind deconvolution procedure, and

is not an acronym or abbreviation).

Additional applications include monitoring of image sharpness and imaging performance in imaging systems, evaluation of image reconstruction software quality, detection of abnormal fine structure in biomedical images, and

monitoring of surface finish in industrial applications.

Subject terms image roughness; image metrology; Lipschitz exponents;

image restoration; APEX deconvolution; Hubble telescope; MRI and PET

brain scans; nanoscale microscopy.

2
1 Introduction

This paper discusses an effective mathematical framework for quantifying the

unsmoothness of images, and then applies this methodology to some significant questions in image restoration.∗ An elegant computational technique,

based on Fast Fourier Transform (FFT) algorithms, translates this theory

into an important and practical image metrology tool. Interesting real images from such fields as astronomy, electron microscopy, and brain research

provide valuable illustrative examples. In all figures, the axis label ’log u’

refers to the natural logarithm of u.

Most natural images f(x, y) display edges, localized sharp features, and

other significant fine scale details or texture, and cannot be modeled as

smoothly differentiable functions of x and y. In many digital image processing tasks, it is necessary to provide prior information that specifies the

degree of unsmoothness in the unknown desired true image. If an image is

incorrectly postulated to be too smooth, the processing algorithm may produce a smoothed out version of the true image, in which critical diagnostic

information has been lost. The class BV (R2) of functions of bounded variation includes functions with sharp edges, as well as smooth continuously

differentiable functions. That class has been used extensively to model images over the last fifteen years or so.1−4 However, it has become increasingly

evident that the class BV (R2) does not allow for the type of fine scale sharp

structures commonly found in natural images, in addition to edges. For this

reason, image deblurring based on the BV (R2) assumption notoriously produces imagery in which important fine scale texture has been severely eroded.

This is known as the staircase effect.

3−7

It is now recognized that significantly

wider function classes, such as the so-called Lipschitz classes Λ(α, p, ∞), are

necessary to adequately describe natural images.

The L

1 Lipschitz exponent α, where 0 < α ≤ 1, is a mathematical index

that can capture the fine-scale content and degree of unsmoothness in an image, provided that image is relatively noise free. Images that are of bounded

variation (including smoothly differentiable images) have α = 1. The value

of α decreases systematically with increasing roughness. Images with significant non differentiable small scale structures typically have α ≪ 1. One may

∗The mathematical and computational techniques described in this paper are covered

by U.S. patent application number 10/928,308. Licensing inquiries should be directed to

Terry Lynch, jtlynch@nist.gov

3
also consider L

p

spaces with p > 1.

In a recent paper,7 a new computational approach for estimating image

Lipschitz exponents was developed, based on blurring the image with successively narrower Lorentzian (Cauchy) probability density functions. This

theory was then used in Ref. 7 to construct a new image deblurring method,

the so-called Poisson Singular Integral (PSI) method. That method is based

on prior knowledge of the L

2 Lipschitz exponent in the unknown true image.

An educated guess for that value, based on knowledge of exponent values

for similar images, is often sufficient to obtain useful sharpening. The major

result in Ref. 7 is the demonstration that the PSI method can recover texture

in cases where the BV (R2) approach fails completely.

The present paper is based on the use of technically more advantageous

Gaussian densities rather than Lorentzian densities. It complements the results in Ref. 7 by exploring further applications of Lipschitz exponents, in

particular the quantitative changes in α that occur when an image is either

degraded or restored. Of significant interest are the documented changes

in α that accompany APEX blind deconvolution of real images.6 Examples

include images from the Hubble space telescope, from magnetic resonance

imaging (MRI) brain scans, from functional cerebral positron emission tomography (PET) imaging, and from state of the art nanoscale scanning

electron microscopy (SEM). We also consider synthetically degraded images

and compare the changes in Lipschitz exponent that result from competing

restoration algorithms.

The above examples indirectly bear upon several other contexts where a

simple practical method of quantifying image roughness can be useful. Many

imaging systems suffer performance degradation over time and require periodic maintenance. In scanning electron microscopes, the shape of the electron beam changes with time, often without the user’s knowledge. Periodic

performance testing can be accomplished by monitoring sharpness degradation in the micrograph of a specially designed test object.8 Such degradation

can be quantified by measuring the increase in Lipschitz exponent. In some

imaging systems, a decrease in α might indicate an increase in system noise.

Measuring the α of a test image can also be used to evaluate and compare

the performance of competing image reconstruction packages. Automatic

measurement of α might be useful in monitoring the smoothness of surface

finishes in certain industrial applications. In some medical applications, measurement of α, along with other appropriate processing, might be helpful in

4
prescreening for certain types of abnormalities.

Sections 2-4 develop the mathematical preliminaries, culminating in Eq.

(15) which can be evaluated using FFT algorithms. Sections 5-7 illustrate

the use of that technique in estimating Lipschitz exponents. Section 8 studies

the behavior of Lipschitz exponents as images are synthetically degraded and

restored. Section 9 demonstrates the value of that tool in confirming APEX

blind sharpening of real blurred imagery. Finally, Section 10 contains the

conclusion.

2 Lipschitz classes and image roughness

A function f(x, y) ∈ BV (R2) may not be differentiable at a given set of

points, but it is constrained to satisfy

Z

R2

|f(x + h1, y + h2) − f(x, y)|dxdy ≤ Const |h|, |h| → 0, (1)

where

|h| = (h

2

1 + h

2

2

)

1/2

. (2)

However, from the standpoint of modeling texture, it is advantageous to

consider functions f(x, y) satisfying a weaker constraint, such as

Z

R2

|f(x + h1, y + h2) − f(x, y)|dxdy ≤ Const |h|

α

, |h| → 0, (3)

where α is fixed, and 0 < α < 1. Such an f is said to be of Lipschitz class

Λ(α, 1, ∞). A function f(x, y) satisfying Eq. (1) necessarily satisfies Eq. (3),

but not vice versa in general. Thus BV (R2) ⊂ Λ(α, 1, ∞).

One may also consider Lipschitz classes Λ(α, p, ∞), 1 ≤ p < ∞, consisting of functions f(x, y) satisfying

Z

R2

|f(x + h1, y + h2) − f(x, y)|

p

dxdy1/p

≤ Const |h|

α

, |h| → 0. (4)

The case p = 2 is important and was used in Ref. 7.

An effective method of estimating Lipschitz exponents can be based on

blurring the image by successive convolutions with increasingly narrower

Gaussians and evaluating the discrete L

p norm of the difference between

5
the blurred and original images. This L

p norm tends to zero as the Gaussian

approaches the Dirac δ-function. An important mathematical theorem relates the L

p Lipschitz exponent of the image to the rate at which the above

L

p norm tends to zero. Moreover, that theorem remains valid for the periodized image problem, in which case the convolution can be accomplished

quite easily using commonly available FFT algorithms.

3 The spaces Λ(α, p, ∞) and the Gaussian kernel

Define the Fourier transform hˆ(ξ, η) of h(x, y) ∈ L

1

(R2) by

F {h} = hˆ(ξ, η) ≡

Z

R2

h(x, y) exp{−2πi(ξx + ηy)}dxdy. (5)

For each fixed t > 0, consider the Gaussian kernel in R2

γ(x, y, t) = exp{−(x

2 + y2

)/4t}

4πt

, (x, y) ∈ R

2

. (6)

We have

γˆ(ξ, η, t) = exp{−t(ξ

2 + η2

)}. (7)

For each t > 0, define the linear operator Gt on L

p

(R2), 1 ≤ p < ∞, by

G

t

f =

Z

R2

γ(u, v, t)f(x − u, y − v)dudv. (8)

Since γ(x, y, t) approaches the Dirac δ-function as t tends to zero, it is not

surprising that k Gtf − f kp→ 0 as t tends to zero. However, the rate at

which this happens depends on the smoothness (or lack thereof) of f(x, y),

and this rate can be used to characterize f(x, y). There is a large literature

on that subject.9,10 We have from Theorem 4 in Ref. 9,

Theorem 1 Let Gt, t > 0, be the Gauss integral operator in (8), and let

0 < α < 1, 1 ≤ p < ∞. Then, f ∈ Λ(α, p, ∞) if and only if there is a

positive constant Kp such that

k G

t

f − f kp ≤ Kp t

α/2

, 0 < t ≤ 1. (9)

6
The above result can be used to fashion an image analysis tool. Theoretically,

given any image f(x, y) in L

1

(R2), one could use the Fourier transform in

Eq. (5) to form

F

n

G

t

f

o

= exp{−t(ξ

2 + η2

)}

ˆf(ξ, η), (10)

for sequences of positive t-values tending to zero. Inverse transformation is

always possible on account of the factor exp{−t(ξ

2 + η2

)}, and this can be

used to produce an infinite sequence of positive numbers

μn =

n

k G

tn f − f k1 / k f k1

o

, tn ↓ 0. (11)

If every such sequence (tn, μn), ultimately lies below the curve μ(t) =

C1 t

α/2

, with 0 < t ≤ t, and suitably chosen constants C1 > 0 and 0 < α < 1,

then

k G

t

f − f k1 ≤ C1 k f k1 t

α/2

, t ↓ 0, (12)

and f(x, y) ∈ Λ(α, 1, ∞) by Theorem 1. However, this requires handling

infinite domain Fourier integrals and does not lead to a practical procedure.

4 Periodized problems and FFT algorithms

A practical procedure can be realized by using the fact that Theorem 1

remains valid in the periodic case.9 We now consider the periodized image

problem, 10−15 and obtain the Fourier series analog of Eq. (10). Let Ω denote

the unit square −1/2 < x, y ≤ 1/2 in R2. The image f(x, y) is now viewed

as originally defined on Ω from which it is extended by periodicity to all of

R2. We redefine the image Fourier transform ˆf(ξ, η) by

ˆf(ξ, η) = Z

Ω

f(x, y) exp{−2πi(ξx + ηy)}dxdy. (13)

where ξ, η are now integers running from −∞ to +∞. The Poisson Summation Formula,10,16 can now be used to construct the periodized Gauss

operator Gt

∗

. That operator is defined by specifying its action on any given

image f(x, y). As shown in Ref. 7, Gt

∗

f is given by the complex Fourier

series

G

t

∗

f =

X∞

ξ,η=−∞

exp{−t(ξ

2 + η2

)}

ˆf(ξ, η) exp{2πi(xξ + yη)}. (14)

7
where ξ, η are integers running from −∞ to +∞. The factor exp{−t(ξ

2+η2

)}

assures uniform convergence of the Fourier series in Eq. (14).

Theorem 1 remains valid with Gt

∗

replacing Gt, and using Gt

∗f in Eq. (14)

in lieu of Gtf in Eq. (8). Next, consider the partial sum

SN (x, y) = X

N

ξ,η=−N

exp{−t(ξ

2 + η2

)}

ˆf(ξ, η) exp{2πi(xξ + yη)}. (15)

Because of uniform convergence, k Gt

∗f − SN kp can be made arbitrarily

small by choosing N large enough in Eq. (15), and Gt

∗

f in Eq. (14) can

be approximated by the finite sum SN . Given the 2J × 2J digitized image

f(x, y), the discrete Fourier transform17 is now the appropriate numerical tool

for approximating the finite Fourier series SJ . One can use FFT algorithms

to approximate the Fourier coefficients ˆf(ξ, η), − J ≤ ξ, η ≤ J, and then

apply the filter exp{−t(ξ

2 + η2

)}. An inverse FFT then yields an accurate

approximation to Gt

∗

f at each of the 2J ×2J pixels, for each small t > 0. We

may then examine the discrete L

p

relative error in Gaussian approximation

as t ↓ 0, and locate constants Cp and α such that

k G

t

∗

f − f kp≤ Cp k f kp t

α/2

, 0 < t ≤ t. (16)

In summary, the results of this section lead to an accurate numerical procedure, based on correct mathematical analysis, for assessing membership in

any Λ(α, p, ∞) space.

Remark 1: Pitfall at very small t > 0. We deal with discretely defined

high resolution 8-bit images f(x, y), typically of size 512×512 or 1024×1024

pixels. Such an f(x, y) may be viewed as a piecewise constant or trigonometric polynomial approximation to the original continuously defined image

intensity field f

∞(x, y), or as some other kind of finite dimensional representation of the hypothetical infinite dimensional object f

∞. The Lipschitz

exponent is predicated on a continuously defined image, and is, in fact, a

property of f

∞(x, y). On the other hand, all norms are equivalent on a finite

dimensional space. Hence, even if f

∞(x, y) is highly non smooth and not of

bounded variation, the discrete total variation norm for f(x, y) is always finite, though it may be very large. To estimate the non smoothness properties

of f

∞(x, y) by examination of the finite dimensional representation f(x, y)

will require some sagacity. As explained theoretically in Ref. 7, there is a

8
finite dimensionality pitfall in the above Gauss singular integral methodology that requires the exclusion of very small values of t > 0. As will be

amply demonstrated in the examples below, the behavior of k Gt

∗f − f kp

at very small t is a spurious artifact that must be ignored. This behavior

is disconnected from the true smoothness properties in the image intensity

field f

∞(x, y), and it gives a false reading for the Lipschitz exponent α.

5 Application to real images; the characteristic elbow

Our first example, in Figure 1, is the 512 × 512 Giza Pyramids image. The

FFT procedure discussed in Section 4 was used to obtain the L

1

relative

errors in Gauss approximation

μ(t) =k G

t

∗

f − f k1 / k f k1, (17)

at 400 values of t given by tn = 0.5(0.95)n, n = 1, 400, where Gt

∗

f is given by

Eq.(14). A plot of μ(t) versus t on a log-log scale produced the solid curve A

in Figure 1. Least squares fitting was used to find the two distinct majorizing

dashed straight lines Γ and Σ. For each dashed line, the y-axis intercept value

obtained by least squares was slightly increased so as to make each line lie visibly above the solid curve A; however, the slope of each line remains the same

as that obtained from least squares. The curve A exhibits a characteristic elbow shape. It consists of a straight line segment with slope ≈ 1 beginning near

log t = −20 and continuing to near log t = −10. There is then a transition

to a different regime, one that is more slowly increasing and that continues

to near log t = 0. The line Γ, defined by log μ(t) = 7.72 + 0.988 log t,

accurately captures the straight line trend in Eq. (17) for very small values

of t, while being grossly inaccurate at larger values of t. It was obtained by

excluding data corresponding to log t > −12 from the least squares fit. The

line Γ implies that k Gt

∗

f − f k1 < 2253 k f k1 t

0.988 for all t > 0. As

stressed in Remark 1, this correct statement primarily reflects the fact that

the 512 × 512 Giza Pyramids image lies in a finite dimensional space,7 but

does not reflect the smoothness properties of the intensity field f

∞(x, y) that

gave rise to the digitized Pyramids image. The slope of Γ gives a false reading

for α. Such Γ line behavior is a feature of every example in this paper. The

9
majorizing dashed straight line Σ, defined by log μ(t) = −1.31 + 0.159 log t,

accurately reflects behavior of Eq. (17) for −10 ≤ log t ≤ −1, while being grossly inaccurate at very small values of t. The line Σ was obtained

by excluding all data corresponding to log t < −10 from the least squares

fit. Note that this still leaves over 160 data points remaining. The behavior

along Σ indicates that k Gt

∗

f −f k1≤ 0.27 k f k1 t

0.159

, 0 < t ≤ 0.1, and this

is taken to be the true behavior in the Pyramids image. From Eq. (9), this

implies α = 2.0 × 0.159 = 0.318. Thus the Pyramids image lies in the space

Λ(0.318, 1, ∞), and is not of bounded variation, since this requires α = 1.

Estimates of α in any other discrete L

p norm can be obtained similarly. It

is recommended that data for very small values of t always be included in

log-log plots of μ(t), so as to enable clear identification of the spurious linear

trend, prior to rejecting that part of the data. This procedure of identifying

the Σ line by using least squares fitting on the ’slowly varying’ interval near

log t = 0, will be used throughout this paper. Blurring with Gaussians leads

to better defined elbows than does blurring with Lorentzians.

Remark 2: Parallel computation. The above blurring process with successively narrower Gaussians can obviously be implemented in parallel. For

large images, such parallel computation would be significantly more efficient.

6 Lipschitz exponents of some typical images

An interesting collection of 512 × 512 8-bit images is shown in Figure 2.

These images include natural as well as man-made objects, extending from

the nanoscale to the planetary scale. As in Figure 1, least squares fitting was

used in the slowly increasing region of the graph of log μ(t) vs log t, to find

the majorizing Σ line in each case. This leads to positive constants C1 and

α for each image such that

k G

t

∗

f − f k1≤ C1 k f k1 t

α/2

, 0 < t ≤ 0.1, (18)

where Gt

∗

is the periodized Gauss operator in Eq.(14). As seen in Table 1,

these images have L

1 Lipschitz exponents in the range 0.18 < α < 0.54.

10
7 Low Lipschitz exponents

Clusters of galaxies are examples of objects whose images can have very low

Lipschitz exponents. In Figure 3, two 512 × 512 images of the Abell 426

(Perseus) cluster are shown side by side. These images were obtained from

the NASA Space Telescope Science Institute Digitized Sky Survey. Each

image is centered on the exact same point in the sky. However, image (a)

covers a 1/2 Mpc square of the sky, while image (b) covers a 1 Mpc square

and thus has a higher density of bright pixels. The solid curves A and B

are the corresponding plots of μ(t) =k Gt

∗f − f k1 / k f k1 versus t, on a

log-log scale. Using least squares fitting on −6 ≤ log t ≤ 0, we find that

image (a) has an L

1 Lipschitz exponent α = 0.116, while image (b) has

α = 0.042. The corresponding Σ lines were not plotted in Figure 3 to avoid

clutter. Notice that the solid curves A and B have identical slopes of near

unity for log t < −10, even though images (a) and (b) have sharply distinct

smoothness properties. This confirms the observation in Remark 1 that the

slope of the Γ line is a finite dimensionality artifact, and is not connected to

image smoothness.

8 Degrading and restoring images; the Lipschitz exponent as an image metrology tool

This section deals with synthetically degraded images, and studies the behavior of Lipschitz exponents as images are degraded and restored. All Lipschitz

exponents in Figures 4 and 5 were estimated using least squares fitting of

the L

1

traces on the interval −7 ≤ log t < 0. As in Figure 3, corresponding

Σ lines were not plotted to avoid clutter. Notice that all traces in Figures 4

and 5 have identical slopes for log t ≤ −10, irrespective of the parent image.

Once again, this is spurious Γ line behavior that must be ignored.

Our first example, in Figure 4, involves noising and denoising the 512×512

Marilyn Monroe image. The original image in Fig. 4(a) has trace A and an

L

1

exponent α = 0.591. Salt and pepper noise with density 0.1 was added

to Fig. 4(a). This produced Fig. 4(b) with trace B. There is a very sharp

difference between traces A and B, and the Lipschitz exponent in Fig. 4(b)

has a substantially smaller value, α = 0.302. Obviously, the addition of noise

can significantly lower Lipschitz exponents, and this leads to two important

11
observations. First, as mentioned in the introduction, the Lipschitz exponent

reflects true image smoothness only when the image is relatively noise free.

Secondly, the sensitivity of α to the presence of noise may be used to monitor

and detect noise in imaging systems, by periodically reevaluating α for a

suitable test image.

We now consider two distinct denoising methods. The first method used

was the ’total variation’ method.1 This method assumes that the original

Marilyn Monroe image is of class BV (R2), and it is based on nonlinear

partial differential equations. The Marquina-Osher scheme18 was used with

∆t = 0.1(∆x)

2

, β = 0.0001, and regularization parameter λ = 2.0. Forward

integration was pursued for 300 time steps ∆t. This resulted in Fig. 4(c)

and the L

1 Gaussian trace C. Evidently, the total variation method has

eliminated a good deal of texture along with the noise in Fig. 4(c), and

trace C lies well below the original trace A. Fig. 5(c) has an L

1 Lipschitz

exponent α = 0.714, larger than the original value of 0.591. The second

denoising method used was 2D median filtering with a 3 × 3 neighborhood.

This produced Fig. 4(d) and L

1

trace D. Fig. 4(d) is a good approximation

to the original, and trace D is closer to trace A than is trace C. The L

1

Lipschitz exponent in Fig. 4(d) is α = 0.645.

Our second example, summarized in Figure 5, involves blurring and deblurring the Marilyn Monroe image. Fig. 5(a) results from synthetic blurring of Fig. 4(a) by convolution with a Lorentzian density. The solid curve

A in Figure 5 is the L

1 Gaussian trace for that blurred image, and it has

α = 0.887. Blurring without adding noise increases α, as blurring is generally

a smoothing operation.

Three mathematically distinct methods of deblurring Fig. 5(a) are examined under the present ideal conditions of perfect knowledge of the point

spread function (psf) and no added noise. Such a study highlights the intrinsic reconstructive ability of each scheme, and the accompanying changes

in α are of interest. The first method used was again the Marquina-Osher

total variation partial differential equation procedure.18 Here, we used ∆t =

0.1(∆x)

2

, β = 0.0001, and regularization parameter λ = 1000. Forward integration was pursued for 100 time steps ∆t. This resulted in Fig. 5(b) and

the L

1 Gaussian trace B. The second method used was the Lucy-Richardson

iterative procedure.19 That procedure was terminated after 200 iterations,

resulting in Fig. 5(c) and trace C. The third method used was the ”slow evolution from the continuation boundary” (SECB) procedure,19 with s = 0.01

12
and K = 1000. This produced Fig. 5(d) and trace D.

In Fig. 5(b), the BV (R2) assumption underlying the Marquina-Osher

method has resulted in a smoothed out reconstruction of the original. This is

particularly evident in the poorly resolved hair. As a result, the L

1 Lipschitz

exponent from trace B has the value α = 0.695, higher than the original

value of 0.591. In Fig. 5(c), the Lucy-Richardson method is a Bayesian

approach based on maximum likelihood. That method does not presume

a-priori smoothness in the unknown original image. On the other hand,

the Lucy-Richardson method typically requires large numbers of iterations

(> 2000) to reconstruct fine detail.19 Terminating the procedure after 200

iterations, as is quite commonly done, and was done in Fig. 5(c), produces

a smoothed out result. Nevertheless, Fig. 5(c) appears more life-like than

Fig. 5(b), and the hair is better resolved. Trace C almost coincides with

trace B for log t > −6, but falls below trace B for smaller t. Surprisingly,

the Lipschitz exponent in Fig. 5(c), with a value α = 0.714, is slightly

higher than in Fig. 5(b). The SECB procedure in Fig. 5(d) is a noniterative

direct deblurring method, based on solving an ill-posed fractional diffusion

equation backwards in time. The ”slow evolution” constraint

19 that is used

to regularize that problem does not require smoothness, and allows for non

smooth solutions not in BV (R2). For this reason, Fig. 5(d) appears to be

a higher quality approximation to the original image than the previous two

reconstructions. This is reflected in trace D which lies above traces B and

C and has a well-defined elbow. The L

1 Lipschitz exponent in Fig. 5(d) has

the value α = 0.616, quite close to the original value of 0.591.

9 Real imagery, APEX blind deconvolution,

and Lipschitz exponents

The above experiments on synthetically degraded imagery with perfectly

known psfs are instructive, and they indicate the Lipschitz exponent to be a

useful image metrology tool. However, in many imaging situations that result

in blurred imagery, the system psf is generally only poorly known. Effective

methods of ”blind deconvolution” that do not require knowledge of the psf are

of considerable interest, but it must be emphasized that such procedures are

fraught with serious mathematical and computational difficulties regarding

13
uniqueness and convergence.

The APEX method is a recently developed blind deconvolution technique

that is targeted at a specific class of shift-invariant blurs, in the form of 2D

isotropic, bell-shaped, heavy-tailed, probability density functions. Not all

images can be usefully enhanced with the APEX method, but the method

has been found useful in a variety of applications, including astronomy,

nanoscale electron microscopy, and MRI and PET brain imaging.6,20,21 As

a rule, APEX-detected psfs that successfully sharpen these images turn out

to be low exponent L ́evy stable laws.6,20,21 Such blurs are very far from Gaussians or Lorentzians, and are not generally known in the imaging literature.

A-priori knowledge about the solution is a fundamental aspect of solving

inverse problems in applied mathematics. This is especially the case with

blind deconvolution where severe ill-conditioning is compounded with non

uniqueness of solutions. The plausibility of APEX reconstructions must be

gauged by experienced analysts using independent considerations. Successful

APEX deconvolution should produce an image that is visibly sharper than

the original, yet one where the newly enhanced features can be traced back

to the original. Such sharpening should be accompanied by a measurable

drop in Lipschitz exponent. A decrease in Lipschitz exponent reflects the

extent to which small scale information has been reconstructed. Quantitative confirmation of sharpening is one of several elements that bear on the

plausibility of the APEX image. We now give several examples of Lipschitz

behavior before and after APEX processing. These examples involve real

images originating from significant areas of application.

Our first example, in Fig. 6, involves an iconic Hubble Space Telescope Whirlpool galaxy (M51) image. The original full resolution image was

taken with NASA’s ultra sophisticated Advanced Camera for Surveys (ACS).

APEX processing of the stepped down 710×1024 color image is discussed in

detail in Ref. 21. Here, we examine Lipschitz behavior in the blue component

of that image, before and after APEX processing. Using least squares fitting

on the interval −11 ≤ log t < 0, we find that the original image 6(a) has the

L

1 Gaussian trace A with Lipschitz exponent α = 0.416. The APEX image

6(b) again shows substantial additional detail, and it has the L

1 Gaussian

trace B with Lipschitz exponent α = 0.191. This is a remarkable 54% drop

in α.

We next consider brain imaging. APEX processing of the 512 × 512

sagittal MRI brain image in Fig. 7 was previously discussed in Ref. 6.

14
Using least squares fitting on the interval −7 ≤ log t < 0, we find that the

original MRI image 7(a) has the L

1 Gaussian trace A with Lipschitz exponent

α = 0.447. The APEX image 7(b) is significantly sharper visually. It has the

L

1 Gaussian trace B with Lipschitz exponent α = 0.296. This represents a

34% decrease in α.

In functional PET imaging, a positron emitting radionuclide is injected

into the patient, and used to tag glucose molecules in their course through the

brain. Performing specific mental tasks activates various parts of the brain,

causing increased glucose uptake and increased positron emission. Centers

of activity translate into relatively bright spots in the PET image. However,

blurring tends to attenuate such relative differences, resulting in a loss of

contrast. APEX processing of the transverse PET image in Fig. 8 was

previously discussed in Ref. 6. The 512×512 image in Fig. 8(a) is a stepped

up version of an original 128 × 128 image. Using least squares fitting on the

interval −6 ≤ log t < 0, we find that image 8(a) has the L

1 Gaussian trace A

with Lipschitz exponent α = 0.502. The APEX-processed image 8(b) shows

identical features as in 8(a), but contrast has been substantially increased.

Significantly, several bright spots appear in 8(b) that were not apparent

in 8(a). Image 8(b) has the L

1 Gaussian trace B with Lipschitz exponent

α = 0.394. This is a 21% decrease in α. The low resolution initial 128 × 128

PET format is a serious limitation in this example. Higher resolution PET

imagery would enable more substantial APEX sharpening.

The last example involves state of the art nanoscale SEM. Previous work

on applying APEX blind deconvolution to SEM imagery is discussed in Ref.

20. Since that time, significantly more powerful equipment has become

available.†It is an interesting fact that the high quality imagery produced

by this new instrumentation can still be usefully enhanced by the APEX

method. Moreover, Lipschitz analysis can be used to confirm such sharpening.

Fig. 9 deals with a Au/Pd decorated magnetic tape sample, taken with a

Hitachi S-4800 SEM using an accelerating voltage of 1 kV. The field of view is

1 μm. The original 1024×768 micrograph in 9(a) has the L

1 Gaussian trace A

with Lipschitz exponent α = 0.303. The APEX-processed micrograph 9(b),

†Certain commercial equipment is identified in this paper to adequately describe experimental procedures. Such identification does not imply recommendation or endorsement by

the National Institute of Standards and Technology, nor does it imply that the equipment

so identified is necessarily the best available for the purpose.

15
which is noticeably sharper, has the L

1 Gaussian trace B with Lipschitz

exponent α = 0.234. This is a 23% decrease in α. Least squares fitting on

the interval −7 ≤ log t < 0 was used.

Sharpening in Figs. 9(b) was confirmed independently, using a proprietary FEI software package.

Remark 3: Uncertainty in Lipschitz exponents. Estimated Lipschitz

exponent values in all examples in this paper, depend on the choice that was

made for the ’slowly varying’ interval in the trace of μ(t), as discussed in

Section 5. Least squares fitting on that interval produces the Σ line whose

slope is α/2. The particular choice of interval is a matter of judgment. A

somewhat larger or smaller interval might have been used in any given case,

resulting in a slightly different value for α. However, as was always the case

here, when Lipschitz exponents are used to evaluate the effectiveness of an

image processing algorithm on a given image, the same fitting interval should

be used in ’before and after’ comparisons.

10 Conclusion

A mathematical framework for quantifying the unsmoothness of images, leading to an effective computational tool for estimating image Lipschitz exponents, has been presented. The L

1 Lipschitz exponent α, where 0 < α ≤ 1, is

a useful image metrology tool that measures fine scale content, provided the

image is relatively noise free. It is found that natural images have L

1 values

of α lying between 0.2 and 0.5, typically. Such images are not of bounded

variation, since this requires α = 1.0.

The use of Lipschitz exponents in the quantitative evaluation of image

reconstruction procedures was studied. Using synthetically degraded images, it was shown that Lipschitz exponents can measure the extent to which

competing denoising algorithms remove texture along with the noise. Likewise, Lipschitz exponents can measure the ability of competing deblurring

algorithms to recover texture. Quantitative confirmation of fine structure

recovery is particularly important in blind image deconvolution, where the

cause of the blur is unknown. This was illustrated in the case APEX blind

deconvolution applied to Hubble space telescope imagery, scanning electron

microscopy, and MRI and PET brain scans.

16
References

[1] L. Rudin, S. Osher, and E. Fatemi, ”Nonlinear total variation based

noise removal algorithms,” Physica D. 60, 259–268 (1992).

[2] G. Aubert and P. Kornprobst, Mathematical Problems in Image Processing, Springer, New York, 2002.

[3] V. Caselles, J. M. Morel, G. Sapiro, and A. Tannenbaum Eds., Special

Issue on Partial Differential Equations and Geometry-Driven Diffusion

in Image Processing, IEEE Transactions on Image Processing 7, March

1998.

[4] J. Weickert, Anisotropic Diffusion in Image Processing, B. G. Teubner,

Stuttgart, Germany 1998.

[5] Y. Gousseau and J. M. Morel, ”Are natural images of bounded variation

?” SIAM J. Math. Anal. 33, 634–648 (2001).

[6] A. S. Carasso, ”The APEX method in image sharpening and the use

of low exponent L ́evy stable laws”, SIAM J. Appl. Math. 63, 593–618

(2002).

[7] A. S. Carasso, ”Singular integrals, image smoothness, and the recovery

of texture in image deblurring,” SIAM J. Appl. Math. 64, 1749–1774

(2004).

[8] A. E. Vlad ́ar, M. T. Postek, and M. P. Davidson, ”Image sharpness

measurement in scanning electron microscopy”, Scanning 20, pp. 24–34

(1998).

[9] M. H. Taibleson, ”On the theory of Lipschitz spaces of distributions on

Euclidean n-space. I. Principal properties”, Journal of Mathematics and

Mechanics 13, 407–478 ( 1964).

[10] P. L. Butzer and R. J. Nessel, Fourier Analysis and Approximation,

Academic Press, New York, 1971.

[11] J. W. Goodman, Introduction to Fourier Optics, Mc Graw-Hill, New

York, 1968.

17
[12] J. D. Gaskill, Linear Systems, Fourier Transforms, and Optics, Wiley,

New York, 1978.

[13] H. C. Andrews and B. R. Hunt, Digital Image Restoration, Prentice

Hall, Englewood Cliffs, NJ 1977.

[14] W. K. Pratt, Digital Image Processing, Second ed., Wiley, New York,

1991.

[15] R. C. Gonzalez and P. Wintz, Digital Image Processing, Second ed.,

Addison-Wesley, Reading, MA 1987.

[16] W. Feller, An Introduction to Probability Theory and its Applications,

Vol. 2, Second ed., Wiley, New York, 1971.

[17] E. O. Brigham, The Fast Fourier Transform and its Applications,

Prentice-Hall, Englewood Cliffs, New Jersey, 1992.

[18] A. Marquina and S. Osher, ”Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise

removal”, SIAM J. Sci. Comput. 22, 387–405 (2000).

[19] A. S. Carasso, ”Linear and nonlinear image deblurring: a documented

study”, SIAM J. Numer. Anal. 36, 1659–1689 (1999).

[20] A. S. Carasso, D. S. Bright, and A. E. Vlad ́ar, ”The APEX method and

real-time blind deconvolution of scanning electron microscope imagery”,

Optical Engineering 41, 2499-2514 (2002).

[21] A. S. Carasso, ”APEX blind deconvolution of color Hubble space telescope imagery and other astronomical data”, Optical Engineering 45,

107004 (2006).

18
Biography Alfred S. Carasso received the Ph.D degree in mathematics at

the University of Wisconsin in 1968. He was a professor of mathematics at

the University of New Mexico, and a visiting staff member at the Los

Alamos National Laboratory, prior to joining the National Institute of

Standards and Technology in 1982. His major research interests lie in the

theoretical and computational analyses of ill-posed problems in partial

differential equations, together with their application in image

reconstruction and computer vision. He is the author of original theoretical

papers, is a patentee in the field of image processing, and is an active

speaker at national and international conferences in applied mathematics.

19
Biography Dr. Andr ́as E. Vlad ́ar is the leader of the Scanning Electron

Microscope Metrology Project in the Nanoscale Metrology Group at the

National Institute of Standards and Technology. He holds M.S. (1977) and

Ph.D. (1984) degrees in Electrical Engineering from the Technical

University of Budapest, Hungary. Until 1991 he worked as a Research

Fellow at the Research Institute for Technical Physics of the Hungarian

Academy of Sciences, Budapest. Dr. Vlad ́ar is an expert in scanning

electron microscope (SEM) critical dimension (CD) metrology; he has

developed several metrology systems based on SEMs. Prior to joining NIST

in 1999, he was a member of the technical staff at the HewlettPackard ULSI

Research Laboratory and primarily worked on accurate, fast, and effective

dimensional measurement methods for 100 nm silicon integrated circuit

technology. At HewlettPackard he established sound metrology methods,

procedures and practices that significantly improved the speed and quality

of dimensional measurements, and surface inspection and film thickness

20
measurements in the Class 1 integrated circuit fabrication facility. He

identified the main sources of measurement errors and developed several

successful methods to eliminate many of them.

21
List of Figure Captions

Fig. 1. 512 × 512 Giza Pyramids image is not of bounded variation, but

has an L

1 Lipschitz exponent α ≈ 0.32. This follows from the graphical

use of Theorem 1, using the FFT techniques discussed in Section 4. Solid

curve A is a plot of μ(t) =k Gt

∗f − f k1 / k f k1 versus t, on a log-log

scale, where Gt

∗

f is defined in Eq.(14). Majorizing dashed straight line Γ,

defined by log μ(t) = 7.72 + 0.988 log t, accurately captures linear behavior

in Eq. (17) for very small values of t, but is grossly inaccurate at larger

values of t. Linear behavior at very small t is misleading, and is unrelated

to true image smoothness. (See Remark 1). Majorizing dashed straight line

Σ, defined by log μ(t) = −1.31 + 0.159 log t, accurately reflects behavior for

−10 ≤ log t ≤ −1, while being grossly inaccurate at very small t. Behavior

along Σ is taken to be true behavior in Giza Pyramids image, and indicates

k Gt

∗

f − f k1≤ 0.27 k f k1 t

0.159

, 0 < t ≤ 0.1. From Eq. (9), this implies

image ∈ Λ(0.318, 1, ∞).

Fig. 2. The above 512 × 512 8-bit images have L

1 Lipschitz exponents α

in the range 0.18 < α < 0.54, and are not of bounded variation. See Table 1.

Fig. 3. Galaxy clusters can have very low Lipschitz exponents. Above

images of Abell 426 (Perseus) cluster obtained from NASA STScI Digitized

Sky Survey. Images (a) and (b) are of size 512 × 512 and each is centered

on the exact same point in the sky. Image (a) covers a 1/2 Mpc square of

the sky, while image (b) covers a 1 Mpc square. Solid curves A and B are

the corresponding plots of μ(t) =k Gt

∗f − f k1 / k f k1 versus t, on log-log

scale. Using least squares fitting on −6 ≤ log t ≤ 0, we find that image

(a) has an L

1 Lipschitz exponent α = 0.116, while image (b) has α = 0.042.

Corresponding Σ lines were not plotted to avoid clutter. Solid curves A and

B have identical slopes ≈ 1 for log t < −10. This confirms the observation

in Remark 1 that the slope of the Γ line is a finite dimensionality artifact,

and is not connected to image smoothness.

Fig. 4. Image Denoising. Noise addition can artificially lower Lipschitz

exponent α, while some noise removal algorithms can eliminate texture and

22
increase α. Gaussian traces A, B, C, and D correspond to images (a), (b), (c),

and (d) respectively. Least squares fitting on (−7, 0) produced the following

results. (Σ lines not plotted to avoid clutter). (a) Original Marilyn Monroe

image has L

1

exponent α = 0.591. (b) Adding salt and pepper noise to image

(a) leads to α = 0.302, a very noticeable reduction. (c) Nonlinear partial differential equation ‘total variation’ noise removal algorithm eliminates noise,

as well as texture, and results in α = 0.714, higher than in original. (d) 2D

median filtering produces closer approximation to original, with α = 0.645.

Fig. 5. Image Deblurring. Not all deblurring algorithms are equally effective. Gaussian traces A, B, C, and D correspond to images (a), (b), (c),

and (d) respectively. Least squares fitting on (−7, 0) produced the following results. (Σ lines not plotted to avoid clutter). (a) Syntheticaly blurred

Marilyn Monroe image has L

1

exponent α = 0.887. (b) Nonlinear partial differential equation ‘total variation’ deblurring of (a) erodes texture and leads

to α = 0.695. (c) Lucy-Richardson deblurring has α = 0.714. (d) SECB

deblurring gives α = 0.616. This is closest to original image value α = 0.591.

Fig. 6. Evaluating the effectiveness of APEX blind deconvolution of Whirlpool

galaxy (M51) image. See Ref. 25. (a) Original Hubble Space Telescope

710×1024 image has Gaussian trace A and L

1 Lipschitz exponent α = 0.416.

(b) APEX processed image has Gaussian trace B and L

1 Lipschitz exponent

α = 0.191, a 54% decrease. Least squares fitting on (−11, 0) was used.

Fig. 7. Evaluating the effectiveness of APEX blind deconvolution of sagittal

brain MRI image. See Ref. 6. (a) Original 512 × 512 image has Gaussian

trace A and L

1 Lipschitz exponent α = 0.447. (b) APEX processed image

has Gaussian trace B and L

1 Lipschitz exponent α = 0.296, a 34% decrease.

Least squares fitting on (−7, 0) was used.

Fig. 8. Evaluating the effectiveness of APEX blind deconvolution of transverse functional PET brain slice. See Ref. 6. (a) Original 512 × 512 image

has Gaussian trace A and L

1 Lipschitz exponent α = 0.502. (b) APEX processed image has Gaussian trace B and L

1 Lipschitz exponent α = 0.394, a

21% decrease. Least squares fitting on (−6, 0) was used.

Fig. 9. Evaluating the effectiveness of APEX blind deconvolution in nanoscale

23
Scanning Electron Microscopy. (a) Original 1024×768 micrograph of Au/Pd

decorated magnetic tape sample has Gaussian trace A and L

1 Lipschitz exponent α = 0.303. (b) APEX processed image has Gaussian trace B and

L

1 Lipschitz exponent α = 0.234, a 23% decrease. Least squares fitting on

(−7, 0) was used.

24
TABLE 1.

Values of (C1, α) in Eq.(18) for each image f(x, y) in Figure 2.

Image Size L

1 norm values of C1 and α

Sydney Opera House 5122 C1 = 0.42, α = 0.476

USS Theodore Roosevelt 5122 C1 = 0.69, α = 0.535

Washington DC Landsat 5122 C1 = 0.80, α = 0.236

Heart ventricle SPECT 5122 C1 = 0.70, α = 0.330

Transverse brain PET 5122 C1 = 1.07, α = 0.502

Sagittal brain MRI 5122 C1 = 0.83, α = 0.447

Weather satellite Hurricane 5122 C1 = 0.82, α = 0.411

Nanoscale electron micrograph 5122 C1 = 0.69, α = 0.183

USAF resolution chart 5122 C1 = 2.08, α = 0.203

25
A

Figure 1: 512 × 512 Giza Pyramids image is not of bounded variation....
Figure 2: The above 512 × 512 8-bit images have L

1 Lipschitz exponents......
(a) (b)

Figure 3: Galaxy clusters can have very low Lipschitz exponents.....
(a) (b)

(c) (d)

Figure 4: Image Denoising. Noise addition can artificially ......
(a) (b)

(c) (d)

Figure 5: Image Deblurring. Not all deblurring algorithms ........
(a) (b)

Figure 6: Evaluating the effectiveness.....Whirlpool galaxy
(a) (b)

Figure 7: Evaluating the effectiveness......sagittal brain MRI ....
(a) (b)

Figure 8: Evaluating the effectiveness......transverse functional PET...
(a) (b)

Figure 9: Evaluating the effectiveness.....magnetic tape sample....