Chapter 2

Fundamentals of Neural
Network Image Restoration

2.1 Image Distortions

Images are often recorded under a wide variety of circumstances. As imaging
technology is rapidly advancing, our interest in recording unusual or irrepro-
ducible phenomena is increasing as well. We often push imaging technology to
its very limits. For this reason we will always have to handle images suffering
from some form of degradation.
Since our imaging technology is not perfect, every recorded image is a de-
graded version of the scene in some sense. Every imaging system has a limit to
its available resolution and the speed at which images can be recorded. Often
the problems of finite resolution and speed are not crucial to the applications of
the images produced, but there are always cases where this is not so. There ex-
ists a large number of possible degradations that an image can suffer. Common
degradations are blurring, motion and noise. Blurring can be caused when an
object in the image is outside the cameras depth of field some time during the
exposure. For example, a foreground tree might be blurred when we have set up
a camera with a telephoto lens to take a photograph of a distant mountain. A
blurred object loses some small scale detail and the blurring process can be mod-
eled as if high frequency components have been attenuated in some manner in
the image [4, 21]. If an imaging system internally attenuates the high frequency
components in the image, the result will again appear blurry, despite the fact
that all objects in the image were in the camera’s field of view. Another com-
monly encountered image degradation is motion blur. Motion blur can be caused
when a object moves relative to the camera during an exposure, such as a car
driving along a highway in an image. In the resultant image, the object appears
to be smeared in one direction. Motion blur can also result when the camera
moves during the exposure. Noise is generally a distortion due to the imaging
system rather than the scene recorded. Noise results in random variations to

©2002 CRC Press LLC

pixel values in the image. This could be caused by the imaging system itself, or
the recording or transmission medium. Sometimes the definitions are not clear as
in the case where an image is distorted by atmospheric turbulence, such as heat
haze. In this case, the image appears blurry because the atmospheric distortion
has caused sections of the object to be imaged to move about randomly. This
distortion could be described as random motion blur, but can often be modeled
as a standard blurring process. Some types of image distortions, such as certain
types of atmospheric degradations [72, 73, 74, 75, 76], can be best described as
distortions in the phase of the signal. Whatever the degrading process, image
distortions may be placed into two categories [4, 21].

• Some distortions may be described as spatially invariant or space invariant.

In a space invariant distortion, the parameters of the distortion function
are kept unchanged for all regions of the image and all pixels suffer the
same form of distortion. This is generally caused by problems with the
imaging system such as distortions in the optical system, global lack of
focus or camera motion.

• General distortions are what is called spatially variant or space variant. In

a space variant distortion, the degradation suffered by a pixel in the image
depends upon its location in the image. This can be caused by internal
factors, such as distortions in the optical system, or by external factors,
such as object motion.

In addition, image degradations can be described as linear or non-linear [21].

In this book, we consider only those distortions which may be described by a
linear model.
All linear image degradations can be described by their impulse response. A
two-dimensional impulse response is often called a Point Spread Function (PSF).
It is a two-dimensional function that smears a pixel at its center with some of the
pixel’s neighbors. The size and shape of the neighborhood used by the PSF is
called the PSF’s Region of Support. Unless explicitly stated, we will from now on
consider PSFs with square shaped neighborhoods. The larger the neighborhood,
the more smearing occurs and the worse the degradation to the image. Here is
an example of a 3 by 3 discrete PSF.
 
0.5 0.5 0.5
1
0.5 1.0 0.5
5
0.5 0.5 0.5

where the factor 51 ensures energy conservation.

The final value of the pixel acted upon by this PSF is the sum of the values
of each pixel under the PSF mask, each multiplied by the matching entry in the
PSF mask.
Consider a PSF of size P by P acting on an image of size N by M. In the
case of a two-dimensional image, the PSF may be written as h(x, y; α, β). The
four sets of indices indicate that the PSF may be spatially variant hence the PSF

©2002 CRC Press LLC

will be a different function for pixels in different locations of an image. When
noise is also present in the degraded image, as is often the case in real-world
applications, the image degradation model in the discrete case becomes [4]:


N 
M
g(x, y) = f (α, β)h(x, y; α, β) + n(x, y) (2.1)
α β

where f (x, y) and g(x, y) are the original and degraded images, respectively, and
n(x, y) is the additive noise component of the degraded image. If h(x, y; α, β) is
a linear function then (2.1) may be restated by lexicographically ordering g(x, y),
f (x, y) and n(x, y) into column vectors of size NM. To lexicographically order
an image, we simply scan each pixel in the image row by row and stack them
one after another to form a single column vector. Alternately, we may scan the
image column by column to form the vector. For example, assume the image
f (x, y) looks like:  
11 12 13 14
21 22 23 24
f (x, y) = 
31 32 33 34


41 42 43 44
After lexicographic ordering the following column vector results:

f = [11 12 13 14 21 22 23 24 31 32 33 34 41 42 43 44]T

If we are consistent and order g(x, y), f (x, y) and n(x, y) and in the same
way, we may restate (2.1) as a matrix operation [4, 21]:

g = Hf + n (2.2)
where g and f are the lexicographically organized degraded and original image
vectors, n is the additive noise component and H is a matrix operator whose
elements are an arrangement of the elements of h(x, y; α, β) such that the matrix
multiplication of f with H performs the same operation as convolving f (x, y)
with h(x, y; α, β). In general, H may take any form. However, if h(x, y; α, β) is
spatially invariant with P
min(N, M ) then h(x, y; α, β) becomes h(x−α, y−β)
in (2.1) and H takes the form of a block-Toeplitz matrix.
A Toeplitz matrix [2] is a matrix where every element lying on the same
diagonal line has the same value. Here is an example of a Toeplitz matrix:
 
1 2 3 4 5
2 1 2 3 4
 
3 2 1 2 3
 
4 3 2 1 2
5 4 3 2 1

A block-Toeplitz matrix is a matrix that can be divided into a number of

equal sized blocks. Each block is a Toeplitz matrix, and blocks lying on the
same block diagonal are identical. Here is an example of a 6 by 6 block-Toeplitz

©2002 CRC Press LLC

matrix:
 
1 2 3 4 5 6
2 1 4 3 6 5  
  H11 H22 H33
3 4 1 2 3 4
  = H22 H11 H22 
4 3 2 1 4 3
  H33 H22 H11
5 6 3 4 1 2
6 5 4 3 2 1
where:
  
1 2 3 4 5 6
H11 = , H22 = , H33 =
2 1 4 3 6 5
Notice that a Toeplitz matrix is also a block-Toeplitz matrix with a block site
of 1 by 1, but a block Toeplitz matrix is usually not Toeplitz. The block-Toeplitz
structure of H comes about due to the block structure of f , g and n created by
the lexicographic ordering. If h(x, y; α, β) has a simple form of space variance
then H may have a simple form, resembling a block-Toeplitz matrix.

2.2 Image Restoration

When an image is recorded suffering some type of degradation, such as mentioned
above, it may not always be possible to take another, undistorted, image of the
interesting phenomena or object. The situation may not recur, like the image
of a planet taken by a space probe, or the image of a crime in progress. On
the other hand, the imaging system used may introduce inherent distortions to
the image which cannot be avoided, for example, a Magnetic Resonance Imaging
system. To restore an image degraded by a linear distortion, a restoration cost
function can be developed. The cost function is created using knowledge about
the degraded image and an estimate of the degradation, and possibly noise,
suffered by the original image to produce the degraded image. The free variable
in the cost function is an image, that we will denote by f̂ , and the cost function is
designed such that the f̂ which minimizes the cost function is an estimate of the
original image. A common class of cost functions is based on the mean square
error (MSE) between the original image and the estimate image. Cost functions
based on the MSE often have a quadratic nature.

2.2.1 Degradation Measure

In this work, the degradation measure we consider minimizing starts with the
constrained least square error measure [4]:
1 1
E= g − Hf̂ 2 + λDf̂ 2 (2.3)
2 2
where f̂ is the restored image estimate, λ is a constant, and D is a smoothness
constraint operator. Since H is often a low-pass distortion, D will be chosen

©2002 CRC Press LLC

to be a high-pass filter. The second term in (2.3) is the regularization term.
The more noise that exists in an image, the greater the second term in (2.3)
should be, hence minimizing the second term will involve reducing the noise in
the image at the expense of restoration sharpness.
Choosing λ becomes an important consideration when restoring an image.
Too great a value of λ will oversmooth the restored image, whereas too small a
value of λ will not properly suppress noise. At their essence, neural networks
minimize cost functions such as that above. It is not unexpected that there exist
neural network models to restore degraded imagery.

2.2.2 Neural Network Restoration

Neural network restoration approaches are designed to minimize a quadratic
programming problem [46, 105, 106, 107, 108]. The generalized Hopfield Network
can be applied to this case [35]. The general form of a quadratic programming
problem can be stated as:
Minimize the energy function associated with a neural network given by:
1 T
E=− f̂ Wf̂ − bT f̂ + c (2.4)
2
Comparing this with (2.3), W, b and c are functions of H, D, λ and n, and
other problem related constraints. In terms of a neural network energy function,
the (i, j)th element of W corresponds to the interconnection strength between
neurons (pixels) i and j in the network. Similarly, vector b corresponds to the
bias input to each neuron.
Equating the formula for the energy of a neural network with equation (2.3),
the bias inputs and interconnection strengths can be found such that as the
neural network minimizes its energy function, the image will be restored.

©2002 CRC Press LLC

 Expanding (2.3) we get:

1  1 
L L L L
E= (gp − hpi fˆi )2 + λ ( dpi fˆi )2
2 p=1 i=1
2 p=1 i=1

1   1  
L L L L L L
= (gp − hpi fˆi )(gp − hpj fˆj ) + λ ( dpi fˆi )( dpj fˆj )
2 p=1 i=1 j=1
2 p=1 i=1 j=1

1   
L L L L
= ((gp )2 − 2gp hpi fˆi + hpi fˆi hpj fˆj )
2 p=1 i=1 i=1 j=1

1  
L L L
+ λ ( dpi fˆi )( dpj fˆj )
2 p=1 i=1 j=1

1 L 
1  
L L L L L
= (gp )2 − gp hpi fˆi + hpi fˆi hpj fˆj
2 p=1 p=1 i=1
2 p=1 i=1 j=1

1  
L L L
+ λ dpi fˆi dpj fˆj
2 p=1 i=1 j=1

1  1 
L L L L L L
= hpj fˆj hpi fˆi + λ dpj fˆj dpi fˆi
2 p=1 i=1 j=1 2 p=1 i=1 j=1


L 
L
1
L
− gp hpi fˆi + (gp )2
p=1 i=1
2 p=1

Hence

1    
L 
1
L L L L L L
E= hpj hpi + λ dpj dpi fˆi fˆj − gp hpi fˆi + (gp )2
2 i=1 j=1 p=1 p=1 i=1 p=1
2 p=1
(2.5)
Expanding (2.4) we get:

1  
L L L
E=− wij fˆi fˆj − bi fˆi + c (2.6)
2 i=1 j=1 i=1

By equating the terms in equations (2.5) and (2.6) we find that the neural
network model can be matched to the constrained least square error cost function
by ignoring the constant, c, and setting:


L 
L
wij = − hpi hpj − λ dpi dpj (2.7)
p=1 p=1

©2002 CRC Press LLC

and


L
bi = gp hpi (2.8)
p=1

where wij is the interconnection strength between pixels i and j, and bi is the
bias input to neuron (pixel) i. In addition, hij is the (i, j)th element of matrix
H from equation (2.3) and dij is the (i, j)th element of matrix D.
Now let’s look at some neural networks in the literature to solve this problem.

2.3 Neural Network Restoration Algorithms

in the Literature
In the network described by Zhou et al. For an image with S + 1 gray levels,
each pixel is represented by S + 1 neurons [46]. Each neuron can have a value
of 0 or 1. The value of the ith pixel is then given by:


S
fˆi = vi,k (2.9)
k=0

where vi,k is the state of the kth neuron of the ith pixel. Each neuron is visited
sequentially and has its input calculated according to:


L
ui = b i + wij fˆj (2.10)
j=1

where ui is the input to neuron i, and fˆi is the state of the jth neuron. Based
on ui , the neuron’s state is updated according to the following rule:

∆fˆi = G(ui )

where


1, u>0
G(u) = 0, u=0 (2.11)


−1, u < 0

The change in energy resulting from a change in neuron state of ∆fˆi is given
by:
1
∆E = − wii (∆fˆi )2 − ui ∆fˆi (2.12)
2
If ∆E < 0, then the neuron’s state is updated. This algorithm may be sum-
marized as:

©2002 CRC Press LLC

Algorithm 2.1:

repeat
{
For i = 1, . . . , L do
{
For k = 0, . . . , S do
{
L
ui = bi + j=1 wij fˆj
∆fˆi = G(ui ) 

1, u>0
where G(u) = 0, u=0


−1, u < 0
∆E = − 21 wii (∆fˆi )2 − ui ∆fˆi
If ∆E < 0, then vi,k = vi,k + ∆fˆi
S
fˆi = k=0 vi,k
}
}
t=t+1
}
until fˆi (t) = fˆi (t − 1)∀i = 1, . . . , L)

In the paper by Paik and Katsaggelos, Algorithm 2.1 was enhanced to re-
move the step where the energy reduction is checked following the calculation of
∆fˆi [105]. Paik and Katsaggelos presented an algorithm which made use of a
more complicated neuron. In their model, each pixel was represented by a single
neuron which takes discrete values between 0 and S, and is capable of updating
its value by ±1, or keeping the same value during a single step. A new method
for calculating ∆fˆi was also presented:

∆fˆi = Ǵi (ui )

where


−1, u < −θi
Ǵi = 0, −θi ≤ u ≤ θi (2.13)


1, u > θi
where θi = − 21 wii > 0.

This algorithm may be presented as:

©2002 CRC Press LLC

Algorithm 2.2:

repeat
{
For i = 1, . . . , L do
{
L
ui = bi + j=1 wij fˆj
∆fˆi = Ǵi (ui ) 
−1, u < −θi

where Ǵi (u) = 0, −θi ≤ u ≤ θi


1, u > θi
where θi = − 21 wii > 0
fˆi (t + 1) = K(fˆi (t) + ∆fˆi )

0, u < 0
where K(u) = u, 0 ≤ u ≤ S


S, u ≥ S
}
t=t+1
}
until fˆi (t) = fˆi (t − 1)∀i = 1, . . . , L)

Algorithm 2.2 makes no specific check that energy has decreased during each
iteration and so in [105] they proved that Algorithm 2.2 would result in a decrease
of the energy function at each iteration. Note that in Algorithm 2.2, each pixel
only changes its value by ±1 during an iteration. In Algorithm 2.1, the pixel’s
value would change by any amount between 0 and S during an iteration since
each pixel was represented by S + 1 neurons. Although Algorithm 2.2 is much
more efficient in terms of the number of neurons used, it may take many more
iterations than Algorithm 2.1 to converge to a solution (although the time taken
may still be faster than Algorithm 2.1). If we consider that the value of each pixel
represents a dimension of the L dimensional energy function to be minimized,
then we can see that Algorithms 2.1 and 2.2 have slightly different approaches
to finding a local minimum. In Algorithm 2.1, the energy function is minimized
along each dimension in turn. The image can be considered to represent a single
point in the solution space. In Algorithm 2.1, this point moves to the function
minimum along each of the L axes of the problem until it eventually reaches
a local minimum of the energy function. In Algorithm 2.2, for each pixel, the
point takes a unit step in a direction that reduces the network energy along that
dimension. If the weight matrix is negative definite (−W is positive definite),
however, regardless of how these algorithms work, the end results must be similar
(if each algorithm ends at a minimum). The reason for this is that when the
weight matrix is negative definite, there is only the global minimum. That is,
the function has only one minimum. In this case the matrix W is invertible and
taking (2.4) we see that:

©2002 CRC Press LLC

 δE
= −Wf̂ − b (2.14)
δ f̂
Hence the solution is given by:

f̂  = −W−1 b (2.15)
(assuming that W−1 exists).
The f̂  is the only minimum and the only stationary point of this cost func-
tion, so we can state that if W is negative definite and Algorithm 2.1 and Algo-
rithm 2.2 both terminate at a local minimum, the resultant image must be close
to f̂  for both algorithms. Algorithm 2.1 approaches the minimum in a zigzag
fashion, whereas Algorithm 2.2 approaches the minimum with a smooth curve.
If W is not negative definite, then local minimum may exist and Algorithms 2.1
and 2.2 may not produce the same results. If Algorithm 2.2 is altered so that
instead of changing each neuron’s value by ±1 before going to the next neuron,
the current neuron is iterated until the input to that neuron is zero, then Algo-
rithms 2.1 and 2.2 will produce identical results. Each algorithm will terminate
in the same local minimum.

2.4 An Improved Algorithm

Although Algorithm 2.2 is an improvement on Algorithm 2.1, it is not optimal.
From iteration to iteration, neurons often oscillate about their final value, and
during the initial iterations of Algorithm 2.1 a neuron may require 100 or more
state changes in order to minimize its energy contribution. A faster method to
minimize the energy contribution of each neuron being considered is suggested
by examination of the mathematics involved. For an image where each pixel is
able to take on any discrete integer intensity between 0 and S, we assign each
pixel in the image to a single neuron able to take any discrete value between 0
and S. Since the formula for the energy reduction resulting from a change in
neuron state ∆fˆi is a simple quadratic, it is possible to solve for the ∆fˆi which
produces the maximum energy reduction. Theorem 2.1 states that this approach
will result in the same energy minimum as Algorithm 2.1 and hence the same
final state of each neuron after it is updated.

Theorem 2.1: For each neuron i in the network during each iteration, there
exists a state change ∆fˆi∗ such that the energy contribution of neuron i is mini-
mized.

Proof:

Let ui be the input to neuron i which is calculated by:

L
ui = b i + j=1 wij fˆj

©2002 CRC Press LLC

 Let ∆E be the resulting energy change due to ∆fˆi .

1
∆E = − wii (∆fˆi )2 − ui ∆fˆi (2.16)
2
Differentiating ∆E with respect to ∆fˆi gives us:

δ∆E
= −wii ∆fˆi − ui
δ fˆi
The value of ∆fˆi which minimizes (2.16) is given by:

0 = −wii ∆fˆi∗ − ui

Therefore,
−ui
∆fˆi∗ = (2.17)
wii
QED.

Based on Theorem 2.1, an improved algorithm is presented below.

Algorithm 2.3.

repeat
{
For i = 1, . . . , L do
{
L
ui = bi + j=1 wij fˆj
∆fˆi = G(ui ) 

−1, u < 0
where G(u) = 0, u=0


1, u>0

1
∆Ess = − wii (∆fˆi )2 − ui ∆fˆi (2.18)
2

−ui
If ∆Ess < 0 then ∆fˆi∗ =
wii
ˆi (t) + ∆fˆ∗ )
fˆi (t + 1) = K(f i

0, u < 0
where K(u) = u, 0 ≤ u ≤ S


S, u ≥ S
}

©2002 CRC Press LLC

 t=t+1
}
until fˆi (t) = fˆi (t − 1)∀i = 1, . . . , L)

Each neuron is visited sequentially and has its input calculated. Using the
input value, the state change needed to minimize the neuron’s energy contribu-
tion to the network is calculated. Note that since ∆fˆi ∈ {−1, 0, 1} and ∆fˆi and
∆fˆi∗ must be the same sign as ui , step (2.18) is equivalent to checking that at
least a unit step can be taken which will reduce the energy of the network. If
∆Ess < 0, then

− 21 wii − ui ∆fˆi < 0

− 21 wii − |ui | < 0
−wii < 2|ui |

Substituting this result into the formula for ∆fˆi∗ we get:

−ui ui 1
∆fˆi∗ = > = ∆fˆi
wii 2|ui | 2
Since ∆fˆi∗ and ∆fˆi have the same sign and ∆fˆi = ±1 we obtain:
1
|∆fˆi∗ | > (2.19)
2
In this way, ∆fˆi∗ will always be large enough to alter the neuron’s discrete
value.
Algorithm 2.3 makes use of concepts from both Algorithm 2.1 and Algorithm
2.2. Like Algorithm 2.1 the energy function is minimized in solution space along
each dimension in turn until a local minimum is reached. In addition, the effi-
cient use of space by Algorithm 2.2 is utilized. Note that the above algorithm
is much faster than either Algorithm 2.1 or 2.2 due to the fact that this algo-
rithm minimizes the current neuron’s energy contribution in one step rather than
through numerous iterations as did Algorithms 2.1 and 2.2.

2.5 Analysis
In the paper by Paik and Katsaggelos, it was shown that Algorithm 2.2 would
converge to a fixed point after a finite number of iterations and that the fixed
point would be a local minimum of E in (2.3) in the case of a sequential algo-
rithm [105]. Here we will show that Algorithm 2.3 will also converge to a fixed
point which is a local minimum of E in (2.3).
Algorithm 2.2 makes no specific check that energy has decreased during each
iteration and so in [105] they proved that Algorithm 2.2 would result in a decrease
of the energy function at each iteration. Algorithm 2.3, however, changes the
current neuron’s state if and only if an energy reduction will occur and |∆fˆi | = 1.
For this reason Algorithm 2.3 can only reduce the energy function and never

©2002 CRC Press LLC

increase it. From this we can observe that each iteration of Algorithm 2.3 brings
the network closer to a local minimum of the function. The next question is
Does Algorithm 2.3 ever reach a local minimum and terminate? Note that the
gradient of the function is given by:

δE
− Wf̂ − b = −u (2.20)
δ f̂
where u is a vector whose ith element contains the current input to neuron i.
Note that during any iteration, u will always point in a direction that reduces
the energy function. If f̂ = f̂  then for at least one neuron a change in state must
be possible which would reduce the energy function. For this neuron, ui = 0.
The algorithm will then compute the change in state for this neuron to move
closer to the solution. If |∆fˆi∗ | > 21 the neuron’s state will be changed. In this
case we assume that no boundary conditions have been activated to stop neuron
i from changing value. Due to the discrete nature of the neuron states we see
that the step size taken by the network is never less than 1.
To restate the facts obtained so far:

• During each iteration Algorithm 2.3 will reduce the energy of the network.

• A reduction in the energy of the network implies that the network has
moved closer to a local minimum of the energy function.

• There is a lower bound to the step size taken by the network and a finite
range of neuron states. Since the network is restricted to changing state
only when an energy reduction is possible, the network cannot iterate for-
ever.

From these observations we can conclude that the network reaches a local
minimum in a finite number of iterations, and that the solution given by Al-
gorithm 2.3 will be close to the solution given by Algorithm 2.1 for the same
problem. The reason Algorithms 2.1 and 2.3 must approach the same local
minimum is the fact that they operate on the pixel in an identical manner. In
Algorithm 2.1 each of the S + 1 neurons associated with pixel i is adjusted to
reduce its contribution to the energy function. The sum of the contributions
of the S + 1 neurons associated with pixel i in Algorithm 2.1 equals the final
grayscale value of that neuron. Hence during any iteration of Algorithm 2.1
the current pixel can change to any allowable value. There are S + 1 possible
output values of pixel i and only one of these values results when the algorithm
minimizes the contribution of that pixel. Hence whether the pixel is represented
by S + 1 neurons or just a single neuron, the output grayscale value that occurs
when the energy contribution of that pixel is minimized during the current itera-
tion remains the same. Algorithms 2.1 and 2.3 both minimize the current pixel’s
energy contribution; hence they must both produce the same results. In practice
the authors have found that all three algorithms generally produce identical re-
sults, which suggests that for reasonable values of the parameter λ, only a single
global minimum is present.

©2002 CRC Press LLC

 Note that in the discussion so far, we have not made any assumptions re-
garding the nature of the weighting matrix, W, or the bias vector, b. W and
b determine where the solution is in the solution space, but as long as they
are constant during the restoration procedure the algorithm will still terminate
after a finite number of iterations. This is an important result, and implies that
even if the degradation suffered by the image is space variant or if we assign a
different value of λ to each pixel in the image, the algorithm will still converge to
a result. Even if W and b are such that the solution lies outside of the bounds
on the values of the neurons, we would still expect that there exists a point or
points which minimize E within the bounds. In practice we would not expect
the solution to lie entirely out of the range of neuron values. If we assume that
Algorithm 2.3 has terminated at a position where no boundary conditions have
been activated, then the condition:
   
 ˆ∗   −ui  1
∆fi  =  < , ∀i ∈ {0, 1, . . . , L}
wii  2
must have been met. This implies that:

1
|ui | < wii , ∀i ∈ {0, 1, . . . , L} (2.21)
2
In [105], Paik and Katsaggelos noticed this feature as well, since the same
termination conditions apply to Algorithm 2.2. The self-connection weight, wii ,
controls how close to the ideal solution the algorithm will approach before ter-
minating. Since increasing the value of λ increases the value of wii , we would
expect also that the algorithm would terminate more quickly and yet be less
accurate for larger values of λ. This is found to occur in practice. When λ is
increased, the number of iterations before termination drops rapidly.

2.6 Implementation Considerations

Despite the increase in efficiency and speed of Algorithm 2.3 when compared to
Algorithms 2.1 and 2.2, there are still a number of ways that the algorithm can
be made more efficient. The ith row of the weighting matrix describes the inter-
connection strengths between neuron i and every other neuron in the network
from the viewpoint of neuron i. The weighting matrix is N M by N M , which is
clearly a prohibitively large amount of data that requires storage. However, the
mathematical discussion in the previous sections implies a shortcut.
By examining (2.7) we observe that in the case of P
min(M, N ), it can be
seen that when calculating the input to each neuron, only pixels within a certain
rectangular neighborhood of the current pixel contribute non-zero components
to the neuron input. In addition, it can be seen that the non-zero interconnection
strengths between any given pixel and a neighboring pixel depend only on the
position of the pixels relative to each other in the case of spatially invariant
distortion. Using the above observations, the input to any neuron (pixel) in
the image can be calculated by applying a mask to the image centered on the

©2002 CRC Press LLC

pixel being examined. For a P by P distortion, each weighting mask contains
only (2P − 1)2 terms. A 5 by 5 degrading PSF acting on a 250 by 250 image
requires a weight matrix containing 3.9 × 109 elements, yet a weighting mask of
only 81 elements. In addition, by considering the finite regions of support of the
degrading and filtering impulse functions represented by H and D, the weighting
masks and bias inputs to each neuron may be calculated without storing matrices
H and D at all. They may be calculated using only the impulse responses of the
above matrices.

2.7 A Numerical Study of the Algorithms

To compare the three algorithms, let us look at two examples. In the first
example, the efficiency of Algorithms 2.1, 2.2 and 2.3 will be compared to one
another. In the second example, a practical example of the use of this method
will be given.

2.7.1 Setup
In both these examples, the images were blurred using a Gaussian PSF with the
impulse response:
  2 
1 x y2
h(x, y) = exp − + 2 (2.22)
2πσx σy 2σx2 2σy
where σx and σy are the standard deviations of the PSF in the x and y directions,
respectively.

2.7.2 Efficiency
The time taken to restore an image was compared among Algorithms 2.1, 2.2
and 2.3. A degraded image was created by blurring a 256 by 256 image with a
Gaussian blur of size 5 by 5 and standard deviation 2.0. Noise of variance 4.22
was added to the blurred image. Each algorithm was run until at least 85% of the
pixels in the image no longer changed value or many iterations had passed with
the same number of pixels changing value during each iteration. Algorithm 2.1
was stopped after the sixth iteration when no further improvement was possible,
and took 6067 seconds to run on a SUN Ultra SPARC 1. Algorithm 2.2 was
stopped after the 30th iteration with 89% of pixels having converged to their
stable states and took 126 seconds to run. Algorithm 2.3 was stopped after the
18th iteration with 90% of pixels stable and took 78 seconds to run. Algorithm
2.3 is much faster than Algorithms 2.1 and 2.2, despite the fact that algorithms
2.1 and 2.3 approach the same local minimum and hence give the same results.
The computation time of Algorithm 2.3 can be expected to increase linearly with
the number of pixels in the image, as can the computation times of Algorithms 2.1
and 2.2. The single step neuron energy minimization technique of Algorithm 2.3
provides its superior speed and this trend holds for any size image. Various types

©2002 CRC Press LLC