91
4
IMAGE SAMPLING AND
RECONSTRUCTION
In digital image processing systems, one usually deals with arrays of numbers
obtained by spatially sampling points of a physical image. After processing, another
array of numbers is produced, and these numbers are then used to reconstruct a con-
tinuous image for viewing. Image samples nominally represent some physical mea-
surements of a continuous image field, for example, measurements of the image
intensity or photographic density. Measurement uncertainties exist in any physical
measurement apparatus. It is important to be able to model these measurement
errors in order to specify the validity of the measurements and to design processes
for compensation of the measurement errors. Also, it is often not possible to mea-
sure an image field directly. Instead, measurements are made of some function
related to the desired image field, and this function is then inverted to obtain the
desired image field. Inversion operations of this nature are discussed in the sections
on image restoration. In this chapter the image sampling and reconstruction process
is considered for both theoretically exact and practical systems.
4.1. IMAGE SAMPLING AND RECONSTRUCTION CONCEPTS
In the design and analysis of image sampling and reconstruction systems, input
images are usually regarded as deterministic fields (1–5). However, in some
situations it is advantageous to consider the input to an image processing system,
especially a noise input, as a sample of a two-dimensional random process (5–7).
Both viewpoints are developed here for the analysis of image sampling and
reconstruction methods.
Digital Image Processing: PIKS Inside, Third Edition. William K. Pratt
Copyright © 2001 John Wiley & Sons, Inc.
ISBNs: 0-471-37407-5 (Hardback); 0-471-22132-5 (Electronic)
92
IMAGE SAMPLING AND RECONSTRUCTION
4.1.1. Sampling Deterministic Fields

P
xy,()F
I
xy,()Sxy,() F
I
j ∆xk∆y,()δxj∆x– yk∆y–,()
k ∞
–=
∞
∑
j ∞
–=
∞
∑
==
F
I
xy,()
j ∆xk∆y,()
F
P
ω
x
ω
y
,()
F
P
ω
x

some bounds such that for and . Performing the
convolution of Eq. 4.1-4 yields
(4.1-6)
Upon changing the order of summation and integration and invoking the sifting
property of the delta function, the sampled image spectrum becomes
(4.1-7)
As can be seen from Figure 4.1-2, the spectrum of the sampled image consists of the
spectrum of the ideal image infinitely repeated over the frequency plane in a grid of
resolution . It should be noted that if and are chosen too
large with respect to the spatial frequency limits of , the individual spectra
will overlap.
A continuous image field may be obtained from the image samples of
by linear spatial interpolation or by linear spatial filtering of the sampled image. Let
denote the continuous domain impulse response of an interpolation filter and
represent its transfer function. Then the reconstructed image is obtained
F
I
ω
x
ω
y
,() S ω
x
ω
y
,()
F
P
ω
x

j ω
xs
– ω
y
k ω
ys
–,()
k ∞
–=
∞
∑
j ∞
–=
∞
∑
=
ω
xs
2π∆x⁄= ω
ys
2π∆y⁄=
F
I
ω
x
ω
y
,()0= ω
x
ω

=
δω
x
j ω
xs
– ω
y
k ω
ys
–,()
k ∞
–=
∞
∑
j ∞
–=
∞
∑
× dα dβ
F
P
ω
x
ω
y
,()
1
∆x ∆y
--------------- F
I

R ω
x
ω
y
,()
94
IMAGE SAMPLING AND RECONSTRUCTION
by a convolution of the samples with the reconstruction filter impulse response. The
reconstructed image then becomes
(4.1-8)
Upon substituting for from Eq. 4.1-2 and performing the convolution, one
obtains
(4.1-9)
Thus it is seen that the impulse response function acts as a two-dimensional
interpolation waveform for the image samples. The spatial frequency spectrum of
the reconstructed image obtained from Eq. 4.1-8 is equal to the product of the recon-
struction filter transform and the spectrum of the sampled image,
(4.1-10)
or, from Eq. 4.1-7,
(4.1-11)
FIGURE 4.1-2. Typical sampled image spectra.
(
a
) Original image
(
b
) Sampled image
2p
∆x
2p

–=
∞
∑
=
Rxy,()
F
R
ω
x
ω
y
,()F
P
ω
x
ω
y
,()R ω
x
ω
y
,()=
F
R
ω
x
ω
y
,()
1

made identical. The first condition is met for a bandlimited image if the sampling
period is chosen such that the rectangular region bounded by the image cutoff
frequencies lies within a rectangular region defined by one-half the sam-
pling frequency. Hence
(4.1-12a)
or, equivalently,
(4.1-12b)
In physical terms, the sampling period must be equal to or smaller than one-half the
period of the finest detail within the image. This sampling condition is equivalent to
the one-dimensional sampling theorem constraint for time-varying signals that
requires a time-varying signal to be sampled at a rate of at least twice its highest-fre-
quency component. If equality holds in Eq. 4.1-12, the image is said to be sampled
at its Nyquist rate; if and are smaller than required by the Nyquist criterion,
the image is called oversampled; and if the opposite case holds, the image is under-
sampled.
If the original image is sampled at a spatial rate sufficient to prevent spectral
overlap in the sampled image, exact reconstruction of the ideal image can be
achieved by spatial filtering the samples with an appropriate filter. For example, as
shown in Figure 4.1-3, a filter with a transfer function of the form
for and (4.1-13a)
otherwise (4.1-13b)
where K is a scaling constant, satisfies the condition of exact reconstruction if
and . The point-spread function or impulse response of this
reconstruction filter is
(4.1-14)
R ω
x
ω
y
,()

y
,()
K
0





=
ω
x
ω
xL
≤ ω
y
ω
yL
≤
ω
xL
ω
xc
>ω
yL
ω
yc
>
Rxy,()
Kω

R ω
x
ω
y
,()
K
0





=
ω
x
2
ω
y
2
+ ω
0
≤
ω
0
2
ω
xc
2
ω
yc

K
J
1
ω
0
x
2
y
2
+



x
2
y
2
+
----------------------------------------=
J
1
·
{}
F
I
xy,()
η
F
I
R

y
2
–=
F
P
xy,()F
I
xy,()Sxy,()F
I
xy,() δxj∆x– yk∆y–,()
k ∞
–=
∞
∑
j ∞
–=
∞
∑
==
R
F
P
τ
x
τ
y
,()EF
P
x
1

2
,()=
98
IMAGE SAMPLING AND RECONSTRUCTION
(4.1-20)
Hence the sampled random field is also stationary with an autocorrelation function
(4.1-21)
Taking the two-dimensional Fourier transform of Eq. 4.1-21 yields the power spec-
trum of the sampled random field. By the Fourier transform convolution theorem
(4.1-22)
where and represent the power spectral densities of the
ideal image and sampled ideal image, respectively, and is the Fourier
transform of the Dirac sampling array. Then, by the derivation leading to Eq. 4.1-7,
it is found that the spectrum of the sampled field can be written as
(4.1-23)
Thus the sampled image power spectrum is composed of the power spectrum of the
continuous ideal image field replicated over the spatial frequency domain at integer
multiples of the sampling spatial frequency . If the power spectrum
of the continuous ideal image field is bandlimited such that for
and , where and are cutoff frequencies, the individual
spectra of Eq. 4.1-23 will not overlap if the spatial sampling periods are chosen such
that and . A continuous random field may be recon-
structed from samples of the random ideal image field by the interpolation formula
(4.1-24)
where is the deterministic interpolation function. The reconstructed field and
the ideal image field can be made equivalent in the mean-square sense (5, p. 284),
that is,
(4.1-25)
if the Nyquist sampling criteria are met and if suitable interpolation functions, such
as the sinc function or Bessel function of Eqs. 4.1-14 and 4.1-16, are utilized.

I
τ
x
τ
y
,()S τ
x
τ
y
,()=
W
F
P
ω
x
ω
y
,()
1
4π
2
---------W
F
I
ω
x
ω
y
,()
᭺

ω
x
ω
y
,()
1
∆x ∆y
--------------- W
F
I
ω
x
j ω
xs
– ω
y
k ω
ys
–,()
k ∞
–=
∞
∑
j ∞
–=
∞
∑
=
2π∆x 2π∆y⁄,⁄()
W

j ∆xk∆y,()Rx j∆x– yk∆y–,()
k ∞
–=
∞
∑
j ∞
–=
∞
∑
=
Rxy,()
EF
I
xy,()F
R
xy,()–
2
{}0=
IMAGE SAMPLING SYSTEMS
99
The preceding results are directly applicable to the practical problem of sampling
a deterministic image field plus additive noise, which is modeled as a random field.
Figure 4.1-4 shows the spectrum of a sampled noisy image. This sketch indicates a
significant potential problem. The spectrum of the noise may be wider than the ideal
image spectrum, and if the noise process is undersampled, its tails will overlap into
the passband of the image reconstruction filter, leading to additional noise artifacts.
A solution to this problem is to prefilter the noisy image before sampling to reduce
the noise bandwidth.
4.2. IMAGE SAMPLING SYSTEMS
In a physical image sampling system, the sampling array will be of finite extent, the

xy,()F
I
xy,()Sxy,()=
Sxy,() Px j∆x– yk∆y–,()
kK
–=
K
∑
jJ
–=
J
∑
=
Pxy,()
∆x ∆y,
Pxy,()xdyd
∞
–
∞
∫
1=
∞
–
∞
∫
D
T
xy,()
Pxy,(
)

–=
J
∑
=
F
P
xy,() F
I
j ∆xk∆ y,()Px j∆x– yk∆y–,()
kK
–=
K
∑
jJ
–=
J
∑
=
F
P
ω
x
ω
y
,()
1
4π
2
---------F
I

ω
x
J
1
2
---
+()∆x



sin
ω
x
∆x 2⁄{}sin
---------------------------------------------
ω
y
K
1
2
---
+()∆y



sin
ω
y
∆ y 2⁄{}sin
-----------------------------------------------=

sequentially overlapped samples in time.
By a simple change of variables, Eq. 4.2-9 may be rewritten as
(4.2-10)
Because only a single sampling pulse is assumed to occur during the integration
period, the limits of Eq. 4.2-10 can be extended infinitely . In this formulation, Eq.
4.2-10 is recognized to be equivalent to a convolution of the ideal continuous image
with an impulse response function with reversed coordinates,
followed by sampling over a finite area with Dirac delta functions. Thus, neglecting
the effects of the finite size of the sampling array, the model for finite extent pulse
sampling becomes
(4.2-11)
In most sampling systems, the sampling pulse is symmetric, so that .
Equation 4.2-11 provides a simple relation that is useful in assessing the effect
of finite extent pulse sampling. If the ideal image is bandlimited and A
x
and A
y
sat-
isfy the Nyquist criterion, the finite extent of the sample pulse represents an equiv-
alent linear spatial degradation (an image blur) that occurs before ideal sampling.
Part 4 considers methods of compensating for this degradation. A finite-extent
sampling pulse is not always a detriment, however. Consider the situation in which
F
S
xy,()
F
S
j ∆xk∆y,() F
I
xy,()Px j∆x– yk∆y–,()xdyd

x
∫
=
F
I
xy,() Pxy–,–()
F
S
j ∆xk∆y,()F
I
xy,()
᭺
*
Px– y–,()[]δxj∆x– yk∆y–,()=
Pxy–,–()Pxy,()=