NISTIR 6570
 

Direct Blind Deconvolution II. 

Substitute Images and the BEAK 

Method 

Alfred S. Carasso 

U. S. DEPARTMENT OF COMMERCE 

Technology Administration 

Mathematical and Computational 

Sciences Division 

National Institute of Standards 

and Technology 

Gaithersburg, MD 20899 

h) 

o 

o 

CENTENNIAL. 

r , 

NISr
 

National 1....1tu•• of S.anelarels 

anel Technology 

Technology Administration 

U.S. Department of Commerce 
NISTIR 6570
 

Direct Blind Deconvolution II. 

Substitute Images and the BEAK 

Method 

Alfred S. Carasso 

u. s. DEPARTMENT OF COMMERCE 

Technology Administration 

Mathematical and Computational 

Sciences Division
 

National Institute of Standards
 

and Technology
 

Gaithersburg, MD 20899
 

November 1, 2000 

u.S. DEPARTMENT OF COMMERCE 

Norman Y. Mineta, Secretary 

TECHNOLOGY ADMINISTRATION 

Dr. Cheryl L. Shavers, Under Secretary 

of Commerce for Technology 

NATIONAL INSTITUTE OF STANDARDS 

AND TECHNOLOGY 

Raymond G. Kammer, Director 
DIRECT BLIND DECONVOLUTION II.

SUBSTITUTE IMAGES AND THE BEAK METHOD

ALFRED S. CARASSO∗

Abstract. The BEAK method is an FFT-based direct blind deconvolution technique previously

introduced by the author, and applied to a limited but significant class of blurs that can be expressed

as convolutions of two-dimensional radially symmetric L ́evy probability density functions. This class

includes and generalizes Gaussian and Lorentzian distributions, but does not include defocus blurs.

The method requires a-priori information on the Fourier transform ˆfe(ξ, η) of the unknown exact

image fe(x, y), namely, the gross behavior of log |

ˆfe(ξ, η)| along a single line through the origin in

the (ξ, η) plane. The present paper significantly extends the applicability of the BEAK method. It

is shown that images of similar objects often display approximately equal gross behavior, and that

gross behavior in such substitute images can be used successfully in numerous practical contexts.

Next, using substitute images, a variant of the BEAK method is developed that can handle defocus

blurs. The paper is illustrated with several examples of blind deconvolution of 512 × 512 images in

the presence of noise, and includes a detailed discussion of an example where the BEAK method

fails.

Key words. image deblurring, blind deconvolution, direct methods, L ́evy density functions,

defocusing, substitute images, BEAK method, SECB method.

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

1. Introduction. This paper is a sequel to [6]. As was the case there, the

procedure described below will generally not be useful for severely blurred images

at high levels of noise, nor for images degraded by arbitrary or unknown processes.

Rather, the method is limited to a narrow class of deblurring problems involving

restricted types of blur at low to moderate intensities, in the presence of low levels of

noise. With the type of blur assumed known, these conditions enable identification of

the parameters in the system point spread function. The discussion below includes

an example where the method fails. Other examples of failure are easily found. All

images in this paper are of size 512 × 512. Unless otherwise indicated, these images

are quantized at 8-bits per pixel, i.e., each pixel value is an integer lying between 0

and 255. Such quantization introduces 8-bit rounding noise which plays a significant

role. Additional noise processes are applied in some cases.

Blind deconvolution seeks to deblur an image without knowing the point spread

function describing the blur. In [6], two methods were developed for direct (i.e., noniterative) blind deconvolution, the BEAK method and the APEX method. These

are Fourier domain techniques for detecting the signature of the system point spread

function from one-dimensional (1-D) analysis of the blurred image. A separate direct

Fourier domain deblurring technique, the SECB method [4], uses this detected point

spread function to deblur the image. Direct blind deconvolution of 512 × 512 images

can be accomplished in minutes on current desktop workstations. Other approaches

to that problem are iterative in nature, but the iterative process is not always wellbehaved [6]. When that process is stable, several thousand iterations may be necessary

to achieve convergence [11], [15]. The discussion in [4] and [6] is focused exclusively on

so-called class G point spread functions. This is a significant but limited class of blurs

that can be expressed as convolutions of 2-D radially symmetric L ́evy ‘stable’ probability density functions. This class includes and generalizes Gaussian and Lorentzian

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

1
2 ALFRED S. CARASSO

distributions, but does not include defocus blurs. The BEAK detection method uses

a novel type of a-priori information about the exact sharp image fe(x, y), namely,

the gross behavior of its Fourier transform ˆfe(ξ, η) along a single line through the

origin in the (ξ, η) plane. Such gross behavior is represented by a smooth function

that provides a least squares fit to the highly oscillatory trace of log |

ˆfe(ξ, η)| along

that line. By pre-rotating the image if necessary, one may assume that when such

gross behavior is known, it is known along the horizontal axis (ξ, 0). A large class of

sharp images, the class W, is exhibited in [6] having the property that gross behavior

can be summarized by the two positive numbers a, b in the expression

log |

ˆfe

∗

(ξ, 0)| ≈ −a |ξ|

b

(1) , a, b > 0, fe(x, y) ∈ W.

where

|

ˆfe

∗

(ξ, 0)| = |

ˆfe(ξ, 0)| /ˆ (2) fe(0, 0).

The second detection method discussed in [6], the APEX method, does not require

prior information on fe(x, y) but assumes that image to be a recognizable object.

Several (fast) interactive trials are necessary before locating a suitable point spread

function using that method.

In the present paper, the APEX method is not discussed further. Rather, the

range of applicability of the BEAK method is extended in two major ways. First,

we observe that gross behavior in class W images provides only superficial prior

information about fe(x, y), and that images of similar objects are often found to

display approximately equal gross behavior. Therefore, in using the BEAK method

for identifying system point spread functions, it is generally not necessary to know the

gross behavior in the original image fe(x, y), but only that in fs(x, y), an appropriate

sharp image of a similar object. There are numerous practical contexts where such

similar images are available. By substituting the gross behavior in fs(x, y) for that

in the unknown fe(x, y), the BEAK method becomes applicable in a wide variety of

situations. This is illustrated with examples of class W images blurred by class G

point spread functions.

A second aim of this paper is to develop a detection method for defocus blurs.

As will be seen below, deconvolution of defocus blurs is strikingly different from class

G blurs. Using substitute images as described above, a variant of the BEAK method

is presented that can approximately identify defocus point spread functions, provided

the defocus is not too severe. This detected psf can then be used to deblur the image.

Since the SECB deblurring procedure used in [6] is restricted to class G blurs, a

modification of that procedure is necessary for defocus and other shift-invariant blurs

not in G. This modification is discussed in the Appendix.

2. Gross behavior in class W images. A-priori information is paramount in

the solution of ill-posed inverse problems. Such information becomes ever more critical

in the case of blind deconvolution, where severe non-uniqueness is compounded with

discontinuous data dependence. In iterative blind deconvolution algorithms, positivity

and support constraints on the convolution components are generally imposed in order

to reduce the multiplicity of solutions. However, such constraints are not always

effective in avoiding traps at local minima (stagnation points), or divergence of the

iterative procedure. There is therefore considerable interest in finding types of prior

information that are widely applicable, and lead to reliable algorithms. Gross behavior

in Fourier space is one example of useful prior knowledge.
BLIND DECONVOLUTION II 3

A B

Fig. 1. log |

ˆfe

∗

(ξ, 0)| ≈ −a |ξ|

b and rough equivalence of gross behavior in two different

spacecraft images. (A) Mariner 5 has a = 2.81, b = 0.190 . (B) Mariner 10 has a = 2.98, b =

0.183.

A B

Fig. 2. log |fˆe

∗

(ξ, 0)| ≈ −a |ξ|

b and rough equivalence of gross behavior in two different

aircraft carrier images. (A) USS Eisenhower has a = 2.91, b = 0.157 . (B) USS Kittyhawk has

a = 3.04, b = 0.158.

In the simplest case, image deblurring is associated with the solution of twodimensional convolution equations

Hf ≡

Z

R2

(3) h(x − u, y − v)f(u, v)dudv ≡ h(x, y) ⊗ f(x, y) = g(x, y),

where g(x, y) is the recorded blurred image, f(x, y) is the desired unblurred image,

h(x, y) is the blurring kernel or point spread function (psf) of the imaging process,

and ⊗ denotes convolution. The psf h(x, y) represents the cumulative effects of all

distortions caused by the media through which signals propagate, as well as all optical

and electronic aberrations produced by imperfect sensing and recording equipment.

It is assumed that h(x, y) is such that the linear problem Hf = g has at most one
4 ALFRED S. CARASSO

A B

Fig. 3. log |

ˆfe

∗

(ξ, 0)| ≈ −a |ξ|

b and rough equivalence of gross behavior in two different face

images. (A) Marilyn Monroe has a = 3.57, b = 0.158 . (B) Ingrid Bergman has a = 3.11, b =

0.183. Solid curve in (A) deviates only slightly from solid curve in (B), despite differences in

corresponding values for a and b.

solution. This is the case for class G point spread functions defined in Fourier space

by

hˆ(ξ, η) = e

−

PJ

i=1

αi(ξ

2+η2

)

βi

(4) , αi ≥ 0, 0 < βi ≤ 1.

The function hˆ(ξ, η) is known as the optical transfer function (otf) of the imaging

process. The blurred image g(x, y) includes noise, which is viewed as a separate

additional degradation,

(5) g(x, y) = ge(x, y) + n(x, y),

where ge(x, y) is the blurred image that would have been recorded in the absence of

noise, and n(x, y) represents the cumulative effects of all noise processes and other

errors affecting final acquisition of the digitized array g(x, y). This includes nonlinear

noise processes where n(x, y) may be a function of f(x, y). Both ge(x, y) and n(x, y)

are unknown, but n(x, y) may be presumed small. The unique solution of (3) when

the right hand side is ge(x, y), is the exact sharp image denoted by fe(x, y). Thus

(6) h(x, y) ⊗ fe(x, y) = ge(x, y).

Since fe(x, y) ≥ 0

|

ˆfe(ξ, η)| ≤ Z

R2

fe(x, y)dxdy = ˆ (7) fe(0, 0) = σ > 0.

Also, since h(x, y) is a probability density,

gˆe(0, 0) = Z

R2

ge(x, y)dxdy =

Z

R2

fe(x, y)dxdy = ˆ (8) fe(0, 0) = σ > 0.

Using σ as a normalizing constant, we may normalize Fourier transform quantities

qˆ(ξ, η) by dividing by σ. Let

qˆ

∗

(9) (ξ, η) = ˆq(ξ, η)/σ,
BLIND DECONVOLUTION II 5

A B

C D

Fig. 4. log |

ˆfe

∗

(ξ, 0)| ≈ −a |ξ|

b and gross behavior in three different MRI brain images.

(A) Brain1 has a = 3.12, b = 0.155. (B) Brain2 has a = 3.28, b = 0.154. (C) Brain3 has

a = 4.03, b = 0.135. (D) Gross behavior roughly equivalent in Brain1 and Brain2, but markedly

different in Brain1 and Brain3. Use of Brain2 as substitute for Brain1 in BEAK method leads to

significantly better results than does use of Brain3. See Figure 8.

denote the normalized quantity. The function |

ˆfe

∗

(ξ, η)| is highly oscillatory, and

0 ≤ | ˆfe

∗

| ≤ 1. Since fe(x, y) is real, its Fourier transform is conjugate symmetric.

Therefore, the function |

ˆfe

∗

(ξ, η)| is symmetric about the origin on any line through

the origin, η = ξ tan θ, in the (ξ, η) plane.

As previously noted, all images in this paper are of size 512×512. For each sharp

image fe(x, y), the discrete Fourier transform is a 512 × 512 array of complex numbers, which we again denote by ˆfe(ξ, η) for simplicity. The ‘frequencies’ ξ, η are now

integers lying between −256 and 256, and the zero frequency is at the center of the

transform array. This ordering is achieved by pre-multiplying fe(x, y) by (−1)x+y.

We shall be interested in values of such transforms along single lines through the

origin in the (ξ, η) plane. The discrete transforms ˆfe

∗

(ξ, 0), and ˆfe

∗

(0, η) are immediately available. Image rotation may be used to obtain discrete transforms along other

directions. All 1-D Fourier domain plots shown in this paper are taken along the axis

η = 0 in the (ξ, η) plane. In these plots, the zero frequency is at the center of the
6 ALFRED S. CARASSO

horizontal axis and the graphs are symmetric about the vertical line ξ = 0.

The class W of well-behaved sharp images fe(x, y) introduced in [6] may be loosely

characterized as follows:

1. |

ˆfe

∗

(ξ, η)| has at most isolated zeroes in the (ξ, η) plane.

2. Neglecting isolated singularities, the global behavior of log |

ˆfe

∗

(ξ, η)| on any line

η = ξ tan θ is roughly monotone decreasing with increasing r = (ξ

2 + η2

)

1/2

; i.e., on

the ray reiθ, a least squares fit to log |

ˆfe

∗

(ξ, η)| with an appropriate monotone decreasing function, provides a fair representation of gross data behavior as r increases.

3. While the rate of decay may vary between rays, this decay is relatively slow, i.e., of

a general order of magnitude comparable to that found in the sharp Mariner images

in Fig 1.

4. The gross behavior of log |

ˆfe

∗

(ξ, η)| along the ray reiθ is defined to be the function

v(r) = −a(θ) r

b(θ)

, a, b > 0, that best fits log |

ˆfe

∗

(ξ, ξ tan θ)| in the least squares

sense.

Nine examples of class W images are displayed in Figures 1 through 4, together with

their gross behavior along the axis η = 0 in the (ξ, η) plane. Several other examples

are given in [6]. What is of interest here, is the fact that ‘similar’ objects appear to

have roughly equivalent gross behavior. In Figures 1 and 2, the values for a and b are

approximately equal in the two spacecraft and in the two carrier images. Hence, the

corresponding gross behavior traces are roughly equivalent. In Figure 3, the values of

a and b for Marilyn Monroe and Ingrid Bergman are substantially different. However,

the corresponding gross behavior traces almost coincide in that case. In Figure 4,

gross behavior in Brain3 is markedly different from that in Brain1 and Brain2. This

shows that similar objects need not always have equivalent gross behavior. In practice, it may be necessary to sample several similar objects before arriving at a good

substitute image for use in the BEAK method.

3. The BEAK method for class G psfs. As discussed in [6], this is a Fourier

domain technique for detecting class G point spread functions acting on class W

images. The method uses 1-D Fourier analysis of the blurred image data g(x, y)

in (3), and requires prior knowledge of gross behavior in the unknown sharp image

fe(x, y) in (6). The detected optical transfer function is then used to solve the ill-posed

problem (3). This is accomplished using another direct Fourier domain procedure, the

SECB method [2, 3, 4].

The BEAK method is based on the following observations. In the basic relation

(10) g(x, y) = h(x, y) ⊗ fe(x, y) + n(x, y),

we may safely assume that the noise n(x, y) satisfies

Z

R2

|n(x, y)|dxdy ≪

Z

R2

(11) fe(x, y)dxdy = σ > 0,

so that,

|nˆ

∗

(12) (ξ, η)| ≪ 1.
BLIND DECONVOLUTION II 7

A B

Fig. 5. Application of BEAK method to 8-bit blurred Eisenhower image, using gross behavior

from sharp Kittyhawk image. (A) 8-bit blurred Eisenhower image g(x, y) obtained by convolution

of sharp Eisenhower image fe(x, y) with ‘long exposure turbulence’ L ́evy psf with α = 0.003, β =

5/6. Trace of log |gˆ

∗

(ξ, 0)| above 8-bit noise level only on small interval −60 ≤ ξ ≤ 60. (B)

log |fˆs

∗

(ξ, 0)| ≈ −3.04 |ξ|

0.158 in sharp Kittyhawk image will be used to obtain BEAK fit to blurred

Eisenhower trace.

Consider the case where the otf is a pure L ́evy density hˆ(ξ, η) = e

−α(ξ

2+η2

)

β

. Since

g = ge + n

log |gˆ

∗

(ξ, η)| = log |e

−α(ξ

2+η2

)

β

ˆfe

∗

(ξ, η) + ˆn

∗

(13) (ξ, η)|.

Let Ω = {(ξ, η) | ξ

2 + η2 ≤ ω2} be a neighborhood of the origin where

e

−α(ξ

2+η2

)

β

|

ˆfe

∗

(ξ, η)| ≫ |nˆ

∗

(14) (ξ, η)|.

Such an Ω exists since (14) is true for ξ = η = 0 in view of (12). The radius ω > 0 of

Ω decreases as α and n increase. For (ξ, η) ∈ Ω we have

log |gˆ

∗

(ξ, η)| ≈ −α(ξ

2 + η2

)

β + log |

ˆfe

∗

(15) (ξ, η)|.

Hence, for |ξ| ≤ ω,

log |gˆ

∗

(ξ, 0)| ≈ −α|ξ|

2β + log |

ˆfe

∗

(16) (ξ, 0)|.

The idea in [6] is to replace the unknown log |

ˆfe

∗

(ξ, 0)| in (16) by its gross

behavior, i.e., its least squares approximation v(ξ) = −a |ξ|

b

, which is assumed

known. Extensive numerical experiments with a wide assortment of images and class

G psfs, indicate this to be a successful strategy. In the present paper, we replace the

unknown log |

ˆfe

∗

(ξ, 0)| in (16) by the least squares approximation vs(ξ) = −a |ξ|

b

to log |

ˆfs

∗

(ξ, 0)|, where fs(x, y) is a substitute image, i.e., an image of a similar object

that is expected to have roughly equivalent gross behavior. Several candidate images

may need to be tried before achieving optimal results. The detection procedure is

then the following. With a and b given a-priori, find positive numbers  ̃α, β,  ̃

so that the function u(ξ) = −α ̃ |ξ|

2β ̃ − a |ξ|b best fits log |gˆ∗

(ξ, 0)| on |ξ| ≤ ω.
8 ALFRED S. CARASSO

A B

C D

Fig. 6. Blind deblurring of aircraft carrier image using BEAK psf detection method. (A) BEAK

fit to 8-bit blurred Eisenhower trace on −60 ≤ ξ ≤ 60, using gross behavior in sharp Kittyhawk trace.

Fit returns detected parameters α = 0.000931, β = 0.97, differing from true values α = 0.003, β =

5/6. (B) Detected psf (dashed) roughly approximates true psf (solid). (C) Blurred Eisenhower

image. (D) SECB deblurring of image C using s = 0.01, K = 1.27, and detected psf parameters

α = 0.000931, β = 0.97. Good restoration achieved even though detected psf not equivalent to true

psf.

This may be accomplished interactively using nonlinear least squares algorithms in

DATAPLOT [10]. The returned values for  ̃α and β ̃ are subsequently used for hˆ(ξ, η)

in the SECB deblurring procedure.

For more general class G otfs where hˆ(ξ, η) = e

−Σ

N

i=1αi(ξ

2+η2

)

βi

, we again seek the

best fit to log |gˆ

∗

(ξ, 0)| on |ξ| ≤ ω, with a function u(ξ) = −α ̃ |ξ|

2β ̃ − a |ξ|b

. Here,

the returned values for  ̃α and β ̃ may be considered average values for the αi, βi, and

are expected to generate a pure L ́evy density that well-approximates the composite

psf.

We illustrate this procedure in Figures 5 and 6. In Figure 5(A), the sharp USS

Eisenhower image shown in Figure 2(A), was artificially blurred by convolution with a

pure L ́evy density with α = 0.003 and β = 5/6. This simulates long exposure imaging

in the presence of turbulence [9]. The result of that numerical convolution was then

truncated to 8-bits. The trace of log |gˆ

∗

(ξ, 0)| in Figure 5(A) lies above 8-bit noise level

only on the small interval −60 ≤ ξ ≤ 60, and is obviously quite different from that

in the original image in Figure 2(A). The sharp USS Kittyhawk image in Figure 5(B)
BLIND DECONVOLUTION II 9

A B C

Fig. 7. Application of BEAK method to 8-bit Brain1 image blurred with L ́evy density where

α = 0.05 and β = 0.6. Blurred image trace lies above 8-bit noise level on interval −50 ≤ ξ ≤ 50. (A)

BEAK fit using gross behavior in sharp Brain2 image returns detected parameters α = 0.0325, β =

0.643, differing from true values. (B) BEAK fit using gross behavior in sharp Brain3 image develops

more curvature, and returns α = 0.00316, β = 0.921, markedly different from true values. (C)

Detected psf in (A) (dotted curve) is much closer to true psf (solid curve), than is detected psf in

(B) (dashed curve).

A B C

Fig. 8. Blind deblurring of Brain1 image using two distinct BEAK detected psfs. (A) 8-bit

blurred Brain1 image. (B) SECB deblurring of image A using s = 0.001, K = 1.27, and Brain2

detected psf parameters α = 0.0325, β = 0.643. Good restoration achieved even though detected psf

not equivalent to true psf. (C) SECB deblurring of image A using s = 0.001, K = 1.27, and Brain3

detected psf parameters α = 0.00316, β = 0.921, results in comparatively weak restoration.

will now be used as the substitute image in the BEAK method. Using a = 3.04 and

b = 0.158, we best fit log |gˆ

∗

(ξ, 0)| in Figure 5(A) on the interval |ξ| ≤ 60, with the

function u(ξ) = −α ̃ |ξ|

2β ̃ − a |ξ|b

, using the fit command in DATAPLOT [10]. That

fit is indicated by the solid curve in Figure 6(A), and the beak or nib in that curve near

ξ = 0 is a characteristic feature of this detection procedure. The fit returns detected

parameters  ̃α = 0.000931, β ̃ = 0.97, that differ substantially from the correct values

α = 0.003, β = 5/6. However, as indicated in Figure 6(B), the detected psf (dashed

curve) roughly approximates the true psf (solid curve). Using this detected psf in

the SECB method with s = 0.001 and K = 1.27, produces a quite useful restoration,

as shown in Figure 6(D). While ringing artifacts are visible in the deblurred image

against the light sky background near the top of the image, considerable sharpening

of the planes on deck, as well as the ‘island’ structure containing the bridge and mast,

has been achieved.

The next example illustrates the importance of locating a good substitute image.
10 ALFRED S. CARASSO

The sharp MRI Brain1 image in Figure 4(A) was blurred synthetically by convolution

with a pure L ́evy density with α = 0.05 and β = 0.6, and the numerical convolution

was rounded to 8-bits. This time, the trace of log |gˆ

∗

(ξ, 0)| lies above 8-bit noise level

on the interval −50 ≤ ξ ≤ 50. In Figure 7(A), the BEAK method is applied on |ξ| ≤ 50

using Brain2 in Figure 4(B) as a substitute image. This produces detected parameters

α = 0.0325, β = 0.643. In Figure 7(B), the BEAK method is applied using Brain3

in Figure 4(C) as the substitute image. The resulting fit develops more curvature

than in Figure 7(A), and returns detected parameters α = 0.00316, β = 0.921. These

values are markedly different from the true values α = 0.05, β = 0.6. As shown in

Figure 7(C), and as may be expected from Figure 4(D), the detected psf using Brain2

is much closer to the true psf, than is the detected psf using Brain3. In Figure 8,

we see that the Brain2 detected psf produces a quite good restoration, whereas the

narrower Brain3 detected psf results in a comparatively weak reconstruction.

Numerous other experiments with class G psfs using the BEAK method with

substitute images, confirm the pattern of behavior typified by the above two examples.

It should be noted that in addition to locating a good substitute image, successful

detection requires having a sufficiently wide interval about ξ = 0 wherein the trace of

log |gˆ

∗

(ξ, 0)| lies above noise level. This was the case in Figures 6(A) and 7(A). Severe

blurring caused by large values of α when β is near unity, and/or high levels of noise,

considerably reduce the width of that interval and can preclude useful detection.

4. Uniform defocus blur. With class G blurs, one can often deblur the image

quite satisfactorily with a psf that is only a rough approximation to the true psf. This

was the case in the USS Eisenhower image in Figure 6(D), the MRI brain image in

Figure 8(B), as well as in several other examples in [6]. We shall find defocus blurs

to be less forgiving.

Defocus blurs are discussed in [14]. Important early work on blind deconvolution

of defocus blurs can be found in [1]. Additional references may be found in [12]. If

R > 0 is the radius of the ‘circle of confusion’, the psf for uniform defocus blur is

given by

h(x, y) =







(πR2)

−1

, x2 + y

2 ≤ R2

,

0, x2 + y

2 > R2

.

(17)

This has a Fourier transform given by the ‘sombrero function’ [8, p. 72]

hˆ(ξ, η) = 2J1(Rθ)/(Rθ), θ =

p

ξ

2 + η2 (18) ,

where J1(x) is the Bessel function of the first kind of order 1. A 1-D cross section

of the sombrero function when R = 0.08 is shown in Figure 9(A) and compared

with a Gaussian. Whereas in class G psfs the severity of the blur is determined

by how rapidly the optical transfer function |hˆ(θ)| decreases as θ ↑ 256, in defocus

blurs, severity is determined by the number of zeroes1in |hˆ(θ)| on 0 < θ ≤ 256.

Taking logarithms of absolute values in Figure 9(B) emphasizes the difference in the

otf signatures of the two types of blur.

Defocus otfs can be generated in Fourier space by using (18) and letting the

frequencies ξ, η be integers with −256 ≤ ξ, η ≤ 256. Specifying R uniquely determines

the defocus otf. Synthetically defocused images can be generated via FFT algorithms

1The first five positive zeroes of J1(x) are 3.83171, 7.01559, 10.17347, 13.32369, and 16.47063.
BLIND DECONVOLUTION II 11

A B

Fig. 9. Fourier domain behavior of Gaussian and defocus blurs. (A) Comparison of e−aξ2

(dashed) with 2J1(Rξ)/(Rξ) (solid), on −250 ≤ ξ ≤ 250, when a = 0.001, R = 0.08, and defocus

otf has 6 zeroes in 0 < ξ ≤ 256. (B) Comparison of corresponding logarithms of absolute values.

Zeroes in defocus otf appear as sharp minima.

A B C

Fig. 10. Deblurring with incorrect defocus parameter and honeycomb artifacts. (A) Noiseless

defocused Sydney image when R = 0.1. Blurred image computed and stored in 64-bit precision. (B)

Successful deblurring of image A using correct defocus value R = 0.1, and modified SECB procedure

with qˆ(ξ, η) as in (34), s = 0.001 and K = 500. (C) Deblurring with incorrect value R = 0.11

produces ‘honeycomb’ artifacts when K = 0.5, despite absence of noise in blurred image. With

larger errors in R and/or higher values of K, artifacts become more severe, eventually obscuring

image entirely.

and multiplication with 2J1(Rθ)/(Rθ). Deblurring such defocused images when the

defocus parameter R is known, can be done using the modified SECB procedure

discussed in the Appendix. This is a fast FFT-based direct method.

It should be noted that the defocus variable R used in this paper is a Fourier

domain variable that enters all calculations only as the argument of the Bessel function

J1 in (18). Our numerical values for R are not comparable to those used in some

studies where R denotes the radius or diameter of the defocus circle, and is measured

in pixels. In comparing the severity of the blurs discussed here with those in other

studies, the relevant measure should be the number of zeroes of ˆh(θ) on 0 < θ ≤ 256.

Our first example is the defocused Sydney image with R = 0.1 in Figure 10(A).
12 ALFRED S. CARASSO

This was calculated and stored in 64-bit precision and may be assumed noiseless.

With the otf known exactly, excellent deblurring is accomplished in Figure 10(B),

using the modified SECB procedure in (33) with s = 0.001, K = 500, and ˆq(ξ, η) as

in (34). In particular, the umbrellas and people on the Opera House deck, and the

logos on the sailboats, are nicely resolved at this high value of K. When the blurred

image contains noise, lower values of K must be used, resulting in lower resolution.

Of great significance to the present paper is the image in Figure 10(C) that results

when the same procedure is used to deblur Figure 10(A), but with the incorrect value

R = 0.11. The honeycomb artifacts obscuring the image in Figure 10(C) are not

typical noise artifacts, as the blurred image is noiseless and K was reduced to the

value K = 0.5. Rather, this phenomenon may be understood as follows. The SECB

formula in (33) may be written as

ˆf

†(ξ, η) = ˆk(ξ, η)ˆg(ξ, η),

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

|hˆ(ξ, η)|

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

(19) .

where ˆk(ξ, η) is the Fourier space representation of the regularized inverse to the

defocus blur operator. The blurred image |gˆ(ξ, η)| is zero or has very small values at

the points in the (ξ, η) plane where |hˆ(ξ, η)| is zero or has very small values. This set of

points, Λ0.1, is determined by the value R = 0.1. With that value of R, the regularized

inverse |

ˆk(ξ, η)| develops appropriate spikes on the point set Λ0.1

, so that the correct

deblurred image results upon multiplication of ˆk(ξ, η) with ˆg(ξ, η). However, with

R = 0.11, |

ˆk(ξ, η)| develops its spikes on the dislocated point set Λ0.11, where the

blurred image |gˆ(ξ, η)| need not always be small or zero. Inevitably, there is false

overamplification of power at selected frequencies in ˆg(ξ, η). This ‘resonance’ effect

produces the artifacts in Figure 10(C). With larger errors in R, and/or larger values

of K, these perturbations change character and become more intense, eventually

obscuring the image entirely.

Several techniques for blind deconvolution of defocus and motion blurs seek to

exploit the regular patterns of zero crossings that occur in the optical transfer functions describing such blurs. The main idea is illustrated in Figure 11. Using R = 0.08,

the sharp Mariner 10 image in Figure 1(B) was synthetically blurred by Fourier multiplication with (18), and the resulting image was stored in 32-bit precision. At

this low level of noise, the trace of log |gˆ

∗

(ξ, 0)| in Figure 11(A) (solid curve) displays a regular pattern of sharp minima, at the same locations as those found in

log |2J1(0.08 θ)/(0.08 θ)| (dashed curve). The defocus parameter R can easily be

determined from blurred image data alone in that case. However, when the blurred

image is stored in 8-bit precision, as in Figure 11(B), log |gˆ

∗

(ξ, 0)| is seriously contaminated by 8-bit rounding noise, and the regular pattern of sharp minima is no

longer evident. The first three local minima in Figure 11(B), at ξ = 47, ξ = 90,

and ξ = 107, respectively yield R = 0.0815, R = 0.0780, and R = 0.095. At higher

noise levels, reliable determination of R from blurred image data alone is generally

not possible.

Using additional a-priori information, a more sophisticated and computationally

intensive iterative procedure is developed in [13], [12]. The method is based on viewing the blurred image as a noisy observation of an autoregressive moving average

(ARMA) Markov random field. Blur detection and image reconstruction are reformulated as an ARMA parameter identification problem. A maximum likelihood principle

characterizes the desired parameter values as those values having most likely resulted
BLIND DECONVOLUTION II 13

A B

Fig. 11. Influence of noise on defocus parameter identification from blurred image trace. (A)

Trace of defocused Mariner 10 image with R = 0.08, computed and stored in 32-bit precision (solid

line), together with trace of defocus otf (dashed line). Sharp local minima in blurred image trace

correctly determine defocus parameter R. (B) Trace of defocused Mariner 10 image with R = 0.08,

stored in 8-bit precision. Pattern of sharp local minima, now seriously affected by 8-bit rounding

noise, does not reliably determine R. First three local minima, at ξ = 47, ξ = 90, and ξ = 107,

respectively yield R = 0.0815, R = 0.0780, and R = 0.095.

in the observed image. The expectation-maximization (E-M) algorithm is then used

to obtain the maximum likelihood solution. The procedure requires prior estimates

of image and point spread function support sizes, as well as preprocessing of the

boundaries in the blurred image. The authors stress that the E-M algorithm is only

guaranteed to converge to a stationary point, and often converges to one of several

local minima rather than to the desired global minimum. This results in erroneous

parameter estimates and a failed reconstruction. The convergence point is strongly

dependent on the initial guesses for parameter values. Repeated trials with different

initial values are usually required. Good a-priori information about the ideal image

and the blurring kernel is essential for useful results.

The method described below is a fast FFT-based direct method. At the levels of

noise and blur intensities where it produces useful reconstructions, the method may be

used in two different ways. It may be used as a stand-alone direct blind deconvolution

technique, or it may be used in conjunction with the above maximum likelihood

procedure, or some other iterative approach, to provide good initial values and other

pertinent information necessary for convergence to the desired global minimum.

5. A modified BEAK method for defocus blur. This procedure assumes

knowledge that the blur is a defocus blur and requires a substitute image. The idea

is similar to that in section 3, except that logarithms are not used. We again begin

with the basic relation

(20) g(x, y) = h(x, y) ⊗ fe(x, y) + n(x, y),

where the noise n(x, y) is assumed to satisfy

Z

R2

|n(x, y)|dxdy ≪

Z

R2

(21) fe(x, y)dxdy = σ > 0,
14 ALFRED S. CARASSO

A B

C D

Fig. 12. Blind deblurring of spacecraft image using BEAK psf detection method. (A) BEAK

detection of approximate optical transfer function in 8-bit defocused Mariner 10 image (left), using gross behavior in sharp Mariner 5 image (right). (B) Least squares fit of trace in (A) with

2|J1(Rξ)/(Rξ)| on −250 ≤ ξ ≤ 250, returns R = 0.081436, not far from exact value R = 0.08.

(C) Defocused Mariner 10 image. (D) Successful deblurring of image (C) using detected value

R = 0.081436, and modified SECB procedure with qˆ(ξ, η) as in (34), s = 0.001 and K = 0.5. Faint

honeycomb artifacts visible against dark background.

so that

|nˆ

∗

(22) (ξ, η)| ≪ 1.

With R > 0 and θ =

p

ξ

2 + η2, let Ω = {(ξ, η) | θ ≤ ω} be a neighborhood of the

origin where

2 |J1(Rθ)/(Rθ)| | ˆfe

∗

(ξ, η)| ≫ |nˆ

∗

(23) (ξ, η)|.

Such an Ω exists from (22), since the left hand side of (23) is unity at the origin. For

(ξ, η) ∈ Ω we have

|gˆ

∗

(ξ, η)| ≈ |hˆ(ξ, η)| | ˆfe

∗

(ξ, η)| = 2 |J1(Rθ)/(Rθ)| | ˆfe

∗

(24) (ξ, η)|.

Hence, for |ξ| ≤ ω,

|gˆ

∗

(ξ, 0)| ≈ 2 |J1(Rξ)/(Rξ)| | ˆfe

∗

(25) (ξ, 0)|.
BLIND DECONVOLUTION II 15

We now replace the unknown |

ˆfe

∗

(ξ, 0)| in (25) with e

vs(ξ)

, where vs(ξ) = −a|ξ|

b

is the

least squares approximation to log |

ˆfs

∗

(ξ, 0)|, and fs(x, y) is the substitute image for

fe(x, y). The detection procedure is then the following. With positive a and b given

a-priori, find a positive number R such that 2 |J1(Rξ)/(Rξ)| best fits |gˆ

∗

(ξ, 0)| e

a|ξ|

b

on |ξ| ≤ ω. Note that the expression

u(ξ) = |gˆ

∗

(ξ, 0)| e

a|ξ|

b

(26) , |ξ| ≤ ω

which involves gross behavior in a substitute image, is a very rough approximation to

|hˆ(ξ, 0)|, the absolute value of the actual defocus optical transfer function on the line

η = 0. Note also that use of logarithms would require the fitting of log{u(ξ)} with

log{2 |J1(Rξ)/(Rξ)|}. The sharp spikes in the latter function are not helpful in the

fitting algorithm. Finally, fitting u(ξ) produces quite different plots than was the case

in Figure 6(A). In particular, there is no beak at ξ = 0 in this version of the BEAK

method.

We illustrate this procedure with two examples. In Figure 12 the Mariner 10

spacecraft image in Figure 1(B) was blurred by Fourier multiplication with (18) using

R = 0.08. As noted in Figure 11(A), the otf has 6 zeroes on 0 < ξ ≤ 256. The resulting

blurred image, rounded to 8-bits, has the trace shown in Figure 11(B). With the

Mariner 5 image in Figure 1(A) as the substitute image, we have a = 2.81, b = 0.190.

We now form u(ξ) in (26). The plot of this approximate otf is shown in Figure 12(A).

Notice that the maximum value in u(ξ) is almost 5, whereas the correct value in the

true otf is 1. In Figure 12(B) we best fit u(ξ) with 2 |J1(Rξ)/(Rξ)| on |ξ| ≤ 250. This

returns R = 0.081436, not far from the exact value R = 0.08. Little change in the

returned value of R results when the fitting procedure is used on the smaller interval

|ξ| ≤ 100. The vertical scale in Figure 12(B) is restricted to a maximum of 2 so as

to better display the resulting fit (solid curve). Using the modified SECB procedure

with ˆq(ξ, η) as in (34), s = 0.001, K = 0.5, and this detected value of R, produces a

quite good reconstruction. Faint honeycomb artifacts are in fact visible against the

dark background.

Our second example is the Marilyn Monroe image in Figure 3(A), defocused

with R = 0.12. In that case the otf has 9 zeroes on 0 < ξ ≤ 256. The blurred

image was rounded to 8-bits, and this was followed by the addition of 1% uniformly distributed random noise, i.e., each 8-bit pixel value g(x, y) was replaced by

{1 + 0.01 n(x, y)}g(x, y), where n(x, y) is a random number drawn from a uniform

distribution in the range [−1, 1]. The Ingrid Bergman image in Figure 3(B) is now

used as the substitute image. With a = 3.11, b = 0.183, we form u(ξ) in (26). This

approximate otf, displayed in Figure 13(A) on |ξ| ≤ 150, has a maximum value of

about 3. A least squares fit to these data with 2 |J1(Rξ)/(Rξ)| on |ξ| ≤ 150 is shown

in Figure 13(B). This returns R = 0.118816, not far from the exact value R = 0.12.

With this detected value of R, the modified SECB procedure with ˆq(ξ, η) as in (34),

s = 0.001, and K = 0.25, results in a high quality reconstruction. Noise artifacts are

clearly visible in Figure 13(D).

Similar behavior is found in many other examples of moderately defocused images

at low levels of noise, provided a good substitute image is used. However, there are also

many cases where the BEAK method does not provide sufficiently accurate values for

R, and further analysis is required to obtain a useful restoration. The next example

points up the limitations of our procedure, and signals some essential difficulties that

attend blind deconvolution by any method.
16 ALFRED S. CARASSO

A B

C D

Fig. 13. Blind deblurring of noisy face image using BEAK psf detection method. (A) BEAK

detection of approximate optical transfer function in noisy defocused Marilyn Monroe image (left),

using gross behavior in sharp Ingrid Bergman image (right). Exact otf has 9 zeroes on 0 < ξ ≤ 256,

and 1% noise was added to the 8-bit blurred image. (B) Least squares fit of trace in (A) with

2|J1(Rξ)/(Rξ)| on −150 ≤ ξ ≤ 150, returns R = 0.118816, not far from exact value R = 0.12.

(C) Noisy defocused Marilyn Monroe image. (D) Successful deblurring of image (C) using detected

value R = 0.118816, and modified SECB procedure with qˆ(ξ, η) as in (34), s = 0.001 and K = 0.25.

Visible noise artifacts.

6. Failure in the BEAK method. The Mariner 10 experiment in Figure 12

was repeated using R = 0.25 when the defocus otf has 20 zeroes on the interval

0 < ξ ≤ 256. As shown in Figure 14(A), this is a severe blur. As before, the blurred

image was rounded to 8-bits and the Mariner 5 image was used as the substitute

image. A least squares fit of the resulting approximate otf with 2 |J1(Rξ)/(Rξ)| on

|ξ| ≤ 50 returns R = 0.22236. Apparently, this is insufficiently close to the correct

value R = 0.25. Restoration with that detected value of R produces severe honeycomb

artifacts that obscure the image, as shown in Figure 14(B).

While the presence of 8-bit noise in the blurred image is not helpful, noise is not

the primary cause of failure in this case. Indeed, if the R = 0.25 defocused image is

computed and stored in 64-bit precision, and the above procedure is reapplied to this

noiseless blurred image, the maximum detected value of R, over a number of trial

subintervals |ξ| ≤ ω ≤ 256, is R = 0.22468. Deblurring the noiseless image using that

detected value of R produces results very similar to Figure 14(B), even when K is
BLIND DECONVOLUTION II 17

A B

Fig. 14. Failure of BEAK method in severely blurred image. (A) 8-bit defocused Mariner 10

image when R = 0.25 and otf has 20 zeroes on 0 < ξ ≤ 256. BEAK detection of approximate

transfer function using sharp Mariner 5 image, and subsequent least squares fit on −50 ≤ ξ ≤ 50,

returns R = 0.22236, insufficiently close to correct value R = 0.25. (B) Deblurring of image (A)

using detected value R = 0.22236, and modified SECB procedure with qˆ(ξ, η) as in (34), s = 0.001

and K = 0.25, produces severe honeycomb artifacts that obscure image.

A B

Fig. 15. Misleading patterns in noiseless blurred image data give false reading for R. (A)

log |gˆe

∗

(ξ, 0)| in 64-bit defocused Mariner 10 image when R = 0.25 (solid), and log |hˆ(ξ, 0)| (dashed).

Despite absence of noise, sharp local minima in blurred image trace are not well-correlated with

zeroes of J1(0.25ξ). (B) Identifying first 3 sharpest minima in noiseless blurred data with first 3

positive zeroes of J1(x) implies R ≈ 0.15, grossly underestimating true value R = 0.25.

reduced to the value K = 0.1.

Some insight into the reason for this failure is provided in Figure 15. In Figure

15(A), which is analogous to Figure 11(A), the trace of log |gˆe

∗

(ξ, 0)| (solid line) in

the 64-bit defocused Mariner 10 image, is shown on the interval |ξ| ≤ 150, together

with the logarithm of the defocus otf when R = 0.25 (dashed line). Evidently, there

are many more local minima in the defocus otf in Figure 15(A), as compared to

Figure 11(A), even though the ξ-interval in Figure 15(A) has been reduced to improve

visibility. This produces a crowding effect. In addition, the defocus otf and blurred

image data are only available at discrete values of ξ. For a given R, the zeroes of

J1(Rξ) cannot always be well-captured on a preassigned discrete mesh. A zero of J1
18 ALFRED S. CARASSO

A B

Fig. 16. Analysis of failure of BEAK method in severely blurred image. (A) Comparison of

approximate otf u(ξ) (dashed curve) with exact otf ue(ξ) (solid curve) on |ξ| ≤ 30. Minima in

u(ξ) at ξ = 15 and ξ = 28, coincide with those in ue(ξ) and are compatible with R ≈ 0.25. (B)

Comparison of least squares fit to u(ξ) (dashed curve) with ue(ξ) (solid curve) on |ξ| ≤ 50. Minima

in least squares fit dislocated relative to true values. Fit returns R = 0.22236, locating first minimum

at ξ = 17 rather than at ξ = 15.

that does not fall on a mesh point translates into a small value of the defocus otf at

the nearest mesh point, if it is not missed altogether. And, if the otf value at that

mesh point is not sufficiently small, the resulting local minimum in the logarithmic

defocus trace will not be particularly sharp.

This behavior of the defocus otf at preassigned discrete values of ξ is translated

into similar behavior of the blurred image trace at these same ξ values. The blurred

image is the primary information in blind deconvolution, and the discrete mesh was

imposed when the digitized image was acquired. Weak local minima in log |hˆ(ξ, 0)|

become weak local minima in log |gˆe

∗

(ξ, 0)|. As a result, the crowded minima resulting

from zeroes of J1 no longer clearly stand out among other local minima in the highly

oscillatory blurred image trace. Hence, despite the absence of noise in the 64-bit

solid curve in Figure 15(A), there is no convincing regular pattern of sharp minima

analogous to that in Figure 11(A).

This point is more clearly made in Figure 15(B), which displays log |hˆ(ξ, 0)| and

log |gˆe

∗

(ξ, 0)| on the interval 0 ≤ ξ ≤ 70. Without the knowledge provided by the

dashed curve, the relatively sharp minima in the solid curve at ξ = 28, ξ = 47 and

ξ = 66, and respectively labeled 1, 2 and 3, might easily be construed as originating

from the first 3 positive zeroes of J1(Rξ). Such a scenario yields 3 distinct but not

seriously incompatible values for R, namely, R = 0.1368, R = 0.1493 and R = 0.1541.

In fact, very similar considerations are used to initiate iterative blind algorithms,

and any one of these 3 values of R might be a plausible initial guess in the ARMA

maximum likelihood approach in [13], [12], for example. Convergence to the desired

global minimum from such initial values would be doubtful.

Clearly, any detection procedure that examines the solid curve in Figure 15(A) and

concludes that R = 0.25, must do so on the basis of quite good a-priori information.

Very helpful but unavailable a-priori knowledge would be log |

ˆfe

∗

(ξ, 0)|, the trace in

the exact sharp Mariner 10 image, since subtracting the latter from the solid curve in

Figure 15(A) produces log |hˆ(ξ, 0)|, the dashed curve in Figure 15(A). Equivalently,
BLIND DECONVOLUTION II 19

we may exponentiate these traces and form

ue(ξ) = |gˆe

∗

(ξ, 0)|/|

ˆfe

∗

(ξ, 0)| = |hˆ (27) (ξ, 0)|.

The BEAK method seeks to emulate the process leading to (27) using available information. With the noisy data |gˆ

∗

(ξ, 0)| and the gross behavior vs(ξ) = −a|ξ|

b

in the

substitute Mariner 5 image, we form the approximate otf u(ξ) in (26). The trace of

u(ξ) on |ξ| ≤ 30 is shown in Figure 16(A) (dashed curve), and compared to the exact

otf ue(ξ) (solid curve). It is noteworthy that u(ξ) has local minima at ξ = 15 and

ξ = 28. These minima coincide with the first two minima in ue(ξ) and are compatible

with R ≈ 0.25. (The first positive zero in J1(0.25ξ) occurs at ξ = 15.32684 which is

not a mesh point). However, the least squares fit to u(ξ) with 2 |J1(Rξ)/(Rξ)| does

not coincide with ue(ξ) !! Instead, the more objective examination of the u(ξ) data

by the least squares algorithm results in the dashed curve in Figure 16(B), whose

minima are dislocated relative to those in the exact otf (solid curve). The returned

value for R is 0.22236, and the fit locates the first minimum at ξ = 17, rather than

at ξ = 15.

7. Possible recovery from failure. While the behavior in Figure 14(B) is more

common in severely defocused images, similar behavior due to inaccurately detected

R values can occur with more moderate blurs. There are three avenues that can

be explored to improve reconstruction in such cases. One should first compare the

first two or three minima in the fitted curve with the corresponding minima in u(ξ),

provided the latter minima are clearly defined. Comparing logarithms of both curves

is often helpful. If there is dislocation, new values for R corresponding to the minima

in u(ξ) can be tried. A second approach assumes the BEAK returned value RB to be

within striking distance of the exact value Re. For 512 × 512 images, 20 trial SECB

restorations, each with a different value of R, can be obtained in about a minute of

cpu time on current desktop workstations. Starting from RB, an interactive visual

search for Re is therefore quite feasible. In that search, care must be taken to use

low values for K in the SECB method, until a value of R is located where artifacts

disappear. Higher values of K can then be used to improve resolution. Finally, RB

can be used as the initial guess in an iterative blind deconvolution algorithm. In the

case of Figure 14(B), the BEAK returned value RB = 0.222 is a more useful initial

guess than are the values R ≈ 0.15 obtained by inspecting the solid curve in Figure

15(B).

8. APPENDIX. A modified SECB procedure for psfs not in G. The

SECB method applied to the convolution integral equation Hf = g in (3) presupposes

the psf h(x, y) to have a non-vanishing Fourier transform hˆ(ξ, η), so that for any fixed

s > 0, hˆs(ξ, η) can be uniquely defined [7, p. 555]. In addition, hˆs(ξ, η) is required to be

the Fourier transform of a probability density. This is the case for class G point spread

functions in (4), and more generally, for 2-D infinitely divisible probability densities,

[7]. The integral operator Hsis then defined to be the operator of convolution

with the inverse transform of hˆs(ξ, η). Solving (3) is equivalent to solving a timereversed diffusion equation, and can be implemented as a continuation problem using

the operator Hs. Such a continuation approach is an important element of the APEX

method discussed in [6]. The SECB constrained solution f

†

to the ill-posed equation

Hf = g is defined by

f

†

(x, y) = Arg min

f∈L2



k Hf − g k

2 +K−2

k Hsf − f k

2





(28) ,
20 ALFRED S. CARASSO

where the positive constants K, s are previously chosen regularization parameters, as

discussed in [2, 3, 4]. This minimization problem has a closed form solution in Fourier

space. We have, with z denoting the complex conjugate of z,

ˆf

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

|hˆ(ξ, η)|

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

(29) .

The functional on the right of (28) embodies the a-priori information that with suitable

fixed K, s > 0, the correct solution f

0

(x, y) satisfies

k Hf 0 − g k≤ ǫ, k Hsf

0 − f0

(30) k≤ Kǫ,

where ǫ > 0 is an L

2 bound for the noise in g(x, y). Rigorous error bounds for SECB

regularization are developed in [3, 4]. Applications of the SECB constraint to ill-posed

problems other than image deblurring are discussed in [5].

The above diffusion formalism is not applicable to out of focus blurs since the

continuous real-valued function ˆh(ξ, η) in (18) has infinitely many sign changes. The

operator Hsis not well-defined, and solving (18) is no longer equivalent to solving

a diffusion equation backwards in time. For deconvolution problems Hf = g where

hˆ(ξ, η) is not an infinitely divisible characteristic function, we may proceed as follows.

Choose ˆq(ξ, η) = exp{−α(ξ

2 + η2

)

β} with fixed α > 0 and fixed 0 < β ≤ 1, and let

Qs denote convolution with the inverse transform of ˆq

s

(ξ, η), 0 ≤ s ≤ 1. The correct

solution f

0

(x, y) of Hf = g again satisfies

k Hf 0 − g k≤ ǫ, k Q

s

f

0 − f0

(31) k≤ Kǫ,

with suitably chosen K, s > 0. Accordingly, we obtain a constrained solution f

† of

Hf = g by means of

f

†

(x, y) = Arg min

f∈L2



k Hf − g k

2 +K−2

k Q

s

f − f k

2





(32) ,

leading to the Fourier domain solution

ˆf

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

|hˆ(ξ, η)|

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

(33) .

The error analysis in [3] must be modified in order to be applicable in (32). For a

given fixed ˆq(ξ, η), use of (33) on a blurred image g(x, y) requires knowledge of the

regularization parameters K, s. These parameters will depend on the choice of Q

and represent a-priori information about the exact solution, some form of which is

always necessary in the solution of ill-posed problems. In practice, if the unknown

sharp image is an easily recognizable object, we may fix a value of s in the range

0.001 ≤ s ≤ 0.1, and adjust K interactively in (33) so as to achieve optimal results.

Higher levels of noise dictate smaller values of K, and vice-versa. Begining with a

small value of K, increasing K increases sharpness in the restored image, until a

threshold value is reached. Further increase in K brings out noise which eventually

obscures the image. In the defocus examples in Sections 4-6, qˆ(ξ, η) was chosen to be

qˆ(ξ, η) = e

−0.075(ξ

2+η2

)

1/2

(34) ,

with integer ξ, η lying between −256 and 256.
BLIND DECONVOLUTION II 21

REFERENCES

[1] M. Cannon, Blind deconvolution of spatially invariant image blurs with phase, IEEE Trans.

Acoust., Speech, Signal Processing, 24 (1976), pp. 58-63.

[2] A. S. Carasso, Overcoming H ̈older continuity in ill-posed continuation problems, SIAM J.

Numer. Anal., 31 (1994), pp. 1535–1557.

[3] A. S. Carasso, Error bounds in nonsmooth image deblurring, SIAM J. Math. Anal., 28 (1997),

pp. 656-668.

[4] A. S. Carasso, Linear and nonlinear image deblurring: a documented study, SIAM J. Numer.

Anal., 36 (1999), pp. 1659-1689.

[5] A. S. Carasso, Logarithmic convexity and the ‘slow evolution’ constraint in ill-posed initial

value problems, SIAM J. Math. Anal., 30 (1999), pp. 479-496.

[6] A. S. Carasso, Direct blind deconvolution, SIAM J. Appl. Math., 61 (2001), pp. 1980-2007,

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

Wiley, New York, 1971.

[8] J. D. Gaskill, Linear Systems, Fourier Transforms, and Optics, Wiley, New York, (1978).

[9] R. E. Hufnagel and N. R. Stanley, Modulation transfer function associated with image

transmission through turbulent media, J. Opt. Soc. Am., 54 (1964), pp. 52-61.

[10] A. Heckert and J. J. Filliben, DATAPLOT Reference Manual, available on World Wide

Web at http://www.itl.nist.gov/div898/software/dataplot/document.htm

[11] S. M. Jefferies and J. C. Christou, Restoration of astronomical images by iterative blind

deconvolution, Astrophysical Journal, 415 (1993), pp. 862-874.

[12] R. L. Lagendijk and J. Biemond, Iterative Identification and Restoration of Images, Kluwer

Academic Publishers, Norwell, MA, 1991.

[13] R. L. Lagendijk, A. M. Tekalp, and J. Biemond, Maximum likelihood image and blur identification: a unifying approach, J. Optical Engineering, 29 (1990), pp. 422-435.

[14] E. L. O’Neill, Introduction to Statistical Optics, Addison-Wesley, Reading, MA, (1963).

[15] T. J. Schulz, Multiframe blind deconvolution of astronomical images, J. Opt. Soc. Am. A, 10

(1993), pp. 1064-1073.