(x,y)Contents:
With Image Restoration one tries to repair errors or distortions in an image, caused during the image creation process. Here we use knowledge about the manner in which the image recording took place.
In general our starting point is a degradation and noise model:
g(x,y) = H ( f(x,y) ) +
Here f(x,y) is the ideal image, the function H( ) yields the degradation caused by the projecting apparatus and (x,y) is the noise in the gray or color level. With a measured g(x,y), the real image, and some knowledge of H( ) and (x,y), we must determine the best approximation of f(x,y), the ideal image.
The degradation and the noise, and thus also the restoration, are largely determined by the recording apparatus and the conditions while acquiring the image. Restoration methods are generally mathematically complex and take a lot of time to calculate. For practical applications it is worth while to make the apparatus as degradation-free as possible, taking into account the technical and financial limitations. Good and thus usually expensive apparatuses can also limit the acquisition of noise. The better the lens for example, the more light allowed to fall on the image plane and thereby reducing noise, next to the given fact that such a lens yields less geometric distortion. Also, more expensive CCD chips, possibly with cooling, yield less noise.
With the use of CT scanners you must always consider and weigh out the consequences of the amount of radiation that is used against the amount that is harmful to the patient. More radiation yields more exposure damage to the patient but also a less noisy image by which the diagnosis and treatment can be improved. Modern x-ray machines use much less radiation than a few years ago. With MRI one shortens the recording time as much as possible because it is difficult and uneasy for the patients to spend a long time in the magnet-tube.
Of the different mathematical models of noise, see figures 5.2 and 5.4, the Gaussian model is most frequently used because it is easy to calculate with and because many sources of noise yield a Gaussian dispersal.
When the degradation process is linear: H( k1f1 + k2f2 ) = k1 H( f1) + k2 H( f2 )
we can write (we temporarily leave the noise out of consideration):
g(x,y) = H(f(x,y)) = H(
,
)
(x-
=
h(x,,y,) is the
"impulse response" or "point spread function", the indication in the image of an ideal light point. The integral is called the "superposition", the "Fredholm integral" or the "first child". When H is
a spatial invariant:
Hf( x-
we can write h( x-,y-), and the following is true:
g(x,y) =h(x,y)
a convolution integral, and taking into account the noise: G(u,v) = H(u,v)F(u,v) + N(u,v)
Now we formulate a discrete version of the spatial invariant degradation model in 1-D:
g[x] = m=0M-1 ( f[m] h[x-m] ) +
Here the rows of g and h are supplementarily filled with zeros till a length of M = A+B-1 with A the length of g and B the length of h.
In matrix notation: g = H f + n with g, f and n column vectors and H
| h[0] h[M-1] ... h[1] |
| h[1] h[0] ... h[2] |
H = | .. .. ... .. | a M*M "circulant" matrix
| h[M-1] h[M-2} ... h[0] |
Such a "circulant" matrix can be diagonalized:
H = WDW-1 where D is a diagonal matrix: D(i,j)=j
with j the eigenvalues of H that, with exception to a factor, are equal to the
Fourier transformation of h[x]. It also seems that W-1 g, W-1f and W-1n are equal to the Fourier transformation
of g[x], f[x] and [x]. We can now write:
g = WDW-1f + n thus W-1g = DW-1f + W-1n
This coincides with the following, also found above: G(u,v) = H(u,v)F(u,v) + N(u,v)
In 2-D the column vectors f, g and n are formed by M2 terms by placing the rows f[x,y], g[x,y] and [x,y] under each other. H then becomes a M2*M2 matrix, a so-called "block-circulant" matrix, which in turn can be diagonalized again: H
= WDW-1, with the same properties mentioned above.
Now, the restoration problem gets the best approximation k for finding f out of g = Hf + n, where g is the measured image, and H and n express the degradation and the noise in the image. In general there is only a limited amount of information known about this.
When H is known but nothing is known about n, then we can only try to find a solution k such that the criteria for the function J is minimal:
J( k) = || g-Hk||2 = ||n||2 ( ||v||2 = vT v )
Differentiation and equalizing to zero finally yields: (use for example (A+B)T = AT + BT and (AB)T = BT AT )
J( k) /
k = ( HT H ) -1 HT g
For a square matrix H and assuming that the inverse of the matrix exists, the equation is reduced to:
k = H-1 ( HT )-1 HT g = H-1 g
( ( HT H ) -1 HT is called the generalized inverse of H, and this always exists)
If we transpose this to a Fourier space by diagonalizing H, then:
k = ( WDW-1 ) -1 g = WD-1W-1 g
W-1 k = D-1W-1 g = W-1 g / D
thus: K(u,v) = G(u,v) / H(u,v), this is called the inverse filter or "unconstraint restoration".
Imagine that we would place a constraint on f, which is defined by a linear operator Q such that ||Qf||2 is minimal. We can require that ||Qk||2 is minimal. We also want || g-Hf||2 = ||n||2 to remain valid for our approximation of k such that || g-Hk||2 - ||n||2 is minimal. We can combine these constraints by using a Lagrange multiplier:
J( k) = ||Qk||2 +
Differentiating and setting equal to zero finally results in:
k = (HT H +QTQ)-1 HT g with
When the Lagrange multiplier is large then the constraint on Q has no influence, in the last
formula you then see that goes to 0 and that we get the old inverse filter solution back. The Lagrange
multiplier indicated the relative importance of the two minimalizations and must be chosen properly.
Sometimes is established by an iterative method such that || g-Hk||2 =
||n||2 is as valid as possible, of course this requires some knowledge about n. Later on we will focus on different possibilities for choosing Q, which can lead to
so-called "constraint restoration" methods.
The restoration by: K(u,v) = G(u,v) / H(u,v) or K(u,v) = F(u,v) + N(u,v) / H(u,v)
is called the "inverse filter". Problems arise when H(u,v) values are small or even zero, the noise term then becomes large and the image becomes dominated by noise. By restricting the frequency range to a region where H(u,v) is large enough, one can even completely avoid this problem [fig 5.27]. This is called pseudo-inverse filtering.
 |
 A 1-D Gaussian kernel for distortions in the horizontal direction. The intensity of each pixel is spread out over the neighboring pixels according to this kernel, 9 to the left and 9 to the right of the central pixel. |
An example from Digital Image Processing (DIP) with Khoros. To the left you see an image that is distorted in the horizontal direction by smudges or blurring according to the 1-D Gaussian kernel given on the right. In reality such a distortion can be created by a movement of the camera during image acquisition. |
 |
 |
In the power-spectrum of the distorted image shown on the far left, the distortion can be seen as vertical bands. The inverse filter removes the distortion adequately, which can be seen on the right. Take into account that the original image had very little noise. |
 |
 |
In the image on the left are the results of a uniform noise with values between 0 and 1 added to the image. (This was accomplished by converting the used floating-point pixel values to unsigned byte values.) In the power-spectrum on the right the vertical distortions are no longer visible. |
 |
 |
To the left you see that the application of the inverse filter yields nothing useful. The pseudo-inverse filter does yield results, the movement has been filtered out, but the amount of noise has clearly increased. The pseudo-inverse filter was done by multiplying G(u,v) with R(u,v) = 1/H(u,v) if H(u,v) > T; 0 else wise; the threshold value T is not given. |
 (original size) |
 (original size) |
 (original size) |
This example has been borrowed from Image Restoration . The image in the center has been blurred with a 2-D
Gaussian matrix with  = 3.0. From this image, by looking perpendicular to the edges, one tries to
determine the size of the smudging or blurring. Here the value of = 2.73 was found. The inverse filter
accompanied by a Gaussian function with that yields the image on the right. |
There are practical applications where H(u,v) can be obtained analytically, but has the value zero in an interesting frequency region, then there is no point in using an inverse or pseudo-inverse filter. An example of this is the linear movement in the x direction of either the camera or object, such as trying to record a race car speeding by:
g(x) = 0) d
x
And in the Fourier space:
G(u) =ux) du =
= 0
Between the {} we have the FT of a shifted function, now use f(x-b) <=> F(u) exp( -j2ub)
:
G(u) = 0
= F(u) T{ sin(
Because of the sinc function, H(u) will have the value 0 for u = n/a for any integer n.
You can also try to determine f(x) directly from g(x). Because g is an integral of f, the resulting value is the sum of the derivatives of g. Better iterative methods do exist.
We choose a Q depending on the correlation matrices of f and n:
Rf = E { f fT} and Rn = E { n nT}
QTQ = Rf-1 Rn
The element (i,j) of the matrix Rf is given by E{f[i]f[j]}, the correlation between the ith and the jth element of f; the same holds for Rn. Rf and Rn and are thus real symmetric matrices. For most images the correlation between pixels is a function of the distance between the pixels and not of their positions, furthermore the correlation has a limited range. When approximated, Rf and Rn are block-circulant matrices, which can be diagonalized. In turn, the diagonal elements are Fourier components: Sf(u,v) and Sn(u,v). Working out the constraint restoration we obtain the following result:
K(u,v)=G(u,v) H-1(u,v) |H(u,v) |2 / {|H(u,v) |2 +
When = 1 is taken we get the so-called Wiener or Least Mean Square filter. With the parameterized
Wiener filter, the is adjusted so that || g-Hk||2 =
||n||2 is as valid as possible. For = 0 we obviously get the original inverse filter,
as you can see in the formula above. The S-terms are usually not known, and are replaced by a constant K, often chosen as 2n / f~, the dispersal of the noise divided by the average signal. In [fig. 5.29] the inverse and the Wiener filter are compared for images
with different signal-noise proportions.
 |
|
 |
 |
 |
An example from Wiener Filtering, the large image can be seen in wien_smooth.GIF. To the original image on the left, additive noise has been added as can be seen in the image in the middle. The root mean square (rms) error in this image in comparison
to the original image is 18.9. To the right is the image after applying the Wiener filter, the rms error has decreased to 10.5. The rms difference between two images f1 and f2 is defined as:
(
x (f1(x) - f2(x))2 / N), the summation is over all N pixels.
In the constrained least squares restoration the value Q is chosen such that the resulting image f is as smooth as possible, that is to say:
Approximating in 1-D from 2f(x)/x2 through f(x+1) - 2f(x) + f(x+1) yields an extra criteria ||Qf||2 with Q being a band matrix. In 2-D we get a
block-circulant matrix. The solution of the constraint restoration finally becomes:
K(u,v)=G(u,v) H*(u,v) / { |H(u,v) |2 +
P(u,v) follows from Q and can be calculated directly from the image. See fig 5.30 for a result where
is chosen manually.
With an iterative method a can be found such that the constraint || g-Hk||2 =
||n||2 is satisfied as much as possible, with the only a priori knowledge being the average and the dispersal of the noise. These can also be retrieved from the image by
approximation. In [fig. 5.31] an application of this method is exemplified.
With all restoration methods, knowledge about the distortion is mandatory. Sometimes this can be extracted from knowledge about the manner in which the image is created, for example from specifications of the lens or given that the camera or an object have displaced. Often, viewing the Fourier spectrum with a trained eye can give clues about the distortions. In [ fig. 5.16, 5.19 and 5.20] some examples of this is given, there the distortions are restricted to points or a line in the Fourier spectrum. These are easily filtered away manually using a notch filter [see fig 5.18].
 |
The degradation can be of a geometric nature, caused by errors in the lenses of the projection system. Lenses often show a typical pincushion or barrel deviation. |
 |
For this you can use a simple, quick but qualitatively poor method, for example the "nearest neighbor" method: take the value of the distorted pixel where the center of the ideal pixel comes from. |
 Original image |
 Nearest neighbor |
 Bilinear interpolation |
However, you can also take the average of the 4 closest distorted pixels, weighed by the inverse of the distance to it ("bilinear interpolation"). |
The example above is from Geometric Correction of Remotely Sensed Data, the nearest neighbor clearly shows the "stair-stepped" effect, with bilinear interpolation the image is blurred out a bit because of averaging. More complex and slower, but better methods can be used, such as: determine the value of the ideal pixel according to the average weight of the area of the distorted pixels belonging to it. Determining the area can be approximated by dividing an ideal pixel up into 5*5 or 10*10 sub pixels, and using the nearest neighbor method for each of them to determine which distorted pixel belongs to it.
 |
Usually, the projection function must be determined experimentally. This can be done by looking at the projections of a set of known points, by making a grid of easily recognizable points. |
x' = a + bx + cy + dx2 + ey2 + fxy
y' = r + sx + ty + ux2 + vy2 + wxy
If only the first 3 parameters are used for x' and y' (a,b,c,r,s,t), then one speaks of an "affine" transformation, a minimum of 3 points is necessary to calculate this.
An example from Camera Lens Calibration shows that the distortion can be large:
 This image was taken with a wide-angle lens (fish eye lens). To the side the corrected image is shown; straight lines in reality are also straight in the image. |
 |
The following correction was applied:
x' = x + x*(K1*r2 + K2*r4 + K3*r6) + P1*(r2 + 2*x2) + 2*P2*x*y
y' = y + y*(K1*r2 + K2*r4 + K3*r6) + P2*(r2 + 2*y2) + 2*P1*x*y
The form follows from the type of lens that was used, the values of the parameters Ki and Pi were determined from the image itself.
With satellite images people want to correct the distortions created by the curving of the earth's surface. The method above can also be used for that. These methods can also be used to fit two
different pictures over each other, for example a satellite or aerial photo with a Kadaster map of DSA images. One can also use these methods to turn an image over a particular angle 
, the distortion function g is then defined as:
| x'| = | cos sin | | x|
| y'| | -sin cos | | y|
Updated on February 3rd 2002 by Theo Schouten.
To the left is the restoration of a blurred and noisy image using the Wiener filter, shown in 5.4. It is obviously
better than the pseudo-inverse filter result. To the right is the original image.