A comparative study of non-blind and blind deconvolution of ultrasound images

PDF

Download

Open document

Flip pages

Summarize

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

A COMPARATIVE STUDY OF NON-BLIND AND BLIND DECONVOLUTION

OF ULTRASOUND IMAGES

By
Na Ri Choi

A Masters Thesis Submitted to the Faculty of
The University of Southern California
In Partial Fulfillment of the Requirements for the
Masters Degree in Applied Mathematics

!"#$%&% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
ABSTRACT
The issue of restoration of ultrasound images through blind deconvolution has been one of the
chief problems in medical ultrasound imaging. This paper focuses on non-blind and blind
deconvolution methods, and compares those techniques against each other. To be specific, the
non-blind deconvolution method is based on attenuating the low-frequency or high-frequency
components and preserving the other. In low-pass and high-pass filters, the ideal filters exhibit a
sharp transition from passband to stopband, which contributes to the undesirable ringing effect. In
contrast, the Gaussian filter exhibits a smooth transition. Butterworth filters allow the freedom to
control the shape of the transition, thereby coming close to an ideal filter. The blind
deconvolution method discussed in this paper involves a iterative multichannel algorithm that is
derived from the spatially adaptive local polynomial approximation (LPA) and intersection of
confidence levels (ICI) method. The second blind deconvolution method that is discussed is the
maximum likelihood technique, which involves estimating the PSF and then utilizing the Wiener
filter.

!"#$%@% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
TABLE OF CONTENTS

ABSTRACT………………………………………………………………………………….….…2

TABLE OF CONTENTS..................................................................................................................3
INTRODUCTION............................................................................................................................4
Chapter 1 Image Deconvolution Model ......................................................................................5
Section 1.1 Noise Model ......................................................................................................6
Section 1.2 Image Enhancement in Spatial Domain.............................................................6
Section 1.2.1 Smoothing spatial filters..........................................................................7
Section 1.2.2 Sharpening spatial filters.........................................................................9
Section 1.3 Image Enhancement in Frequency Domain.....................................................11
Section 1.3.1 Smoothing frequency-domain Filters....................................................12
Section 1.3.2 Sharpening frequency-domain Filters....................................................15
Chapter 2 Image Deconvolution Techniques………….............................................................17
2.1 Inverse Filter...........................................................................................................17
2.2 Wiener Filter...........................................................................................................17
2.3 Lucy-Richardson.....................................................................................................18
2.4 Blind Deconvolution...............................................................................................19
2.4.1 LPA-ICI blind deconvolution algorithm………………………….........19
2.4.2 Maximum likelihood deconvolution algorithm………………..………21
Chapter 3 RESULTS.................................................................................................................23
3.1 Average and Gaussian filter results........................................................................23
3.2 Lucy-Richardson results.........................................................................................24
3.3 LPA-ICI blind deconvolution results......................................................................25
3.4 Maximum likelihood deconvolution algorithm......................................................27

Chapter 4 CONCLUSION........................................................................................................30
ACKNOWLEDGEMENTS………….…………………………………...………………………30
REFERENCES...............................................................................................................................31
BIBLIOGRAPHY...........................................................................................................................33

!"#$%A% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

INTRODUCTION
Medical ultrasound imaging is advantageous in today’s clinical practice because of its
noninvasive nature, mobility, accessibility, and safety. Unfortunately, the drawback to these
advantages is the reduced quality of ultrasound images compared to alternative imaging systems,
such as X-ray Computed Tomography or Magnetic Resonance Imaging. As a result, there have
been great efforts in the past few decades to improve the quality of medical ultrasound images [7].
To create a 2-D ultrasound image, a focused ultrasound beam is transmitted in a specific
direction in the investigated tissue. The obtained radio-frequency (RF) image!!!!!!!!is the result
of the convolution between the point-spread function (PSF) !!!!!! of the ultrasound probe and a
tissue-reflectivity function !!!!!!! The noise !!!!!! accounts for experimental and model
noises. The pair !!!!!!represents the indices, corresponding to axial and lateral directions,
respectively [1]. The convolution model of ultrasound image formation portrays this linear
relationship between the acoustic field and studied tissues, shown in the following equation [1]:
!!!! ! ! !! !!! !!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
The stationary convolution model does not account for the formation of a whole RF-
image because the PSF has spatial dependency from dispersive attenuation, diffraction effects,
and phase aberrations. However, the low spatial variability of the above phenomenon allows the
division of the whole RF-image into several segments, each modeled using the stationary
convolution model but with a different PSF. The “local” PSF and the corresponding reflectivity
function are unknown, and therefore it becomes a blind deconvolution problem [1].

!"#$%B% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

CHAPTER 1
IMAGE DECONVOLUTION MODE
According to the convolution theorem, a convolution of two spatial functions may be
expressed as a product of their respective Fourier transform. In other words, the deconvolution
model can be written as [6]:
! !!! !! !!! !! !!! !! !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
where !!!!!!!"#!! are the Fourier transforms of !!!!!! and !!!respectively.
The two-dimensional discrete Fourier transform can be expressed as a forward transformation [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! ! ! !!! !! !!!!!!!
!!!
!!!
!!!
!!!

where !!!!! !!!!!!!!!!!!and !!!!! !!!!!!!!!!!!are called transform variables. The
input image is !!!!!!, and!! and!!!are its row and column dimensions. The forward
transformation kernel is represented by ! !!!!!!! ! By applying the inverse transform of
! !!! !! we can recover the original image !!!!!! [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! ! ! !!! !! !!!!!!!
!!!
!!!
!!!
!!!

where m=0, 1,… , M-1, n=0, 1, … , B-1, and the inverse transformation kernel is represented by
! !!!!!!! .
We can further write the two-dimensional discrete Fourier transform of a function image !!!!!!
with size A!! by [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !
!
!"
! !!! !
!!!!!
!"
!
!
!"
!
!
!!!
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!!

!"#$%C% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

1.1 Noise Model
In ultrasound images, there is degradation from the presence of signal dependent noise
called speckle noise, which is multiplicative in nature. This is an inherent characteristic of
medical ultrasound imaging, which decreases the diagnostic value compared to other imaging
modalities. Since the ultrasound imaging noise is multiplicative and non-Gaussian, it is generally
more difficult to remove compared to additive noise [3]. The multiplicative noise model can be
written as [3]:
!!!! ! !!!! !!!!!!!
where the speckle image !!!!!! is the product of the original image !!!!!! and the non-
Gaussian noise !!!!!!!
In order to convert multiplicative noise into additive noise, a logarithmic transformation
is applied and then the noise is approximated as an additive zero mean Gaussian noise [3]:
!"!!!! !!"!!!!!!! !"!!!!!!

1.2 Image Enhancement in Spatial Domain
The spatial filter has a neighborhood, usually a square or rectangle and a predefined
operation performed directly on the image pixels in the neighborhood. It creates a new pixel at
the neighborhood center’s coordinates from the filtering operation. The response !!!!!! of
linear spatial filtering is the sum of products of the filter coefficients and corresponding image
pixels. For linear filtering of an image ! with a mask of size !!!!!, we let !! !!!! and
!! !!!!, where !! ! and !! !!!!"!! With a filter of size!!!!!!, we have [10]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! ! !!!!!!!!!! !!!! !!
!
!!!!
!
!!!!
!
The above equation must be repeated for !! !!!!!!!!!! and !!! !!!!!!!!!!!
to create a complete filtered image [10].!

!"#$%D% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

1.2.1 Smoothing Spatial Filters (Low-Pass Filters)
Smoothing spatial filters are typically used for blurring an image or and reducing noise in
an image. Spatial filters are based on their behavior in the spatial frequency domain. The low-
pass filters are spatial filters whose effect on the output image are the same as preserving low-
frequency components, such as coarser image details, and also have an effect of attenuating high-
frequency components, such as fine details. Blurring includes the removal of small details of an
image prior, object extraction, and bridging small gaps in curves. Noise reduction can be
achieved by blurring a linear filter by a nonlinear filter [6].
Mean Filter
The mean filter, also called neighborhood averaging, is the simplest and most popular
spatial smoothing filter. The response of a smoothing linear spatial filter is the average of pixels
in a neighborhood. By using smoothing filters, the value of every pixel is replaced by the average
of the gray levels in the neighborhood. The convolution mask for a !!!!! mean filter is given by
[6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
! ! !
! ! !
! ! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
In weighted average filtering, each pixel is multiplied by differing coefficients, so more
importance is given to some pixels at the expense of other pixels. An example of a weighted
average mask is given by [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !
!
!"
!
!
!
!"
!
!
!
!
!
!
!
!"
!
!
!
!"
!
!
!"
! ! !
! ! !
! ! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
The general formula to use a weighted averaging filter on an M!N image is given by [6]:
!"#$%E% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
! !!! !
! !!! !!!! !!!! !!
!
!!!!
!
!!!!
! !!!
!
!!!!
!
!!!!
!
Gaussian Filter
The Gaussian filter uses a discrete kernel from a radially symmetric form of continuous Gaussian
function [10]:
! !!! !
!
!!!
!
!
!
!!
!
!!
!
!!
!
!

There are two free parameter’s when specifying discrete approximations to this function [10]:
1. The favorable kernel size (!!! filter mask);
2. The standard deviation !of Gaussian function.
The Gaussian filter has a smoothing effect on an image, and the degree of that effect doesn’t
depend on the absolute value of the kernel size, but instead depends on the standard deviation.
Another nice property of the Gaussian function is that the Fourier transform is a Gaussian
function as well. A large standard deviation value of a Gaussian function is an example of a low-
pass filter where the high spatial frequency feature of an image is minimized. The Gaussian filter
removes some of the existing noise and degrades high frequency detail [10].
Order-Statistics Filters
Order-statistics filters, otherwise known as nonlinear spatial filters, have responses based
on ordering or ranking pixels under the kernel, rather than processing pixel values using the
convolution operator. A popular filter under this category is the median filter. The median filter is
when each pixel value is replaced by the median of gray levels in the neighborhood. Median
filters are very good for filtering impulse noise, or very bright (salt) and very dark (pepper) noise
in an image. In general, median filters are far superior in removing salt-and-pepper noise,
compared to average filtering [6].
Median Filtering
The median filter is a commonly used. The median filter is superior to the mean filter because the
median filter fufills the limitations of the mean filter. To be specific, the median filter is better at
keeping the sharp high-frequency features and removing noise, especially salt and pepper noise.
!"#$%F% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
Each pixel is replaced by a statistical median in the !!! neighborhood, instead of being
replaced by the mean [10].
The median m is the number where it falls in the middle of a set of numbers; it can be
viewed as a midpoint of values. The median is a pixel value that comes from its pixel
neighborhood, so it has strength against outliers and does not produce a new impractical pixel
value. This avoids adge blurring and loss of image fine features. The ordering of values in the
neighbhoordhood at each pixel location is computationally lengthy [10].

1.2.2 Sharpening Spatial Filters (High-Pass Filters)
The main idea behind sharpening spatial filters is to emphasize fine detail in an image or
enhance a detail that was blurred. The high-pass filters are spatial filters whose effect on the
output image is the same as preserving high-frequency components, such as fine details and edges
in an image, and attenuating low-frequency components. In general, first-order derivatives create
thicker edges in an image, whereas second-order derivatives have a stronger response to fine
detail, including isolated points and thin lines. Also, first-order derivatives have a stronger
response to a gray-level step function, whereas second-order derivatives create a double response
at step changes in gray level. In most real applications, the second-order derivative is superior to
the first-order derivative because the former has the ability to enhance fine detail [6].
One use of second-order derivatives is the Laplacian. The Laplacian for an image
!!!!!! is written as [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!! !
!
!
!!!
!!
!
!
!
!
!!!
!!
!

The second derivatives are typically approximated as [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!!!
!!
!
! ! !!!!! !! !!!!! !!!!!!!!
and
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!!!
!!
!
! ! !!!!! !! !!!!! !!!!!!!!
We can then express the Laplacian as a sum of products:
!"#$%GH% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
!!!!!!!!!!
!
!!! ! ! !!!!! !! !!!!! !! !!!!! !! !!!!! !!!!!!!!
This equation can be implemented using the convolution mask, giving an isotropic result for
rotations in !"
!
increments [3]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! !! !!!!
!! !!!! !!
!!!! !! !!!!

An alternate implementation of the Laplacian includes the diagonal neighbors [3]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !! !!
!! !!!! !!
!! !! !!

It should be noted that other common implementations include the reversed signs of each
coefficient in the convolution masks above.
1.3 Image Enhancement in Frequency Domain
The formal discrete convolution of two functions !!!!!! and !!!!!! of size !!! is defined as
the following [10]:
!!!! !!!!! !
!
!"
!
!!!
!!!
!!!
!!!
!!! !!!!!!!!!!
Letting !!!!!! and !!!!!! represent the Fourier transforms of !!!!!! and !!!!!!, respectively,
the convolution theorem states that !!!! !!! and ! !!! !! !!! form a Fourier transform
form. Thus, the result is the following [10]:
!!!! !!!!! !! !!! !! !!! !
Filters in the spatial domains and frequency domains form a Fourier transform pair. In
other words, given a filter in a frequency domain, the corresponding filter in the spatial domain is
obtained by taking the inverse Fourier transform of the former, and the reverse also holds true. In
practice, taking the inverse transform of a filter in the frequency domain to compute an equivalent
spatial domain filter of the same size is not computationally efficient. Instead, it is more
computationally efficient to do filtering in the frequency domain, but we use smaller filters in the
spatial domain. If possible, it is best to specify filters in the frequency domain, take their inverse
!"#$%GG% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
transform, and use the resulting filter in the spatial domain to create smaller spatial filter masks
[6].
Filters based on Gaussian functions are important because their shapes are easily
specified since the forward and inverse Fourier transforms of a Gaussian function produce real
Gaussian functions. To simply the notation, we will discuss with only one variable. Let H[u]
denote a frequency domain. The Gaussian filter function is given as [6]:
! ! !!!
!!
!
!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
where ! denotes the standard deviation of the Gaussian curve.
The corresponding filter in the spatial domain is calculated and written as:
!! ! !!!"!
!!!
!
!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equations (1.5) and (1.6) are important results because they form a Fourier transform pair, where
both components are Gaussian and real. Gaussian curves are intuitive and can be easily
manipulated. Another note is that these functions are reciprocal to one another. The narrower the
frequency domain filter is, the more it will reduce the low frequencies, which will increase
blurring. In spatial domain, this creates a wider filter, which means a larger mask. These
properties help create an understanding of the filtering properties in spatial and frequency
domains. The highpass filter can be constructed as a difference of Gaussian filters [6]:
!
!
! !!!
!!
!
!!!
!
!
!!!
!!
!
!!!
!
!

where !!! and !
!
! !
!
. The corresponding filter in the spatial domain is calculated and given
as follows [6]:!
!
!
! ! !!!
!
!!
!!!
!
!
!
!
!
!
! !!!
!
!!
!!!
!
!
!
!
!
!

!"#$%G&% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
1.3.1 Smoothing Frequency-Domain Filters (Low-Pass Filters)
Low-pass filters (LPF) attenuate the high-frequency components of the Fourier transform
of an image, while preserving the low-frequency components. The basic model for frequency
domain filtering is given by [10]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !! !!! !!!!!!!
where !!!!!! is the Fourier transform of an image that needs to be smoothed. The goal is
choose a filter transfer function H(u, v) that produces G(u, v) by decreasing the high-frequency
components of F(u, v). Three types of lowpass filters are discussed [10].
Figure 1. Types of filters in frequency domains. The perspective plot of transfer function, image
filter, and radial cross sections of filters for (a) ideal lowpass filter (b) Butterworth lowpass
filter, and (c) Gaussian lowpass filter are shown [2].
Type of Filters
in Frequency-
Domain
Perspective plot of
transfer function

Filter as an image Filter radial cross
sections

Ideal Lowpass
Filter

Butterworth
Lowpass Filter

Gaussian
Lowpass Filter

!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
!"#"$%&'()%#*'+,-.*//"0#1'203'*34
www.imageprocessingbook.com
©2002 R. C. Gonzalez & R. E. Woods
Chapter 4
Image Enhancement in the
Frequency Domain
Chapter 4
Image Enhancement in the
Frequency Domain
!"#$%G@% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
Ideal LPF
The first filter is the 2-D ideal lowpass filter. This is a filter that eliminates all high-frequency
components of the Fourier transform that have a distance greater than the designated distance !
!

from the origin of the transform. The ideal lowpass filter has the transfer function [10]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!! !
! !"!!!!!!!!!
!
! !"!! !!! !!
!

where !
!
is a fixed nonnegative quantity, and !!!!!!! !
!
!!
!
! When the value of !
!
is
decreased, the outcome image is not only the blurry versions of the input image, but also displays
ringing effects due to the sharp transition of stop-band and pass-band in ideal LPF [6]. All
frequencies within the circle of radius !
!
are passed without attenuation. However, all
frequencies outside of the circle are attenuated. See Figure 1(a).
Butterworth LPF
The second filter is the Butterworth lowpass filter. It is defined as [10]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!! !
!
!!
! !!!
!
!
!!

where ! !!! ! !
!
!!
!
and has a cutoff frequency !
!
with order n. The Butterworth lowpass
filter doesn’t have a sharp discontinuity like the ideal lowpass filter, which means the Butterworth
establishes a clear cutoff between filtered and passed frequencies. The shape of the frequency
response is controlled by n. When the value of n increases, the filter exhibits steeper transitions
and approaches the ideal filter [10]. As seen in Figure 1(b), the Butterworth transfer function does
not exhibit a sharp transition like the ideal LPF does. Instead, the Butterworth LPF exhibits a
smooth transition between high and low frequencies.
Gaussian LPF
The third filter is called the 2-D Gaussian lowpass filter, which is defined as [6]:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!! ! !
!!
!
!!!!!!!!
!

!"#$%GA% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
Using a transfer function based on a Gaussian curve, the Gaussian low-pass filter attenuates high
frequencies. The bell-shape curve width is controlled using a specific parameter !!!which is
equivalent to the cut-off frequency !
!
! Decreasing the value of !
!
corresponds to stricter
filtering, which results in blurrier results. The Gaussian lowpass filter does not achieve as much
smoothing compared to Butterworth lowpass filter, which is due to Butterworth’s “tight” control
of the transition between low and high frequencies about cutoff frequency, as seen in Figure 1(c)
[6].

1.3.2 Sharpening Frequency-Domain Filters
Sharpening frequency domain filters is achieved by highpass filtering, which involves the
attenuation of low-frequency components without disrupting the high frequency information in
Fourier transform. The highpass filters and lowpass filters can be described in the following
relationship [2]:
!
!!
!!! ! !!!
!"
!!!!!
where !
!!
is the transfer function of the corresponding highpass filters and !
!"
!!!!! is the
transfer function of corresponding lowpass filter. Three types of high-pass filters are discussed [6].
Ideal HPF
The ideal highpass filter (HPF) is characterized as:
! !!! !
!! !"!!!!!!!!!
!
!! !"!!"#$%!

where !
!
is the specified nonnegative quantity of the cutoff frequency and ! !!! ! !
!
!!
!
!!
See Figure 3(a) for a visual representation. An ideal HPF attenuates all frequency components
within a specified radius from the Fourier transform center, while enhancing all other components.
Similar to the ideal low pass filter, the ideal highpass filter is also unreliable. They have the same
ringing properties that cause distorted and thickened boundaries of the original object [2].

!"#$%GB% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

Figure 2. Types of high-pass filters in frequency domain: (a) 2D-frequency response plot and (b)
3D-frequency response plot [3].

Butterworth HPF
The Butterworth highpass filter transfer function with order n of the filter and cutoff frequency
!
!
from the origin is defined as:
! !!! !
!
!!
!
!
! !!!
!!

The performance of the ideal highpass filter and Butterworth highpass filter is comparable
because the center spot sizes of the both filters are similar (See Figure 3). There is a smoother
transition into higher values of cutoff frequencies with the Butterworth highpass filter. [6]

Types of
High-Pass
Filters in
Frequency
Domain
3D frequency
response plot
2D frequency
response plot

Ideal
Highpass
Filter

Butterworth
Highpass
Filter

Gaussian
Highpass
Filter

250 FREQUENCY-DOMAINFILTERING
FIGURE 11.15 Frequency response plot for an ideal HPF: (a) 3D view; (b) 2D view from
the top.
11.4.2 GaussianHPF
AGaussianhigh-passfilterattenuateslowfrequenciesusingatransferfunctionwhose
shape is based on a Gaussian curve. The behavior of the filter can be controlled by
specifying the parameter σ, which is functionally equivalent to the cutoff frequency
D
0
. The Gaussian HPF can be mathematically described as
H
G
(u,v)= 1−e
−(D(u,v)
2
)/2σ
2
(11.24)
Figure 11.16 shows the frequency response plot for a Gaussian HPF.
11.4.3 ButterworthHPF
The Butterworth high-pass filter can be mathematically described as
H
B
(u,v)=
1
1+[D
0
/D(u,v)]
2n
(11.25)
whereD
0
is the cutoff frequency andn is the order of the filter.
FIGURE11.16 FrequencyresponseplotforaGaussianHPF:(a)3Dview;(b)2Dviewfrom
the top.
250 FREQUENCY-DOMAINFILTERING
FIGURE 11.15 Frequency response plot for an ideal HPF: (a) 3D view; (b) 2D view from
the top.
11.4.2 GaussianHPF
AGaussianhigh-passfilterattenuateslowfrequenciesusingatransferfunctionwhose
shape is based on a Gaussian curve. The behavior of the filter can be controlled by
specifying the parameter σ, which is functionally equivalent to the cutoff frequency
D
0
. The Gaussian HPF can be mathematically described as
H
G
(u,v)= 1−e
−(D(u,v)
2
)/2σ
2
(11.24)
Figure 11.16 shows the frequency response plot for a Gaussian HPF.
11.4.3 ButterworthHPF
The Butterworth high-pass filter can be mathematically described as
H
B
(u,v)=
1
1+[D
0
/D(u,v)]
2n
(11.25)
whereD
0
is the cutoff frequency andn is the order of the filter.
FIGURE11.16 FrequencyresponseplotforaGaussianHPF:(a)3Dview;(b)2Dviewfrom
the top.
HIGH-PASSFILTERING 251
FIGURE11.17 FrequencyresponseplotforaButterworthHPFofordern=4:(a)3Dview;
(b) 2D view from the top.
Figure 11.17 shows the frequency response plot for a Butterworth HPF of order
n = 4.
11.4.4 High-FrequencyEmphasis
TheapplicationofaHPFtoanimageusuallycausesacrispeningofthehigh-frequency
contentsoftheimageattheexpenseofthelow-frequencycontents,whichareseverely
attenuated, leading to a loss of information present in large patches of the original
image (its low-frequency components). High-frequency emphasis is a technique that
preserves the low-frequency contents of the input image (while enhancing its high-
frequencycomponents)bymultiplyingthehigh-passfilterfunctionbyaconstantand
adding an offset to the result, that is,
H
hfe
(u,v)= a+bH(u,v) (11.26)
whereH(u,v) is the HPF transfer function,a≥ 0 andb>a. Figure 11.18 shows an
example of high-frequency emphasis.
FIGURE11.18 High-frequency emphasis: (a) input image; (b) result of applying a second-
orderButterworthHPF(withD0 =30)totheinputimage;(c)resultofhigh-frequencyemphasis
witha= 0.5 andb= 1.
HIGH-PASSFILTERING 251
FIGURE11.17 FrequencyresponseplotforaButterworthHPFofordern=4:(a)3Dview;
(b) 2D view from the top.
Figure 11.17 shows the frequency response plot for a Butterworth HPF of order
n = 4.
11.4.4 High-FrequencyEmphasis
TheapplicationofaHPFtoanimageusuallycausesacrispeningofthehigh-frequency
contentsoftheimageattheexpenseofthelow-frequencycontents,whichareseverely
attenuated, leading to a loss of information present in large patches of the original
image (its low-frequency components). High-frequency emphasis is a technique that
preserves the low-frequency contents of the input image (while enhancing its high-
frequencycomponents)bymultiplyingthehigh-passfilterfunctionbyaconstantand
adding an offset to the result, that is,
Hhfe(u,v)= a+bH(u,v) (11.26)
whereH(u,v) is the HPF transfer function,a≥ 0 andb>a. Figure 11.18 shows an
example of high-frequency emphasis.
FIGURE11.18 High-frequency emphasis: (a) input image; (b) result of applying a second-
orderButterworthHPF(withD0 =30)totheinputimage;(c)resultofhigh-frequencyemphasis
witha= 0.5 andb= 1.
250 FREQUENCY-DOMAINFILTERING
FIGURE 11.15 Frequency response plot for an ideal HPF: (a) 3D view; (b) 2D view from
the top.
11.4.2 GaussianHPF
AGaussianhigh-passfilterattenuateslowfrequenciesusingatransferfunctionwhose
shape is based on a Gaussian curve. The behavior of the filter can be controlled by
specifying the parameter σ, which is functionally equivalent to the cutoff frequency
D0. The Gaussian HPF can be mathematically described as
HG(u,v)= 1−e
−(D(u,v)
2
)/2σ
2
(11.24)
Figure 11.16 shows the frequency response plot for a Gaussian HPF.
11.4.3 ButterworthHPF
The Butterworth high-pass filter can be mathematically described as
HB(u,v)=
1
1+[D0/D(u,v)]
2n
(11.25)
whereD0 is the cutoff frequency andn is the order of the filter.
FIGURE11.16 FrequencyresponseplotforaGaussianHPF:(a)3Dview;(b)2Dviewfrom
the top.
250 FREQUENCY-DOMAINFILTERING
FIGURE 11.15 Frequency response plot for an ideal HPF: (a) 3D view; (b) 2D view from
the top.
11.4.2 GaussianHPF
AGaussianhigh-passfilterattenuateslowfrequenciesusingatransferfunctionwhose
shape is based on a Gaussian curve. The behavior of the filter can be controlled by
specifying the parameter σ, which is functionally equivalent to the cutoff frequency
D0. The Gaussian HPF can be mathematically described as
HG(u,v)= 1−e
−(D(u,v)
2
)/2σ
2
(11.24)
Figure 11.16 shows the frequency response plot for a Gaussian HPF.
11.4.3 ButterworthHPF
The Butterworth high-pass filter can be mathematically described as
HB(u,v)=
1
1+[D0/D(u,v)]
2n
(11.25)
whereD0 is the cutoff frequency andn is the order of the filter.
FIGURE11.16 FrequencyresponseplotforaGaussianHPF:(a)3Dview;(b)2Dviewfrom
the top.
!"#$%GC% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
Gaussian HPF
The Gaussian high-pass filter uses a transfer function based on a Gaussian curve to attenuate low
frequencies. The behavior of the Gaussian filter is controlled using a specific parameter !!!which
is equivalent to the cutoff frequency !
!
! The Gaussian highpass filter transfer function in 2-D
with cutoff frequency locus at distance !
!
from the origin is defined as [3]:
! !!! ! !!!
!!
!
!!!!!!!!
!
!

In Figure 2, the ideal highpass filter cuts off all components of the Fourier transform that are
within the distance !
!
from the center. However, the Butterworth and Gaussian highpass filters
do not have sharp discontinuities, and there are no clear cut-off components between the passed
and filtered frequencies. Gonzales found that the performance of the Gaussian highpass filter
achieves smoother results compared to the aforementioned ideal highpass filter and Butterworth
highpass filter [3].

2 Image Deconvolution Techniques
The implementation of image deconvolution techniques will be discussed in this section.
2.1 Inverse Filter
The inverse filter is the inverse of the degradation function. However, upon implantation,
the inverse filter will enlarge the high frequency noise and thereby with the high-pass filter nature,
and it will cause the degraded image to have a magnified high frequency noise. Therefore, the
inverse filter is defined as follows:
! !!! !
!
! !!!
!!!!!! !!! !!!!
!!!!!!!!!!!!!!!! !!! !!

where ! is a threshold use to lessen the effect of zeroes in the degradation function [7], [10].
!"#$%GD% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
2.2 Wiener Filter
With the presence of noise distribution, a more sophisticated approach than the inverse
filter is required. The Wiener filter attempts to model the error through statistical methods,
specifically through minimum mean square estimator. When the error is modeled, the average
error is minimized. If we have the original image !!!!!!, degraded version !!!!!!, and
restored version !!!!!!,then we can compute the sum of squared differences in every pixel
between the original image and corresponding restored image. Then this shows us how well the
restoration algorithm performed [10].
In Fourier domain, the Wiener filter is described as [10]:
!!!!!!!
!
!
!!!!!
!!!!
!
!!
!"
!!!!!

where !
!!
!!!!!!!
!!
!!!!!! is the noise-to-signal ratio. The !
!!
!!!!! denotes the additive noise
power spectrum and !
!!
!!!!! denotes the image power spectrum. In a high frequency region, the
noise-to-signal ratio will be relatively large, so the high frequency response of the deconvolution
filter is suppressed.

If there is an absence of noise, i.e. !
!!
!!! ! !, the Wiener filter is simply the inverse filter [3]:
! !!! !
!
!!
!!! !!
!
!
! !!!
!!!!! !!! ! !!
!!!!!!!!!!!!!!!! !!! ! !

2.3 Lucy-Richardson algorithm
The Lucy Richardson filter is an iterative nonlinear image restoration process. It attempts
to determine the maximum-likelihood solution, given the PSF and Poisson noise assumption.
Consider the discrete form of Equation (1), except that the PSF is a 2-D matrix and the input and
output distributions are vectors. This means that within the input distribution, each ith pixel has
!"#$%GE% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
value!
!
! The observed value of the ith pixel in output !
!
can be written as the following
summation [10]:
!
!
! !
!"
!
!
!

where !
!"
is the observed output pixel. The discrete PSF is usually normalized so that !
!" ! !
!
!!
The iterative Lucy-Richardson algorithm is defined as [10]:
!
!
! !
!
!
!"
!
!
!
!"
!
! !
!

It is difficult to determine the specific number of iterations needed, but it usually depends on the
complexity and size of PSF matrix. It should be noted that increasing the number of iterations not
only slows down the computational process, but also increases noise and creates a ringing effect.
Thus, the number of iterations is found manually depending on the PSF size [9].

2.4 Blind Deconvolution
2.4.1 LPA-ICI blind deconvolution algorithm
The multiframe blind deconvolution of noisy images is based on energy criterion
minimization constructed in the frequency domain and using the recursive gradient-projection
algorithm. The local polynomial approximation (LPA) of image and blur operators is used for
regularization and filtering. The intersection of confidence levels (ICI) is used for the varying
window sizes of the LPA. The LPA-ICI blind deconvolution algorithm is spatially adaptive and
nonlinear with respect to the image and blur operators [5].
The adaptive LPA-ICI filtering algorithm forms directional linear filters with rotated
directional non-symmetric kernel !
!!!
!produced by the LPA. The angle ! represents the filter
!"#$%GF% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
directionality, which is defined by non-symmetric window function used in LPA. The kernel
support length is represented by !, or the kernel scale parameter in the direction. The directional
estimates are computed using the following convolution [5]:
! !
!!!
!
! !! !!
! !!!
!
! !
In the frequency domain, the directional estimates can be computed using the
corresponding Fourier transform of the convolution [5].
! !
!!!
!
!! !
!!!
!
!!
!
!!!
For every estimation pixel !, the ICI algorithm selects an adaptive scale parameter !.
Then the estimates ! !
!!!
!
are computed for a grid of !!"! !!
!
!!
!
!! !!
!
! where!
!
! !
!
!
!
!
!!! !
!
!!The adaptive scale is represented by !
!
, which is determined by the largest scales
in H, where the estimate is similar to the estimates corresponding to the smaller window sizes.
First, we take a sequence of confidence intervals [5]:
!
!
! ! !
!!!
!
!!!
! !!
!
!
!
! !! !
!!!
!
!!!
! !!
!
!
! ! !!! !!!!! !!
where !
! !!
!
!
!
! ! !
!
!
!!
!
!
!!!
!
represents the standard deviation of ! !!!
!
!
!
and !! ! is a
parameter. Then, the algorithm considers the confidence intervals intersection !
!
! !
!
!
!!!
and
let!!
!
represent the largest indices s where !
!
is non-empty. The optimal scale !
!
is expressed as
!
!
! !
!
! and thus the optimal scale estimate is ! !!!
!
!
!!
!
for each direction ! [5].
The fixed !!parameter plays an important part of the algorithm because it reveals whether
the estimate deviations are small or large. If the estimate deviation is too small, then it leads to
signal undersmoothing. If the estimate deviation is large, it leads to oversmoothing. For
convenience, we write the LPA-ICI multidirectional algorithm as an adaptive filter with inputs !
and ! and an output !. It can be written as !! !! !!! where !! is an operator function [5].
The LPA-ICI algorithm is done in the spatial domain and requires both forward and
inverse Fourier transforms of the frequency domain estimates !
!
!!!
!!
!!!
. The projection
!
!
!
!
{!
!
!
} indicates that the PSF !
!
!
maximal size does not exceed a fixed value. This PSF
size reduction also reduces the amount of computations [5]. Katkonick, et al. describes the blind
image multichannel devolution algorithm as follows [5]:
!"#$%&H% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
1. Initialization: Gaussian density is used for the initial estimate !
!
!!!
and the mean of
observed images for initial !
!!!
! !
!
!!!!!
!
!!!
.
2. Image estimation: Compute !
!!!
without projection, where
!
!!!
!!
!
!
!!!
!
!
!!!
!!
!
!
!
!
!
! !!!
!
!
! !
!
!
!!!
!
!
!
!
!
!
!
!
!
3. Image filtering: Filter !
!
using the LPA-ICI algorithm.
3a. Compute the inverse of the Fourier transform: !
!!!
!!
!!
!!
!!!
!.
3b. Compute the standard deviation estimate !
!
!!! of the noise !
!!!
by using the
signal’s differences estimates.
3c. Filter !
!!!
using the algorithm !
!!!
! !! !
!
!!
!
!
4. Image projection: Project the image !
!!!

4a. Project !
!!!
onto the interval segment !!! !!
!!!
!!
!
!
!
!
,
where !
!
!
!!"# !!!"#!!!!!! .
4b. Compute !
!!!
!!!!
!
!.
5. PSF estimation: Compute !
!
!!!
without filtering and projection, repeating this process
!
!"#
many times. The reason for this is to speed up the convergence rate.
6. PSF projection: Project the PSF !
!

6a. Compute !
!
!!!
!!
!!
!
!
!
!
6b. Project the PSF estimates !
!
!!!
!!
!
!
!
!
!
!
,
where !
!
!
!
!
!
!
!
!
!
!
! !
!
!!!!!!
!
! ! and !
!
! ! !!!!!"!!! !
!
! !
!
!!.
7. PSF filtering: Filter !
!
!
using the LPA-ICI algorithm.
7a. Compute the noise standard deviation in !
!
!!!
!
7b. Filter !
!
!!!
using !
!
!!!
! !! !
!
!
!!
!
!
!
8. Increase the value of of k and repeat Step 2 – Step 8 K amount of times.

!"#$%&G% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

2.4.2 Maximum Likelihood Blind Deconvolution
Consider the original deconvolution problem, previously discussed in Equation (1.0), defined on
rectangular coordinates by !!
!
!!
!
! (!! !!!!!!! !!! !!) [8]:
!!
!
!!
!
! !!
!
!!
!
!!!
!
!!
!
! !!!
!!
!!
!!
!!!
!
!!
!
!
!
!!!!"
!!
where we wish to estimate the original image !! Since we do not know the noise !!and the point
spread function !, the problem looks impossible to deconvolve.
However, we can use feasible solutions through iterative procedures, done within the frequency
domain. There are two basic constraints in feasible solution [10]:
1. They have finite support. In other words, the desired input distribution is limited to a
specific spatial area and is zero outside of the said area.
2. Any potential solution is strictly positive.
The original image !!!
!
!!
!
! is assumed to be defined by an autoregressive process of low order
[8]:
!!
!
!!
!
! !!
!
!!
!
!!!
!
!!
!
! !!!
!!
!!
!!
!!!
!
!!
!
!
!
!!!!"
!!
where !!
!!
!!
!!
is the image model coefficient. Next, let us calculate the Fourier transforms of
the above two equations to define the likelihood function [8].
! !
!
!!
!
!! !
!
!!
!
! !
!
!!
!
!!!!
!
!!
!
!
! !
!
!!
!
!! !
!
!!
!
! !
!
!!
!
!!!!
!
!!
!
!
where the Fourier transforms correspond to the lowercase letters.
The Fourier phase is obtained by the Fourier modulus and the distorted image h by using an
iterative algorithm. The maximum likelihood estimation aims to identify the autocorrelation, not
necessarily the PSF itself. The unknown parameters are [8]:
!! !! !
!
!!
!
!! !
!
!!
!
!!
!
!
!!
!
!
!
where r !
!
!!
!
represents the autocorrelation function. It is defined as [8]:
! !
!
!!
!
! !
!
!
!!!!"
!!
!!
!!
!!
!!!
!
!!
!
!!!
!
!!
!
!
!"#$%&&% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
The maximum likelihood estimator of the parameter vector ! is characterized by:
!
!"
!{!"#
!
!!!!}
where !!!! represents the log-likelihood function of !.
This blind deconvolution algorithm attempts to restore the PSF from the degradation, not only the
original input distribution. The main disadvantage of the Lucy-Richardson algorithm is that it has
noise amplification issues, so if the image is too noisy, then the algorithm will perform too many
iterations, which would produce a speckled restored image [8]. In Figure 11, this deconvolution
procedure exhibited some ringing effects, but the restored images are clearly sharper.

!"#$%&@% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
3 Results
3.1. Average and Gaussian Filter Results
The following results were obtained through the MATLAB program, with the code provided by
Solomon and Breckdon [10]. The real ultrasound images were provided by Ultrasonix. Figure 4
provides the result of an average and Gaussian filter applied to an ultrasound image. Although
there are no drastic differences between the two filters, the Average filter did better than the
Gaussian filter. In Figure 4, there is a clearer facial definition, whereas the Gaussian filter caused
the image to oversmooth.

Original
Average

Gaussian

Figure 4. The above image is the ultrasound of twins in 8 week gestation. (a) Original image (b)
Average filter (c) Gaussian filter.
!"#$%&A% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

3.2 Lucy-Richardson Results
The following results were obtained through the MATLAB program, with the code provided by
Solomon and Breckdon [10]. The real ultrasound images were provided by Ultrasonix. After
specifying the PSF and variance of noise, the first observation is produced by 10 iterations of LR
deconvolution. The second observation is an estimate of the overall standard deviation of the
noise. The third observation is when a weighing function, which gives zero weight to the pixels
that are on the image’s border, is applied to eliminate ringing effects in the previous observations.
(See Figure 5). When comparing Lucy-Richardson to the other non-blind deconvolution
counterparts, the Lucy-Richardson did a superior job deconvoluting the distorted ultrasound
image. This can be confirmed in the error analysis comparison table in Figure 6.

Original

1
st
observation

2
nd
observation

3
rd
observation
Figure 5. Results from Lucy-Richardson algorithm.
!"#$%&B% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
Filters MSE PSNR RMSE SNR
Average 17.9392 35.9287 4.2355 30.3427
Gaussian 21.0863 33.9732 4.5920 30.3834
Lucy-Richardson 3.2589 31.2592 1.8052 29.8526
Figure 6. Error values (a) Mean Square Error (b) Peak Signal to Noise Ratio (c) Root Mean
Square Error and (d) Signal to Noise Ratio, for average, Gaussian, and Lucy-Richardson filters

3.3 LPA-ICI blind deconvolution results

The gray scale image that is used for this algorithm is a second trimester fetus. The image source
is Ultrasonix. Through a code derived from LASIP [6], the MATLAB implementation produced
the following results. Figure 8 displays the multiple observations of frames in the image, while
Figure 9 portrays the corresponding PSF. Figure 10 displays the error terms. The numerical
results shows that the restoration quality is good, and shows that this LPA-ICI blind
deconvolution technique is efficient for noisy data in ultrasound images.

Original Image

Deconvoluted Image
Figure 7. Blind deconvolution through Local polynomial approximation (LPA) and intersection
of confidence levels (ICI) algorithm.

!"#$%&C% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

Figure 8. LPA-ICI iterative steps

Figure 9: LPA-ICI blind deconvolution point-spread function estimation corresponding to each
observation. Column (a) PSF for Observation 1 (b) PSF for Observation 2 (c) PSF for
Observation 3

!"#$%&D% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

Estimate Errors Fetus Image
ISNR (Improved Signal to Noise Ratio) 10.9996
SNR (Signal to Noise Ratio) 25.7986
PSNR (Peak Signal to Noise Ratio) 33.4755
MSE (Mean Square Error) 29.2101
RMSE (Root Mean Square Error) 5.4046
MAE (Mean Absolute Error) 4.0122
MAX (Maximum Absolute Error) 34.6602
Figure 10. LPA-ICI blind deconvolution estimate errors computation according to Figure 7.

3.4 Maximum Likelihood Blind Deconvolution results
The maximum likelihood blind deconvolution algorithm can be used without knowledge of the
distortion. The algorithm provides a restoration of the image and the PSF at the same time. A
damped, accelerated, and iterative Lucy-Richardson algorithm restores the image. The Lucy-
Richardson algorithm was originally created for data with discrete events that have a Poisson
distribution. The maximum likelihood maximizes the likelihood that the image, when it is
convolved with the PSF, is an instance of the distorted image, with the assumption of Poisson
noise [4]. The image that is used for this algorithm is a 12 week fetus, with its hand visible.
Through a code derived from Solomon and Breckon [10], the MATLAB implementation
produced the following results. In figure 12(b), the image was produced by smoothing the edges
of the image and performing an initial estimate of the PSF. In figure 12(c), the image was
produced by giving another estimate of the PSF. In figure 12(d), the image was produced by
performing a Wiener filter with the blindly recovered PSF [10]. Figure 13 shows the error values,
and they demonstrate a better restoration quality that the LPA-ICI algorithm.
!"#$%&E% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
Original Image
Deconvoluted Image
Figure 11. Blind deconvolution through maximum likelihood.

Figure 12. Process of maximum likelihood blind deconvolution. Top row: (a) Original Image (b)
First deconvolved image. Bottom row: (c) Deconvolve using initial estimate PSF (d) Wiener filter
with blind recovered PSF.

!"#$%&F% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

Figure 13. Maximum likelihood blind deconvolution estimate errors computation according to
Figure 11.

Image Comparison Between LPA-ICI and Maximum Likelihood Blind Deconvolution

Original Image

LPA-ICI
blind deconvolution

Maximum likelihood
blind deconvolution
Figure 14. (a) Original image (b) Result of blind deconvolution using LPA-ICI algorithm (c)
Result of blind deconvolution using maximum likelihood estimates.

Estimate Errors Fetus Image
ISNR (Improved Signal to Noise Ratio) 9.8642
SNR (Signal to Noise Ratio) 26.0527
PSNR (Peak Signal to Noise Ratio) 34.7321
MSE (Mean Square Error) 28.9375
RMSE (Root Mean Square Error) 5.3531
MAE (Mean Absolute Error) 5.9852
MAX (Maximum Absolute Error) 40.7325
!"#$%@H% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
CONCLUSIONS
When it comes to ultrasound images, the non-blind methods to deconvolve an image were not as
effective as the blind deconvolution methods were. The spatial filter is a technique where the new
value of an image pixel is the function of its original value, and sometimes the neighbors’ values.
Low-pass filters are used to smooth or reduce noise in an image. The mean filter is a popular low-
pass linear filter that is designed to average the pixel values inside a neighborhood. However, as
seen in Figure 1, the mean filter is not effective in reducing noise in ultrasound images. The high-
pass filters include using Laplacian masks, and we saw that it was effective in sharpening the
image, more than the aforementioned low-pass filter.
The blind deconvolution methods produced large estimate errors, but they are clearly better at
deconvolving images than the non-blind methods. The iterative multichannel blind deconvolution
algorithm incorporated denoising and regularization derived from the spatially adaptive LPA-ICI
method. This produced similar results than the second blind deconvolution method, but the
maximum likelihood technique where it incorporates estimating the PSF and then utilizing the
Wiener filter was superior. The maximum likelihood blind deconvolution starts from two
different initiation PSF estimates. Generally, the blind deconvolution results are sensitive to the
initial PSF estimate [10]. The third iteration of the restoration uses Wiener filtering with the
recovered PSF.

ACKNOLWEDGMENTS
I would like to thank my committee members for their support, and especially for the committee
chair Dr. C. Wang for faithfully meeting with me every week, even during his sabbatical period.

!"#$%@G% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
REFERENCES
1. Campisi, P., & Egiazarian, K. (2007). Blind image deconvolution: Theory and applications.
(pp. 169-228). NW: CRC Press.
2. Gonzalez, R. & Woods, R. (2001). Digital image processing. NJ: Prentice Hall.
3. Hiremath, P., & Akkasaligar, P., & Badiger, S. (2008). Speckle Noise Reduction in Medical
Ultrasound Images. Retrieved from http://www.intechopen.com/download/get/type
/pdfs/id/45101
4. Kasturiwala, S., & Ladhake, S. (2010). Superresolution: A novel application to image
restoration. International Journal on Computer Science and Engineering, 2(5), 1659-1664.
Retrieved from http://www.enggjournals.com/ijcse/doc/IJCSE10-02-05-77.pdf
5. Katvonik, V., Paliy, D., Egiazarian, K., & Astola, J. (2006). Frequency domain blind
deconvolution in multiframe imaging using anisotropic spatially-adaptive denoising.
Retrieved from http://www.cs.tut.fi/~lasip/papers/EUSIPCO2006-
FrequencyDomainBlindDeconvolutionInMultiframeImaging.pdf
6. Marques, O. (2001). Practical image and video processing using MATLAB. NJ: John Wiley&
Sons, Inc.
7. Michailovich, O. (2007). Blind deconvolution of medical ultrasound images: A parametric

inverse filtering approach. IEEE Transactions on Image Processing, 16(12).

Retrieved from http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4376242

&isnumber=4376239
8. Nakajima, N. (1993), Blind deconvolution using the maximum likelihood estimation and the

iterative algorithm. Optics Communication, 100(1), 59-66. Retrieved from

http://dx.doi.org.libproxy.usc.edu/10.1016/j.bbr.2011.03.031
!"#$%@&% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%

9. Sharma, S., Shipra, S., & Mehra, R. (2013). Image restoration using modified Lucy-
Richardson algorithm in the presence of Gaussian and motion blur. Advance in
Electronic and Electric Engineering, 3(8), 1063-1070. Retrieved from
http://www.ripublication.com/aeee.htm
10. Solomon, C., & Breckon, T. (2011). Fundamentals of digital image processing: A practical
approach with examples in MATLAB. NJ: Wiley-Blackwell.

!"#$%@@% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
BIBLIOGRAPHY
1. Campisi, P., & Egiazarian, K. (2007). Blind image deconvolution: Theory and applications.
(pp. 169-228). NW: CRC Press.
2. Gonzalez, R. & Woods, R. (2001). Digital image processing. NJ: Prentice Hall.
3. Hiremath, P., & Akkasaligar, P., & Badiger, S. (2008). Speckle Noise Reduction in Medical
Ultrasound Images. Retrieved from http://www.intechopen.com/download/get/type
/pdfs/id/45101
4. Joshi, N. (2008). PSF estimation using sharp edge prediction. In Computer vision and
pattern recognition. Retrieved from http://vision.ucsd.edu/kriegman-grp/research/
psf_estimation/psf_estimation.pdf
5. Kasturiwala, S., & Ladhake, S. (2010). Superresolution: A novel application to image
restoration. International Journal on Computer Science and Engineering, 2(5), 1659-1664.
Retrieved from http://www.enggjournals.com/ijcse/doc/IJCSE10-02-05-77.pdf
6. Katvonik, V., Paliy, D., Egiazarian, K., & Astola, J. (2006). Frequency domain blind
deconvolution in multiframe imaging using anisotropic spatially-adaptive denoising.
Retrieved from http://www.cs.tut.fi/~lasip/papers/EUSIPCO2006-
FrequencyDomainBlindDeconvolutionInMultiframeImaging.pdf
7. Katkovnik, V. (2006). Local approximations in signal and image processing (LASIP).

Retrieved from http://www.cs.tut.fi/~lasip/

8. Kundur, D., & Hatzinakos, D. (1996). Blind image deconvolution. IEEE Signal Process
Magazine, 13(3), 42-64.
9. Marques, O. (2001). Practical image and video processing using MATLAB. NJ: John Wiley&
Sons, Inc.
!"#$%@A% '%!#()*+*,$-.%/,012%#3%4#5678$51%*51%98$51%:.;#5-#80,$#5%#3%<8,+*=#051%>(*?.=% %%
10. Michailovich, O. (2007). Blind deconvolution of medical ultrasound images: A parametric

inverse filtering approach. IEEE Transactions on Image Processing, 16(12).

Retrieved from http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4376242

&isnumber=4376239
11. Nakajima, N. (1993), Blind deconvolution using the maximum likelihood estimation and the

iterative algorithm. Optics Communication, 100(1), 59-66. Retrieved from

http://dx.doi.org.libproxy.usc.edu/10.1016/j.bbr.2011.03.031

12. Sharma, S., Shipra, S., & Mehra, R. (2013). Image restoration using modified Lucy-
Richardson algorithm in the presence of Gaussian and motion blur. Advance in
Electronic and Electric Engineering, 3(8), 1063-1070. Retrieved from
http://www.ripublication.com/aeee.htm
13. Solomon, C., & Breckon, T. (2011). Fundamentals of digital image processing: A practical
approach with examples in MATLAB. NJ: Wiley-Blackwell.

Abstract (if available)

Abstract The issue of restoration of ultrasound images through blind deconvolution has been one of the chief problems in medical ultrasound imaging. This paper focuses on non‐blind and blind deconvolution methods, and compares those techniques against each other. To be specific, the non‐blind deconvolution method is based on attenuating the low‐frequency or high‐frequency components and preserving the other. In low‐pass and high‐pass filters, the ideal filters exhibit a sharp transition from passband to stopband, which contributes to the undesirable ringing effect. In contrast, the Gaussian filter exhibits a smooth transition. Butterworth filters allow the freedom to control the shape of the transition, thereby coming close to an ideal filter. The blind deconvolution method discussed in this paper involves a iterative multichannel algorithm that is derived from the spatially adaptive local polynomial approximation (LPA) and intersection of confidence levels (ICI) method. The second blind deconvolution method that is discussed is the maximum likelihood technique, which involves estimating the PSF and then utilizing the Wiener filter.

Linked assets

University of Southern California Dissertations and Theses

Asset Metadata

Creator Choi, Na Ri (author)

Core Title A comparative study of non-blind and blind deconvolution of ultrasound images

School College of Letters, Arts and Sciences

Degree Master of Arts

Degree Program Applied Mathematics

Publication Date 04/28/2014

Defense Date 03/24/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag blind image deconvolution,digital image processing,image restoration,Lucy Richardson filter,maximum likelihood,OAI-PMH Harvest,ultrasound,Wiener filter

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Wang, Chunming (committee chair), Fulman, Jason (committee member), Ziane, Mohammed (committee member)

Creator Email nachoi@usc.edu,narichoi818@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-404236

Unique identifier UC11295754

Identifier etd-ChoiNaRi-2432.pdf (filename),usctheses-c3-404236 (legacy record id)

Legacy Identifier etd-ChoiNaRi-2432.pdf

Dmrecord 404236

Document Type Thesis

Format application/pdf (imt)

Rights Choi, Na Ri

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA