NISTIR 7632

Bochner Subordination, 

Logarithmic Diffusion Equations, 

and Blind Deconvolution of 

Hubble Space Telescope Imagery 

and Other Scientific Data

Alfred S. Carasso
NISTIR 7632

Bochner Subordination, 

Logarithmic Diffusion Equations, 

and Blind Deconvolution of 

Hubble Space Telescope Imagery 

and Other Scientific Data 

Alfred S. Carasso 

Mathematical and Computational Sciences Division 

Information Technology Laboratory 

September 2009

U.S. Department of Commerce 

Gary Locke, Secretary 

National Institute of Standards and Technology 

Patrick D. Gallagher, Director

 
BOCHNER SUBORDINATION, LOGARITHMIC DIFFUSION

EQUATIONS, AND BLIND DECONVOLUTION OF HUBBLE SPACE

TELESCOPE IMAGERY AND OTHER SCIENTIFIC DATA

ALFRED S. CARASSO∗

Abstract. Generalized Linnik processes and associated logarithmic diffusion equations can be

constructed by appropriate Bochner randomization of the time variable in Brownian motion and the

related heat conduction equation. Remarkably, over a large but finite frequency range, generalized

Linnik characteristic functions can exhibit almost Gaussian behavior near the origin, while behaving

like low exponent isotropic L ́evy stable laws away from the origin. Such behavior matches Fourier

domain behavior in a large class of real blurred images of considerable scientific interest, including

Hubble space telescope imagery and scanning electron micrographs. This paper develops a powerful

blind deconvolution procedure based on postulating system optical transfer functions (otf) in the

form of generalized Linnik characteristic functions. The system otf and ‘true’ sharp image are

then reconstructed by solving a related logarithmic diffusion equation backwards in time, using the

blurred image as data at time t = 1. The present methodology significantly improves upon previous

work based on system otfs in the form of L ́evy stable characteristic functions. Such improvement

is validated by the substantially smaller image Lipschitz exponents that ensue, confirming increased

fine structure recovery. These results resolve the unexplained appearance of exceptionally low L ́evy

stable exponents in previous work on the same class of images. The paper is illustrated with striking

enhancements of gray scale and colored images.

Key words. image deblurring; blind deconvolution; Bochner subordination; generalized Linnik

laws; low exponent L ́evy stable laws; fractional and logarithmic diffusion equations; Lipschitz exponents; Hubble space telescope; scanning electron microscope; Whirlpool galaxy; Starburst galaxy.

AMS subject classifications. 35R25, 35B60, 60E07, 68U10.

1. Introduction. This paper deals with blind deconvolution of a class of real

blurred images of considerable scientific interest, as opposed to synthetically degraded

phantoms. A preselected class of trial optical transfer functions (otfs), not previously

known in image analysis, is shown to produce striking enhancements of Hubble space

telescope images and other astronomical data, as well as useful sharpening of scanning

electron micrographs of interest in Nanotechnology. Such improvements may not be

apparent in the reduced size images in the printed issue of this journal. However,

significant enhancement becomes evident when the on-line version of this paper is

viewed at full size on a modern high resolution device, such as a wide screen, active

matrix, liquid crystal display (LCD) monitor. As noted below, this degree of fine

structure reconstruction is difficult to achieve with existing techniques.

Blind deconvolution seeks to deblur an image without knowing the cause of the

blur. The methodology used here is based on identifying a plausible blurring otf from

within a restricted class of candidate blurs. The choice of trial otfs is governed by the

Fourier domain behavior in the blurred imagery under consideration. Previous work

on the above class of problems [10], [11], was based on candidate otfs in the form

of isotropic L ́evy stable characteristic functions [18], [27], [31], and the use of timereversed diffusion equations involving fractional powers of the negative Laplacian. An

important observation in [7], [10], [11], is that the successful deblurring otfs are characterized by low L ́evy exponents, with typical values less than 0.5, while Gaussian otfs

have exponent 2. Such low values are quite exceptional in applications where L ́evy

stable laws appear. The physical origin, if any, of such exponent values is not known.

∗Computing and Mathematical Sciences Laboratory, National Institute of Standards and Technology, Gaithersburg, MD 20899. (alfred.carasso@nist.gov).

1
2 ALFRED S. CARASSO

The present method is based on otfs in the form of generalized Linnik characteristic

functions [16], [26], [32], and the use of time-reversed diffusion equations involving

the logarithm of the negative Laplacian plus the identity. We show that this results

in higher quality reconstructions than previously obtained. This improvement is validated by the substantially smaller image Lipschitz exponents that ensue, confirming

increased fine structure recovery [8], [13]. Moreover, the high frequency behavior in

the successful Linnik deblurring otfs resolves the unexplained appearance of low L ́evy

exponents in [10] and [11].

While Linnik otfs have distinctly different behavior near the origin than do L ́evy

stable otfs, there is a close connection between these two classes of objects. This connection centers on the notion of Bochner subordination of stochastic processes [5], [6],

[18], [32], whereby a given Markov process X(t), t > 0, is observed in new stochastic

‘operational time’ S(t), rather than in standard clock time t, resulting in the process

X(S(t)). This concept has found fruitful application in several branches of science,

engineering, and finance [1], [15], [19], [20], [35]. In that context, generalized Linnik

characteristic functions describe the family of processes obtained when isotropic L ́evy

stable motions are observed in the randomized operational time Γ(t), where for each

t > 0, Γ(t) obeys a Gamma distribution on u ≥ 0. The important special case

of Brownian motion observed in stochastic time Γ(t) is called the Variance Gamma

Process [28], [29]. This has been used successfully in option pricing.

Bochner subordination also plays a role in operator semigroup theory and evolutionary partial differential equations. By combining subordination with the HilleYosida theorem, a functional calculus for semigroup generators is developed in [30],

whereby entirely new semigroups can be created by randomizing the time variable. In

[12], it is shown that the subordinator Γ(t) discussed above is one of a large family of

stochastic time changes with the remarkable property of converting all C0 semigroups

on a given complex Banach space B, into holomorphic semigroups on B.

In this paper, the use of generalized Linnik densities as trial blurring kernels is

primarily motivated by their mathematical properties. The semigroup property of

infinite divisibility [18], [27], [32], the Fourier domain Gaussian-like behavior near the

origin, together with the monotone convex high frequency behavior which imitates

the blurred image data, all play essential roles. We do not claim an actual physical

basis for Linnik otfs, and the significance of this paper ultimately rests on the quality

of the reconstructions. However, our results strongly suggest a possible physical basis

for such Linnik otfs, mediated by subordination as in previously mentioned applications. While Gaussian otfs, synonymous with Brownian motion, are a standard first

choice in image reconstruction, image acquisition involves the interaction of several

interfacing optical and electronic devices, each producing a small distortion of the

input signal. Additional aberrations may result from the scattering properties of the

medium through which radiation propagates. A significant empirical discovery [21],

[22], [23], is that a large variety of electro-optical imaging devices have otfs given by

L ́evy stable characteristic functions. Very recently [3], an optical material was created

with a specific kind of inhomogeneity, in which the scattering of light waves results

in L ́evy flights rather than Brownian motion. A prescribed L ́evy exponent for this

scattering process can be engineered by proper synthesis of the inhomogeneity. As

is well-known, L ́evy stable motions are related to Brownian motion through Bochner

subordination. A second subordination, involving the Γ(t) process, would result in

generalized Linnik otfs. Conceivably, Linnik otfs may play an unsuspected role in

many imaging situations.
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 3

The blind deconvolution procedure presented here is fundamentally different from

such variational blind formulations as [14] and [24]. These methods aim to solve the

blind deconvolution problem in full generality, by minimizing an appropriate cost

functional. Such techniques do not appear useful for the present class of problems.

The method in [14] is primarily intended for ‘blocky’ images, and point spread functions with edges. Other promising approaches such as [2], exploit prior knowledge

of the edge map in the unknown image. Astronomical images are typically diffuse

and partly amorphous, and are very far from being of bounded variation. Use of [14]

tends to eliminate individual star clusters and other localized bright areas, which are

interpreted as unwanted noise. The method in [24] requires a prior guess for the system otf. However, even with a good guess, this method often returns a questionable

sharp image associated with a new, but physically impossible otf [9]. If and when

comprehensive blind deconvolution methods become available, the methodology presented here is likely to remain useful. It can provide independent reconstructions for

comparison, and it can provide valuable initial guesses that might be refined by more

elaborate procedures.

2. Imagery with monotone convex Fourier decay. All images in this paper

are 8-bit images, with pixel values scaled from 0 to 255. Given an image g(x, y), we

define its Fourier transform by

gˆ(ξ, η) ≡

Z

R2

g(x, y)e

−2πi(ξx+ηy)

(1) dxdy,

and we denote by k g kr the L

r

(R2) norm of g(x, y). Typically, r = 1 or 2. We also

use the ‘total variation’ seminorm

(2) k g kTV ≡k ∇g k1 .

In general, a real blurred image g(x, y) contains noise n(x, y) which may be multiplicative,

(3) g(x, y) = ge(x, y) + n(x, y), k n kr≪k ge kr, r = 1, 2,

where ge(x, y) denotes the exact blurred image that would have been recorded in the

absence of any errors or noise. Neither ge(x, y) nor n(x, y) are known, only their sum

g(x, y).

The example in Figure 1 is illustrative of a large class of images of galaxies and

other astronomical objects. An original 512×512 pixels Pleiades (M45) image g(x, y),

is shown on the left hand side of Figure 1. Of particular interest is the plot of

log |gˆ(ξ, 0)| vs |ξ|, shown as the dashed curve immediately below. While that trace

is locally highly oscilllatory, it is globally monotone decreasing and convex. The

solid curve p(ξ) captures the gross behavior in log |gˆ(ξ, 0)| away from the origin. The

deblurred image on the right hand side is a good representation of a sharp Pleiades

image f(x, y), and we shall consider f(x, y) to be the ‘true’ Pleiades image for the

purpose of this discussion. We note that log |

ˆf(ξ, 0)| exhibits a very similar globally

monotone convex trace, whose gross behavior is captured by the solid curve q(ξ).

The curve p decays faster than the curve q. Both p and q were chosen so as to have

a maximum value of 0, attained at ξ = 0. In addition, as shown in the bottom

drawing in Figure 1, the curves log |gˆ(ξ, 0)| and log |

ˆf(ξ, 0)| are such that away from

the origin, the ratio of these two quantities remains fairly constant, fluctuating around
4 ALFRED S. CARASSO

 Original M45 True M45 ?

Original Pleiades (M45) True Pleiades ?

p q

original

true

original/true

Behavior in Fourier transform domain

Fig. 1. Above Pleiades Fourier transform behavior is representative of a large class of images.

Away from the origin, blurred image data log |gˆ(ξ, 0)| on left, and ‘true’ image data log |

ˆf(ξ, 0)| on

right, exhibit similar highly oscillatory globally convex monotone behavior, but left hand data decays

faster. Moreover, the ratio log|gˆ(ξ, 0)|/ log |

ˆf(ξ, 0)| remains fairly constant, fluctuating around some

mean value larger than unity.
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 5

some mean value larger than 1. Accordingly,

log |gˆ(ξ, 0)|/ log |

ˆf(ξ, 0)| ≈ Constant > 1, |ξ| ≫ 1,

p(ξ)/q(ξ) ≈ Constant > 1, |ξ| ≫ 1.

(4)

A comparison of the left and right hand sides in Figure 1 invites the following,

possibly naive, interpretation. Assume that the original image g(x, y) on the left is a

blurred version of the true image f(x, y) on the right, and results from the convolution

of f(x, y) with an isotropic, shift invariant, point spread function h(x, y). Note that

h(x, y) is a 2D probability density function since it is non-negative and integrates to

unity. The Fourier transform of h(x, y) is the optical transfer function hˆ(ξ, η). That

object is necessarily a 2D characteristic function and must obey Bochner’s theorem

on positive definite functions [4], [17]. Since h(x, y) is isotropic, so is hˆ(ξ, η), and we

write

hˆ(ξ, η) = ˆh(ρ), ρ =

p

ξ

2 + η2 (5) .

Assume further that the characteristic function hˆ(ρ) is positive and monotone decreasing. This is the case for a large class of isotropic otfs, based on Schoenberg’s theorem

on completely monotone functions [17], [33]. With ⊗ denoting convolution in L

1

(R2),

we have

ge(x, y) + n(x, y) ≡ g(x, y) = h(x, y) ⊗ f(x, y),

log |gˆ(ξ, η)| = log hˆ(ρ) + log |

ˆf(ξ, η)|,

log |gˆ(ξ, 0)| ∼ p(ξ), log |

ˆf(ξ, 0)| ∼ q(ξ), |ξ| ≫ 1,

log hˆ(ρ) ∼ p(ρ) − q(ρ), ρ ≫ 1.

(6)

Thus, under the preceding assumptions, behavior in hˆ(ρ) away from the origin can be

inferred from that in p − q. In practice, given an original image g(x, y) with Fourier

domain behavior similar to that on the left in Figure 1, we may find p(ξ) but not q(ξ).

However, for a large class of blurred images g(x, y), the a-priori assumption that the

corresponding true image f(x, y) has Fourier behavior similar to that on the right of

Figure 1, and obeys (4), is found to lead to useful reconstructions. Accordingly, using

(4) and (6), we postulate that with some unknown positive constant c < 1

log hˆ (7) (ρ) ∼ (1 − c)p(ρ), ρ ≫ 1.

This paper seeks to identify the system otf by a least squares fit to the data log |gˆ(ξ, 0)|

with a suitable monotone decreasing function p(ξ) as above, with p(0) = 0. That

function must be such that for any constant b with 0 < b < 1, the function bp(ρ) is

the logarithm of a characteristic function. This is the case if and only if exp{p(ρ)}

is an infinitely divisible characteristic function [18], [27].

3. Infinitely divisible otfs and subordinated diffusion equations. In the

deconvolution problem h(x, y) ⊗ f(x, y) = g(x, y) to be analyzed below, the targeted

class of images directs attention to the following three types of infinitely divisible
6 ALFRED S. CARASSO

isotropic otf

hˆ(ρ) = exp(−σρ2), σ > 0, ρ =

p

ξ

2 + η2, Gaussian

hˆ(ρ) = exp(−σρα), σ > 0, 0 < α ≤ 2, Levy stable  ́

hˆ(ρ) = (1 + σρα)

−λ

, σ, λ > 0, 0 < α ≤ 2. Generalized Linnik

(8)

Each of these three types is associated with a diffusion equation. The Gaussian

case, representing Brownian motion, corresponds to the heat equation as is wellknown. The other two types involve stochastic processes and corresponding diffusion

equations that result from randomizing the time variable in Brownian motion. We

show this by introducing the notion of subordinated semigroup [12], [30], [32].

Consider a family of functions {pt(u)} indexed by t ≥ 0 and defined on u ≥ 0,

and such that for each fixed t > 0, pt(u) is a probability density function on u ≥ 0.

Assume that pt(u) ∗ ∗ps(u) = pt+s(u), where ∗ ∗ denotes convolution on u ≥ 0, and

that pt(u) → δ(u) as t ↓ 0. For fixed t > 0, the Laplace transform of pt(u) is defined

by

L(pt) = Z ∞

0

pt(u)e

−uz (9) du, Re z > 0.

Definition 1. A Bochner subordinator is a family {pt(u)} as defined above, such

that L(pt) = exp{−tψ(z)}, where ψ(z) is holomorphic for Re z > 0 and continuous

for Re z ≥ 0, with Re ψ(z) ≥ 0. Moreover, ψ(0) = 0, and ψ

′

(x) is completely

monotone on x > 0. The function ψ(z) is called the Bernstein function.

The following two subordinator examples, involving the Gamma and Inverse

Gaussian families respectively, yield closed form expressions.

pt(u) = {Γ(λt)}

−1σ−λtuλt−1

e

−u/σ, σ, λ > 0, L(pt) = (1 + σz)−λt

,

pt(u) = {te−t

2/4u}/

√

4πu3, L(pt) = e

−t

√

z

.

(10)

The Inverse Gaussian is the special case β = 1/2 of the important L ́evy subordinator

family of index β, where 0 < β ≤ 1. For β 6= 1/2, this subordinator is not known in

closed form and is defined implicitly by

pt(u) = L

−1

(e

−tzβ

(11) ), Re z > 0.

In the above, whenever multivalued Bernstein functions ψ(z) appear, the particular

branch of ψ(z) such that Re ψ(z) > 0 on Re z > 0, is always understood.

Let U(t) be the operator semigroup e

−tA associated with the well-posed linear

evolution equation wt = −Aw, t > 0, w(0) = f, in a Banach space B. Bochner

subordination involves randomizing the time variable in U(t) as follows. Let {pt(u)}

be a Bochner subordinator with Bernstein function ψ(z), and define a new semigroup

T (t) by

T (t)f =

Z ∞

0

(12) pt(s)U(s)fds, f ∈ B

Formally, T (t) is the semigroup e

−tψ(A) associated with the well-posed linear evolution

equation wt = −ψ(A)w, t > 0, w(0) = f. See [12], [30], [32].
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 7

Next, let A be a linear partial differential operator in the space variables in R2,

with constant coefficients. Let ˆh(ξ, η, t) denote the Fourier transform of the Green’s

function for the linear evolution equation wt = −Aw. The Green’s function hˆ

ψ(ξ, η, t)

for wt = −ψ(A)w can be obtained from hˆ(ξ, η, t), using the composition law for

subordinated characteristic functions [32, Chapter 6],

hˆ

ψ(ξ, η, t) = exp h

−tψ(− log{hˆ(ξ, η, 1)})

i

(13) .

We may now apply the foregoing analysis to discover the diffusion equations associated

with (8). We begin with the well-posed forward heat equation and its Gaussian

Green’s function

wt = a∆w, t > 0, a > 0,

hˆ(ξ, η, t) = exp(−tσρ2), σ = 4aπ2, ρ =

p

ξ

2 + η2.

(14)

L ́evy subordination with ψ(z) = z

β

in (14), leads to a well-posed forward fractional

diffusion equation, with a L ́evy stable Green’s function

wt = −



b(−∆)β



w, t > 0, b > 0, 0 < β ≤ 1,

hˆ(ξ, η, t) = exp(−tδρ2β), δ = b(4π

2

)

β

, ρ =

p

ξ

2 + η2.

(15)

Next, Gamma subordination with ψ(z) = λ log{1+σz} in (15), results in a well-posed

forward logarithmic diffusion equation, with a generalized Linnik Green’s function

wt = −



λ log{1 + c(−∆)β}



w, t > 0, c > 0, 0 < β ≤ 1,

hˆ(ξ, η, t) = (1 + γρ2β)

−λt, γ = c(4π2

)

β

, ρ =

p

ξ

2 + η2.

(16)

The generalized Linnik process may also be viewed as resulting from a single subordination of Brownian motion with the Bernstein function ψ(z) = λ log(1 + σzβ

).

Given the deconvolution problem h(x, y) ⊗ f(x, y) = g(x, y), with a known otf

hˆ(ξ, η) in the form of one the three types in (8), we view hˆ(ξ, η) as the Green’s

function hˆ(ξ, η, t) at time t = 1, in the corresponding forward evolution equation

wt = −Lw, in one of (14), (15), or (16). Deconvolution is mathematically equivalent

to solving wt = −Lw backwards in time, given the noisy blurred image g(x, y) as data

at time t = 1. For 0 < t < 1, w(x, y, t) is a partially deblurred image. The fully

deblurred image f(x, y) is the solution at time t = 0.

Such backwards continuation in diffusion equations is notoriously ill-posed. The

SECB method, (see [25] for an up to date discussion), is a well-regularized continuation

procedure that takes into account the presence of noise in g(x, y) at t = 1. With n(x, y)

as in (3), let constants ǫ, M, be given such that

(17) k w(0) k2≡k f k2≤ M, k w(1) − g k2≡k n k2≤ ǫ, ǫ ≪ M.

For any constant K > 0 such that K ≪ M/ǫ define s

∗ by

s

∗

(18) = log {M/(M − Kǫ)} / log(M/ǫ).

The slow evolution constraint applied to the backwards solution of wt = −Lw requires

that there exist a known small constant K > 0 and a known fixed s ≫ s

∗

, such that

(19) k w(s) − w(0) k2≤ Kǫ.
8 ALFRED S. CARASSO

Given the regularization parameters K and s, the SECB solution is that initial value

w

†

(x, y, 0) = f

†

(x, y) which minimizes

k w(1) − g k

2

2 +K−2

k w(s) − w(0) k

2

2

(20) ,

over all choices of initial values w(x, y, 0) in L

2

.

We may find the SECB deblurred image f

†

(x, y) in closed form in the Fourier

transform domain. We have, with z denoting the complex conjugate of z,

ˆf

†

(ξ, η) = hˆ(ξ, η, 1) ˆg(ξ, η)

|hˆ(ξ, η, 1)|

2 + K−2|1 − hˆs(ξ, η, 1)|2

(21) ,

leading to f

†

(x, y) upon inverse transforming. Of particular interest is the partially

deblurred image w

†

(x, y, t) defined by

wˆ

†

(ξ, η, t) = hˆ(ξ, η, t) hˆ(ξ, η, 1) ˆg(ξ, η)

|hˆ(ξ, η, 1)|

2 + K−2|1 − hˆs(ξ, η, 1)|2

(22) , 0 ≤ t < 1.

This can be efficiently implemented using FFT algorithms.

4. Blind deconvolution by marching backwards in time. As indicated in

Section 2, given an image g(x, y), we seek to identify the system otf by a least squares

fit to the data log |gˆ(ξ, 0)| with a suitable monotone decreasing function p(ξ) with

p(0) = 0, as in Figure 1. The function p(ρ) is chosen so that exp{p(ρ)} is one of the

three characteristic function types in (8). We then define

ˆh(ξ, η, t) = exp{tp(ρ)}, 0 ≤ t ≤ 1, ρ =

p

ξ

2 + η2 (23) .

where hˆ(ξ, η, t) is the Green’s function for the corresponding diffusion equation in

one of (14), (15), or (16). Using this in (22), together with g(x, y) as data at t = 1,

great benefit derives from the ability to perform the deconvolution in slow motion by

marching backwards in time in the diffusion equation wt = −Lw. Recall that from

(7), the system otf satisfies log hˆ(ρ) ∼ (1−c)p(ρ), ρ ≫ 1, for some unknown positive

constant c < 1. Continuation all the way to t = 0 in (22), with hˆ(ξ, η, t) as in (23), is

equivalent to falsely setting c = 0 in (7). This necessarily produces oversharpening.

Displaying the evolution of w

†

(x, y, t) as t decreases from 1 to 0, allows for monitoring

the deblurring process. At first, the partially deblurred images w

†

(x, y, t) become

sharper and slightly noisier as t decreases. However, as expected, ringing, noise,

and other artifacts begin to appear, indicating that continuation has proceeded too

far. Diagnostic statistical information about w

†

(x, y, t) can also be calculated for

selected values of t, as t decreases. Of particular interest are the discrete L

1 norm

k w

†

(., t) k1, and the discrete TV norm k ∇w

†

(., t) k1. In well-behaved deconvolution,

image flux is conserved, and k w

†

(., t) k1 remains constant as t ↓ 0. At the same time,

k ∇w

†

(., t) k1 increases monotonically, reflecting the gradual sharpening of edges and

other localized singularities. Such visual and statistical monitoring enable selection of

an optimal image w

†

(x, y, t), which is found at some t > 0. Terminating continuation

at t > 0, is equivalent to redefining p(ρ) in (23) to be p(ρ) = (1 − t)p(ρ), and

then selecting the image at t = 0 as optimal. Clearly, t corresponds to the unknown

constant c in (7).

In practice, the true system otf is seldom found, and the image L

1 norm may

show a modest increase as t decreases from t = 1 to t = t. Conservation of L

1 norm
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 9

Levy fit

Linnik

 fit

Linnik

Levy

Behavior near the origin

Least squares fit in M51

Fig. 2. Blurred M51 image data log |gˆ(ξ, 0)| can be well-fitted with a L ́evy stable otf, as well

as with a generalized Linnik otf. Both otfs have similar behavior away from the origin, but differ

fundamentally near the origin. This difference will play a crucial role.
10 ALFRED S. CARASSO

Original M51 Levy deblur Linnik deblur

Lipschitz exponents in Whirlpool image

Fig. 3. Logarithmic diffusion equation produces higher quality reconstruction of Whirlpool

galaxy image than does fractional diffusion equation. Such differences in quality are reflected in the

smallness of the image Lipschitz exponent, which equals twice the slope of the corresponding Σ line.

Here, Original image has Lip α = 0.39, L ́evy image has Lip α = 0.21, and Linnik image has Lip

α = 0.13.

in w

†

(x, y, t) can be enforced for any desired t, by rescaling w

†

(x, y, t) to the value

k g k1. Moreover, as emphasized in [7, Figure 1], given a blurred g(x, y), there may

be infinitely many distinct otfs hˆ(ξ, η) that can competently deblur that image.

5. Logarithmic vs fractional diffusion in Whirlpool galaxy image. We

shall now demonstrate the significance of the preceding developments by applying

them to sharpen real blurred images with unknown blurs, subject to real but un-
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 11

known noise processes. Our first example is a 1024 × 1024 pixels image g(x, y) of the

Whirlpool galaxy (M51), taken by Travis Rector and Monica Ramirez at the Kitt

Peak National Observatory (NOAO/AURA/NSF), (National Optical Astronomy Observatory/Associated Universities for Research in Astronomy/National Science Foundation).

In Figure 2, we examine the data log |gˆ(ξ, 0)|. With ρ =

p

ξ

2 + η2, we first consider a L ́evy stable fit using hˆ(ρ) = exp(−δρ2β

), δ > 0, 0 < β ≤ 1. While a Gaussian

otf (β = 1.0) is an allowed possiblity, the monotone convexity property in log |gˆ(ξ, 0)|

results in a least squares fit with δ = 1.07, β = 0.15. Next, with w(x, y, 1) = g(x, y),

we solve the fractional diffusion equation wt = −



b(−∆)β



w, 0 ≤ t ≤ 1, b =

δ(4π

2

)

−β

, backwards in time, using an FFT implementation of (22). With the regularization parameters K = 40.0, s = 0.001, the optimal image was found at t = 0.75.

The L

1 norm was conserved on t ≤ t ≤ 1.0, while the TV norm increased almost

fourfold. The reconstructed image is the middle image at the top of Figure 3. The

value β = 0.15 for the L ́evy exponent, is much lower than what is typically found in

most physical applications of L ́evy stable laws.

The analysis in Section 2 emphasizes the behavior of hˆ(ρ) away from the origin,

and makes it plausible that the system otf might be identifiable from that behavior.

For blurred imagery with monotone convex Fourier decay, low exponent L ́evy stable

laws are natural candidate otfs in blind deconvolution. The fundamental observation

of this paper is that there exist otfs with behavior away from the origin that is

almost identical to that in low exponent L ́evy stable laws, yet with behavior near

the origin that is almost Gaussian. Indeed, consider a generalized Linnik fit to the

same M51 data log |gˆ(ξ, 0)|, using hˆ(ρ) = (1 + γρ2σ)

−λ

, γ, λ > 0, 0 < σ ≤ 1,

as shown in Figure 2. This results in σ = 0.86, γ = 0.09, λ = 0.821. The bottom

drawing in Figure 2 illustrates the major difference between Levy and Linnik otfs

near the origin. We next solve the backwards in time logarithmic diffusion problem,

wt = − [λ log{1 + c(−∆)σ}] w, c = γ(4π

2

)

−σ

, 0 ≤ t ≤ 1, w(x, y, 1) = g(x, y).

Here, with K = 40.0, s = 0.001 in (22), we now find the optimal image at t = 0.65.

Again the L

1 norm was conserved on t ≤ t ≤ 1.0, but the TV norm increased almost

sevenfold. The reconstructed image is the rightmost image at the top of Figure 3.

Clearly, while the L ́evy deblurred image significantly sharpens the original, the

Linnik deblurred image provides higher quality reconstructions of the spiral arms,

dustlanes, and galactic cores. The image Lipschitz exponent, discussed below, is a

valuable image metrology tool that can measure the extent of fine structure recovery

in deblurring, and can quantify any improvement produced using logarithmic rather

than fractional diffusion equations. This tool will be used throughout this paper.

6. Fine structure recovery and image Lipschitz exponents. An image

f(x, y) has L

1 Lipschitz exponent α, if and only if

Z

R2

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

α

(24) , |h| → 0,

where |h| = (h

2

1+h

2

2

)

1/2

, and α is fixed with 0 < α ≤ 1. Most images are not smoothly

differentiable functions of x and y, but display edges, localized sharp features, and

other significant fine scale details or texture. The Lipschitz exponent measures the fine

structure content of an image, provided that image is relatively noise free. The value

of α decreases with increasing fine structure. An image that is of bounded variation,

or smoother, has α = 1. Most natural images have α < 0.6, and are not of bounded

variation. Images of starfields, galaxies, and clusters of galaxies, often have α < 0.3.
12 ALFRED S. CARASSO

In [8], [13], an effective method for estimating image Lipschitz exponents is developed, based on a major theorem in [34]. For fixed τ > 0, define the linear operator

Gτ by means of the Fourier series

G

τ

f =

X∞

ξ,η=−∞

exp{−τ(ξ

2 + η2

)}

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

Let μ(τ) =k Gτ f − f k1 / k f k1, be the L

1

relative error in approximating f with

Gτ f. An image f(x, y) has Lipschitz exponent α if and only if μ(τ) = O(τ

α/2

) as

τ ↓ 0. Using FFT algorithms, (25) can be evaluated for each fixed τn > 0 in a sequence

{τn} tending to zero, together with μ(τn). By plotting μ(τn) vs τn on a log-log scale,

positive constants C and α can be located such that μ(τ) ≤ C τ α/2 as τ ↓ 0.

The three solid curves in the bottom drawing in Figure 3, are plots of μ(τ) vs τ

for each image in Figure 3. Each plot is majorized by a dashed straight line Σ. The

Lipschitz exponent α of each image is equal to twice the slope of the corresponding

Σ line. Accordingly, the original M51 image has Lip α = 0.39. The L ́evy deblurred

image has Lip α = 0.21, indicating significant sharpening of the original. However,

the Linnik deblurred image has Lip α = 0.13, a striking improvement.

7. Tail behavior in physical space and mysterious L ́evy exponents. The

improved reconstruction of M51 using the Linnik otf, is directly traceable to the

difference in behavior near the origin in the bottom drawing in Figure 2. As can

be inferred from that drawing, selecting the L ́evy otf rather than the Linnik otf,

implicitly declares the low frequencies in the blurred image g(x, y) to be significantly

more attenuated than would be implied by the Linnik otf. Consequently, in deblurring

with the fractional diffusion equation, such low frequencies are necessarily amplified

to a larger extent than would occur using the logarithmic diffusion equation. Such

spurious overamplification causes saturation and loss of resolution in the galactic cores

and spiral arms in the L ́evy image. To mitigate this effect, the exit point t in the

L ́evy case must generally be chosen larger than in the Linnik case, which leads to less

sharpening. As shown in Figure 4, such L ́evy image saturation phenomena near the

galactic cores are fairly common. As will be seen below, Linnik otf behavior near the

origin plays a beneficial role in all the blind deconvolution experiments in this paper.

In physical (x, y) space, both the isotropic 2D L ́evy stable and generalized Linnik

probability densities are heavy-tailed densities. However, although the corresponding

characteristic functions can have the almost identical high frequency behavior shown

in Figure 2, the tail behavior in physical (x, y) space is distinctly different. Indeed,

such tail behavior is directly related to the smoothness of the characteristic function

near the origin in Fourier (ξ, η) space. The following important results may be found

in [26].

hˆLevy(ξ, η) = exp(−δρ2β), δ > 0, 0 < β < 1, ρ =

p

(ξ

2 + η2),

hLevy(x, y) = O



R−2−2β



, R ↑ ∞, R =

p

(x

2 + y2).

(26)

ˆhLinn(ξ, η) = (1 + γρ2σ

)

−λ

, γ, λ > 0, 0 < σ < 1, ρ =

p

(ξ

2 + η2),

hLinn(x, y) = O



R−2−2σ



, R ↑ ∞, R =

p

(x

2 + y2).

(27)
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 13

Levy deblur

Orig M61 Orig NGC3310

Levy deblur

Linnik deblur Linnik deblur

Fig. 4. Fractional diffusion processing generally causes saturation and loss of resolution near

galactic cores, due to the behavior of L ́evy stable otfs near the origin, as shown in Figure 2. The

use of logarithmic diffusion avoids this difficulty.
14 ALFRED S. CARASSO

Applied to the M51 image in Figure 3, where δ = 1.07, β = 0.15, σ = 0.86, γ = 0.09,

and λ = 0.821, this leads to

hLevy(x, y) = O



R

−2.3



, hLinn(x, y) = O



R

−3.72 

(28) , R ↑ ∞.

Thus, the Linnik density has a much thinner tail than the L ́evy density in physical

space. In summary, in the Whirlpool galaxy image in Figure 3, the successful generalized Linnik point spread function has the tail behavior of a L ́evy stable density

with exponent β = 0.86 in physical (x, y) space, but the high-frequency behavior of

a stable density with β = 0.15 in Fourier (ξ, η) space. Similar behavior characterizes the successful Linnik otfs in other blind deconvolution experiments in this paper.

This result resolves the unexplained appearance of exceptionally low L ́evy exponents

in the detected successful L ́evy otfs in [10] and [11]. Clearly, such low exponents

are not physically meaningful, but stem from choosing L ́evy stable otfs, rather than

generalized Linnik otfs, to match the blurred image data in [10] and [11].

8. Logarithmic diffusion and Hubble space telescope imagery. It is remarkable that the seemingly naive methodology developed in the previous sections

can be useful in improving Hubble imagery. However, in many cases, faint background

objects can be made more visible, and the structure of foreground objects can become

more clearly defined. This appears to be possible for the third-generation instrument

known as the Advanced Camera for Surveys (ACS), recently repaired in May 2009, as

well as for the ‘workhorse’ Wide Field and Planetary Camera 2, (WFPC2), replaced

with WFPC3 in May 2009. Access to the on-line version of this paper will be helpful

in the subsequent discussion. Credit for the three images to be deconvolved below

includes NASA, (National Aeronautics and Space Administration), ESA, (European

Space Agency), the Hubble Heritage Team (STScI/AURA), (Space Telescope Science

Institute/Association of Universities for Research in Astronomy), A. Riess (STScI),

and K. Noll (STScI).

1. Pinwheel galaxy NGC 1309. The example at the top of Figure 5 is a February

2006 Hubble telescope image of the spiral galaxy NGC 1309. The image was taken

with Hubble’s most powerful camera, the Advanced Camera for Surveys. NASA describes this image as featuring bright bluish clusters of star formation, together with

dust lanes spiraling into a yellowish central core of older-population stars. The image is complemented by far-off background galaxies. This galaxy is also home to a

supernova that can help astronomers measure the expansion rate of the universe.

Considerable enhancement of these features is possible using logarithmic diffusion. The color image was broken up into its constituent red, green, and blue component images (RGB), and each component was treated in turn. For each component, log |gˆ(ξ, 0)| was best-fitted with the generalized Linnik otf expression hˆ(ρ) =

(1 + γρ2σ)

−λ

, γ, λ > 0, 0 < σ ≤ 1. For the blue component, this results in

γ = 0.174, λ = 0.653, and σ = 0.88. Otfs for the other two components were found to

differ only slightly from the blue component otf. The latter was therefore used for all

three components. We next solve the backwards in time logarithmic diffusion problem,

wt = − [λ log{1 + c(−∆)σ}] w, c = γ(4π

2

)

−σ

, 0 ≤ t ≤ 1, w(x, y, 1) = g(x, y), using

(22). In all three cases, with the regularization parameters K = 40.0, s = 0.001, and

with t = 0.7, the L

1 norm was conserved on 1 ≤ t ≤ t, while the TV norm increased

fivefold.

The Linnik deblurred image is shown at the bottom of Figure 5. At the outskirts

of that image, several far-off galaxies have become more visible, and the scythe-shaped
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 15

Original HST NGC1309

 Linnik deblur

Fig. 5. Successful blind deconvolution of Hubble space telescope NGC 1309 image. Original

blue component image has Lip α = 0.25. Linnik blue component has Lip α = 0.09.
16 ALFRED S. CARASSO

Linnik deblur

Original HST NGC6050

Fig. 6. Successful blind deconvolution of Hubble space telescope NGC 6050 image. Original

blue component image has Lip α = 0.28. Linnik blue component has Lip α = 0.15.
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 17

Original HST NGC3310

 Linnik deblur

Fig. 7. Successful blind deconvolution of Hubble space telescope NGC 3310 image. Original red

component image has Lip α = 0.26. Linnik red component has Lip α = 0.14.
18 ALFRED S. CARASSO

Original SEM

Linnik deblur

Fig. 8. Blind deconvolution of scanning electron micrograph of complex crystalline structure

recovers numerous surface particles and other fine detail. Original image has Lip α = 0.37. Linnik

image has Lip α = 0.20.
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 19

Original SEM

Linnik deblur

Fig. 9. Blind deconvolution of scanning electron micrograph of slab like crystalline structure

recovers fine scale detail. Note improvement in corner structures. Original image has Lip α = 0.48.

Linnik image has Lip α = 0.25.
20 ALFRED S. CARASSO

structure at eleven o’clock is now brighter and more clearly defined. Previously barely

visible outer spiral arms have been recovered, along with a multitude of previously

faint stars. Within the galaxy’s main body, the bluish star clusters along the spiral

arms, and the dust lanes, are now much better resolved.

Evidently, significant fine-structure reconstruction was achieved, and this is reflected in Lipschitz exponents. For the blue component, the original image has Lip

α = 0.25, while the Linnik image has Lip α = 0.09. L ́evy deblurring of the blue

component required termination at t = 0.8, and displayed loss of resolution due to

saturation. That image had Lip α = 0.12.

2. Colliding galaxies NGC 6050. The example at the top of Figure 6 involves

a spectacular collision between two spiral galaxies. That image is part of a large collection of similar images released by NASA in April 2008. Acquired using Hubble’s

Wide Field and Planetary Camera 2, the image represents a rarely observed snapshot

of a galactic merger thought to require several hundred million years to complete.

It is believed that the Milky Way and Andromeda galaxies will eventually merge in

similar fashion.

In applying logarithmic diffusion to this image, distinct generalized Linnik otfs

hˆ(ρ) = (1 + γρ2σ)

−λ

, γ, λ > 0, 0 < σ ≤ 1, were obtained for each RGB component.

The (γ, λ, σ) triplet for each component was as follows: red = (0.195, 0.692, 0.85);

green = (0.201, 0.643, 0.88); and blue= (0.151, 0.684, 0.83). We next solve the

backwards in time logarithmic diffusion problem, wt = − [λ log{1 + c(−∆)σ}] w, c =

γ(4π

2

)

−σ

, 0 ≤ t ≤ 1, w(x, y, 1) = g(x, y), using (22). In all three cases, with

K = 100, s = 0.001, and t = 0.75, the L

1 norm was conserved on 1 ≤ t ≤ t, while

the TV norm increased by a factor of about 3.5.

The Linnik image at the bottom of Figure 6 is noticeably brighter, and displays

better resolution of the structural details in the two galaxies. Bluish young star clusters along spiral arms are better defined, together with dust lanes spiraling around the

three cores of older stars. Background galaxies have become more visible. As a result,

significant reduction in Lipschitz exponents was recorded. For the blue component,

the original image has Lip α = 0.28, while the Linnik image has Lip α = 0.15. L ́evy

deblurring of the blue image resulted in Lip α = 0.18.

3. Starburst galaxy NGC 3310. The example at the top of Figure 7 is a September

2001 Hubble telescope image of the Starburst galaxy NGC 3310, acquired using the

Wide Field and Planetary Camera 2. That galaxy is of great interest to astronomers

as it is known to produce clusters of new stars at a prodigious rate. Several hundred

such clusters are visible as bright blue diffuse objects around the galaxy’s spiral arms,

each cluster representing up to a million stars.

In applying logarithmic diffusion to this image, the Linnik otf parameters for the

red component were found to be γ = 0.328, λ = 0.557, and σ = 0.92. The red and

blue otfs almost coincided, while the green otf differed slightly. Accordingly, the red

image otf was used for all three components. We next solve the backwards in time

logarithmic diffusion problem, wt = − [λ log{1 + c(−∆)σ}] w, c = γ(4π

2

)

−σ

, 0 ≤ t ≤

1, w(x, y, 1) = g(x, y), using (22). In all three cases, with K = 40.0, s = 0.001, and

t = 0.75, the L

1 norm increased by about 3% on 1 ≤ t ≤ t, with a threefold increase

in TV norm.

The Linnik deblurred image at the bottom of Figure 7 shows significant recovery

of blue clusters of young stars at the periphery of the image; also, there is noticeable
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 21

enhancement of structural detail in the main body of the galaxy, and in the surrounding dust lanes. In contrast, L ́evy deblurring tends to produce saturation and loss of

resolution near the core, as previously shown in Figure 4. For the red component, the

original image has Lip α = 0.26, while the Linnik image has Lip α = 0.14. The L ́evy

image has Lip α = 0.16.

9. Scanning electron micrographs. There are other important classes of images which display the monotone convex behavior described in section 2, and for

which logarithmic diffusion deconvolution may be useful. One such class is scanning

electron micrographs. In this section, we revisit experiments previously described in

[11], and we verify that Linnik otfs can improve substantially on previous results. The

two images considered here were taken by John Small at NIST. They are micrographs

of a complex multiform crystalline compound of mercury. The author expresses his

gratitude to Dr. David S. Bright (NIST) for making these images available.

The example at the top of Figure 8 has complex crystalline structure and fine scale

surface detail that are of interest. The Linnik otf parameters for log |gˆ(ξ, 0)| were found

to be γ = 0.0266, λ = 0.5496, and σ = 0.89. We next solve the backwards in time

logarithmic diffusion problem, wt = − [λ log{1 + c(−∆)σ}] w, c = γ(4π

2

)

−σ

, 0 ≤

t ≤ 1, w(x, y, 1) = g(x, y), using (22). With K = 50.0, s = 0.001, and t = 0.6, the

L

1 norm was conserved on 1 ≤ t ≤ t, with a threefold increase in TV norm. The

Linnik image at the bottom of Figure 8 shows recovery of numerous surface particles

and other detail, and the complex morphology has become better defined. Deblurred

image quality is particularly noticeable in this example. The original image has Lip

α = 0.37, and the Linnik image has Lip α = 0.20. The L ́evy image was terminated at

t = 0.7 and had Lip α = 0.25.

The slab like structure at the top of Figure 9 also exhibits interesting surface

detail. The Linnik otf parameters were found to be γ = 0.022, λ = 0.776, and

σ = 0.86. With K = 50.0, s = 0.001, and t = 0.7, the L

1 norm was conserved on

1 ≤ t ≤ t, with a threefold increase in TV norm. Here again, slab surface details are

well-recovered, and the two corner structures reveal interesting small scale granularity.

The original image has Lip α = 0.48, and the Linnik image has α = 0.25. The L ́evy

image was terminated at t = 0.75 and had Lip α = 0.34.

Table 1 below summarizes the changes in image Lipschitz exponents recorded in

each of the six deconvolution experiments in this paper. Linnik deblurring consistently

produces smaller Lipschitz exponents than does L ́evy deblurring. The last column in

Table 1 is especially noteworthy.

TABLE 1

Lipschitz exponents before and after deconvolution

Image Original Lip L ́evy Lip Linnik Lip Linnik/Orig

M 51 α = 0.39 α = 0.21 α = 0.13 33 %

NGC 1309 α = 0.25 α = 0.12 α = 0.09 36 %

NGC 6050 α = 0.28 α = 0.18 α = 0.15 53 %

NGC 3310 α = 0.26 α = 0.16 α = 0.14 54 %

SEM Crystal α = 0.37 α = 0.25 α = 0.20 54 %

SEM Slab α = 0.48 α = 0.34 α = 0.25 52 %
22 ALFRED S. CARASSO

10. Concluding remarks. A priori knowledge about the solution is an essential

element in the successful computation of ill-posed inverse problems. Such knowledge

informs the construction of the blind deconvolution procedure discussed in this paper.

This methodology was shown to be capable of producing credible reconstructions in

two distinct classes of real blurred images with real noise, one at the nanoscale, and

the other at the cosmological scale. These two classes are acquired using distinct

imaging modalities and are of considerable scientific interest. The accompanying

sizeable reduction in Lipschitz exponents is highly significant.

This work substantially improves on previous work on the same classes of images

[10], [11], based on isotropic L ́evy stable otfs. Low L ́evy exponents were found necessary to accomodate the Fourier behavior in the blurred image data. As shown in

section 3, L ́evy stable motions are related to Brownian motion through subordination.

The present work is based on a new class of otfs that have the same high frequency

behavior as low exponent stable otfs, but behave like high exponent stable laws near

the origin. Such generalized Linnik otfs result from subordination of L ́evy stable motions by the Gamma process, and their behavior near the origin was shown to play a

vital role in fine structure recovery. The compelling quality of the new reconstructions

strongly suggests that the detected generalized Linnik otfs that produce these results,

must emulate essential aspects of the true system otfs in the two classes of images.

Linnik otfs are not currently known in image analysis.

Sixty years after its introduction in [5], Bochner’s seminal mathematical idea continues to yield a rich harvest of important applications.

REFERENCES

[1] C. Albanese and S. Lawi Time quantization and q-deformations, J. Phys. A: Math. Gen. 37

(2004), pp. 2983–2987.

[2] L. Bar, N. Sochen, and N. Kiryati, Semi-blind image restoration via Mumford-Shah regularization, IEEE Transactions on Image Processing, 15 (2006), pp. 483–493.

[3] P. Barthelemy, J. Bertolotti, and D. S. Wiersma, A L ́evy flight for light, Nature 453

(2008), pp. 495-498.

[4] S. Bochner, Monotone Funktionen, Stieltjes Integrale und harmonische Analyse, Math. Ann.,

108 (1933), pp. 378-410.

[5] S. Bochner, Diffusion equation and stochastic processes, Proc. Nat. Acad. Sci. USA 35 (1949),

pp. 368-370.

[6] S. Bochner, Harmonic Analysis and the Theory of Probability, University of California Press,

Berkeley, CA 1955.

[7] A. S. Carasso, The APEX method in image sharpening and the use of low exponent L ́evy

stable laws, SIAM J. Appl. Math. 63 (2002), pp. 593–618.

[8] A. S. Carasso, Singular integrals, image smoothness, and the recovery of texture in image

deblurring, SIAM J. Appl. Math. 64 (2004), pp. 1749-1774.

[9] A. S. Carasso, False characteristic functions and other pathologies in variational blind deconvolution. A method of recovery, SIAM J. Appl. Math. 70 (2009), pp. 1097–1119.

[10] A. S. Carasso, APEX blind deconvolution of color Hubble space telescope imagery and other

astronomical data, Optical Engineering 45 (2006), 107004.

[11] A. S. Carasso, D. S. Bright, and A. E. Vladar ́ , APEX method and real-time blind deconvolution of scanning electron microscope imagery, Optical Engineering 41 (2002), pp.

2499-2514.

[12] A. S. Carasso and T. Kato, On subordinated holomorphic semigroups, Trans. Amer. Math.

Soc. 327 (1991), pp. 867–878.

[13] A. S. Carasso and A. E. Vladar ́ , Calibrating image roughness by estimating Lipschitz exponents, with appplications to image restoration, Optical Engineering 47 (2008), 037012.

[14] T. F. Chan and C. K. Wong, Total variation blind deconvolution, IEEE Transactions on

Image Processing 7 (1998), pp. 370-375.
LOGARITHMIC DIFFUSION EQUATIONS AND BLIND DECONVOLUTION 23

[15] I. Eliazar and J. Klafter, L ́evy, Ornstein-Uhlenbeck, and subordination: spectral vs jump

description, Journal of Statistical Physics 119 (2005), pp. 165-196.

[16] M. B. Erdogan and I. V. Ostrovskii, Analytic and asymptotic properties of generalized

Linnik probability densities, J. Math. Anal. Appl. 217 (1998), pp. 555–578.

[17] G. Fasshauer, Meshfree Approximation Methods with MATLAB, Interdisciplinary Mathematical Sciences - Vol. 6, World Scientific Publishers, Singapore, 2007.

[18] W. Feller, An Introduction to Probability Theory and its Applications, Vol. 2, Second ed.,

Wiley, New York, 1971.

[19] V. Ganti, A. Singh, P. Passalacqua, and E. Foufoula-Georgiou, Subordinated Brownian

motion model for sediment transport, Physical Review E 80 (2009), 011111 pp 1–9.

[20] H. Geman, D. P. Madan, and M. Yor, Time changes for L ́evy processes, Mathematical

Finance 11 (2001), pp. 79-96.

[21] C. B. Johnson, Classification of electron-optical device modulation transfer function, Advances

in Electronics and Electron Physics 33B (1972), pp. 579–584.

[22] C. B. Johnson, MTFs: a simplified approach, Electro-Optical Systems Design 4 (1972), pp.

22–26.

[23] C. B. Johnson, Point-spread functions, line-spread functions, and edge-response functions

associated with mtfs of the form exp[−(ω/ωc)

n], Applied Optics, 12 (1973), pp. 1031–

1033.

[24] L. Justen and R. Ramlau, A non-iterative regularization approach to blind deconvolution,

Inverse Problems 22 (2006), pp. 771–800.

[25] J. Lee and D. Sheen, F. John’s stability conditions versus A. Carasso’s SECB constraint for

backward parabolic problems, Inverse Problems 25 (2009), 055001.

[26] S. C. Lim and L. P. Teo, Analytic and asymptotic properties of multivariate generalized Linnik

probability densities, arXiv:0903.5344v1 [math.PR] 30 Mar 2009.

[27] E. Lukacs, Characteristic Functions, Second ed., Griffin, London, 1970.

[28] D. P. Madan and E. Seneta, The Variance Gamma model for share market returns, Journal

of Business 63 (1990), pp. 511–524.

[29] D. P. Madan, P. P. Carr, and E. C. Chang, The Variance Gamma process and option

pricing, European Finance Review 2 (1998), pp. 79-105.

[30] R. S. Phillips, On the generation of semigroups of linear operators, Pacific J. Math. 2 (1952),

pp. 343–369.

[31] G. Samorodnitsky and M. S. Taqqu, Stable Non-Gaussian Random Processes: Stochastic

Models with Infinite Variance, Chapman and Hall, New York, 1994.

[32] K. Sato, L ́evy Processes and Infinitely Divisible Distributions, Cambridge University Press,

1999.

[33] I. J. Schoenberg, Metric spaces and completely monotone functions, Ann. of Math. 39 (1938),

pp. 811-841.

[34] M. H. Taibleson, On the theory of Lipschitz spaces of distributions on Euclidean n-space. I.

Principal properties, Journal of Mathematics and Mechanics 13 (1964), pp. 407–478.

[35] A. Weron and M. Magdziarz, Anomalous diffusion and semimartingales, Europhysics Letters

86 (2009), 60010 pp. 1–6.