- Email: [email protected]
Image reconstruction by parametric cubic convolution
COMPU’IER
VISION,
GRAPHICS,
AND
IMAGE
PROCESSING
23,
258-212
Image Reconstruction by Parametric Cubic Convolution STEPHEN K. PARK NASA,
Langley
Reseurch
Center,
Hampton,
Virginia
23465
AND ROBERT A. SCHOWENGERDT Office
of Arid
La&
Studies
and Electrical Tucson,
Engineering Department, Arizona 85721
University
of Arizona,
Received May 2 1. 1982; revised July 27, 1982 A parametric implementation of cubic convolution image reconstruction is presented which is generally superior to the standard algorithm and which can be optimized to the frequency content of the image. I. INTRODUCTION
Cubic convolution is a process originally developed for the reconstruction (i.e., interpolation, resampling) of Landsat digital images [ 11. Since its development, this process has come to be generally acknowledged as providing a good compromise between computational complexity and reconstruction accuracy [2]. Although it is recognized that cubic convolution image reconstruction is actually interpolation implemented by convolving the image samples with a smooth, spatially limited, one-parameter family of functions constructed by joining cubic polynomials in a splinelike fashion [3, 41, traditionally the parametric nature of cubic convolution has not been exploited and one member of this family has received virtually all the attention; it is the interpolation function originally proposed [l] and it remains the standard for comparison [5] and software implementation [6]. In this paper cubic convolution is formulated and analyzed parametrically, i.e., the reconstruction properties of the one-parameter family of cubic convolution interpolation functions are discussed and the image degradation (i.e., sampling and reconstruction blur [7]) associated with reasonable choices of this parameter is analyzed. It is demonstrated, with an analysis in the frequency domain, that in an image-independent sense there is an optimal value for this parameter and it is not the standard value commonly referenced in the literature. This important conclusion is consistent with a recently published result [8], which is based upon a complementary spatial domain analysis. In addition, it is demonstrated that in an image-dependent sense, cubic convolution can be tailored to any class of images characterized by a common energy spectrum, i.e., a different optimal parameter value exists for each different image energy spectrum. Thus, this paper establishes theoretically that a parametric implementation of cubic convolution image reconstruction is possible which is generally superior to the standard and which can be adapted to the frequency content of the image. 258 0734-189X/83
$3.00
Copyright 0 1983 by Academic Press, Inc. All rights of reproduction in any form reserved.
PARAMETRIC 2. IMAGE
CUBIC CONVOLUTION
259
RECONSTRUCTION
In a recent article [7], the authors analyzed the degradation caused by image sampling and interpolative reconstruction. The present paper uses the results of this analysis as well as the same notation and terminology; for simplicity the analysis is one-dimensional. A typical imaging, sampling and reconstruction system is illustrated in Fig. 1. The input scene (or object) is denoted f( x), the PSF of the imaging subsystem is denoted h(x), and the output of the imaging subsystem (prior to sampling) is the image g(x). Thus
g(x)=m*w
(1)
where * denotes convolution. As in [7], the spatial coordinate x is normalized in units of sample interval and the output of the sampling subsystem is the sampled image g,(x) given by g,(x)
= dx)comb(x)
(24
where comb(x)
=
F 8(x-n) n=-CC
(‘W
is the familiar sampling function [9]. Image reconstruction is commonly modeled as [lo]
g,(x) = g,(x)* 44
(3)
where the interpolation function r(x) is the reconstruction subsystem’s response to a sampled input which is one at x =0 and zero at x = n = rtl, t-2,... . The objective of image reconstruction is to choose T(X) so that the reconstructed image g,(x) accurately reproduces the output of the imaging subsystem g(x). Image reconstruction should not be confused with (linear) image restoration [I l] whose objective is to choose r(x) so that g,(x) accurately reproduces the input to the imaging subsystem f(x). Image reconstruction is an interpolation process in which the image samples (Eq. (2 a )) are reproduced exactly and thus, in general, r(0) = 1 r(n)
f(x) v SCENE
FIG.
IMAGING s!(x) SUBSYSTEM +h(x) IMAGE
= 0
n = +l,
+2, . . . .
SAMPLING g$) RECONSTRUCTION q,(x) a SUBSYSTEM SUBSYSTEH .comb(x) SAMPLED ?r(x) RECDNSTRUCTED - IMAGE I HAGF
1. Block diagram of a general imaging, sampling, and reconstruction system.
(4)
260
PARK
AND
SCHOWENGERDT
Equation (3) can be written equivalently as
d-4 = E s(n)rb - 4 n=-CC
(5)
which illustrates that, for a fixed x, the extent to which neighboring image samples contribute to the reconstructed image is determined by the spread of the interpolation function. For digital image reconstruction, computational considerations dictate that this spread should be small; typically, T(X) is identically zero for 1x1larger than two. It is well known [9-121 that if the image g(x) is band-limited and sufficiently sampled, then it can be reconstructed exactly by using the ideal interpolation function r(x)
sin( 7rx ) = sine(x) = 7rx .
In practice, exact digital reconstruction is impossible because the ideal interpolation function is not spatially limited and the sum in Eq. (5) must be truncated by using a truncated form of Eq. (6) as the interpolation function. However, truncating Eq. (6) to a practical width, say 1x1 I 2, introduces undesirable ripples in the spectrum (i.e., Fourier transform) of the resultant interpolator [5] because of the slope discontinuity at x = + 2. One method of suppressing these ripples is to construct the interpolation function as a smooth, spatially limited approximation to sine(x) on the interval - 2 5 x I 2 and eliminate the slope discontinuity at the end points; cubic convolution is an implementation of this method. 3. PARAMETRIC
CUBIC
CONVOLUTION
(PCC)
The interpolation function commonly associated with cubic convolution image reconstruction is the piecewise cubic polynomial illustrated in Fig. 2. By inspection this function (i) is even, (ii) is identically zero for all Ix] 2 2; (iii) is smooth, i.e., continuous with a continuous slope for all x, and (iv) satisfies Eq. (4). In addition, it can be verified that this function satisfies an important fifth property: (v) for UN x f
r(x-n)=
1.
n=-W Equation (7) is satisfied if and only if the image reconstruction exactly a constant (i.e., uniform) image.
FIG.
comparison,
2. The
interpolation function commonly associated with cubic the ideal interpolation function, sinc( x). is also illustrated.
process reproduces
convolution
((Y = - 1). car
PARAMETRIC
CUBIC
261
CONVOLUTION
It is known [3] that the function in Fig. 2 is just one member (a = - 1) of a one-parameter family of piecewise cubic polynomials given by (a + 2)1x13 - (a + 3)1x12 + 1
44 = (Ylx13 - 5a1x1* + thqxl-
1x1< 1 1 I Ix] < 2 otherwise.
4a
0
(8)
The parameter LYhas physical significance; it is the slope of T(X) at x = 1. The common value [3-61 is (r = - 1, which duplicates the slope of the ideal interpolator (Eq. (6)) at x = 1 (see Fig. 2); however, it was recently proven [8] that cx= -0.5 is superior to (Y= - 1 in that it provides an interpolator with better convergence properties. An intermediate choice, cx= -0.75, is also sometimes advocated; this choice results from forcing Eq. (8) to have second derivative continuity at x = 1. It can be shown directly that Eq. (8) satisfies conditions (i) through (v) for all (Y. Various members of this family of interpolators are illustrated in Fig. 3 correspondingtoa= l,O, - 0.5, - 1, and - 2. Since T(X) is an even function, only values for x 2 0 are illustrated. Note that for (Y = 0, T(X) is nonnegative and identically zero for 1x1 2 1. Equation (8) can be written equivalently as
r(x) = To(X)+ q(x) where
cw + Nlxl - l)*
Ix] < 1 otherwise
0
(9b)
and
dx)
=
14*(1x1- 1)
I4 < 1
(1x1 - l)(]x] - 2)* 0
1 I Ix] I 2 otherwise.
I.0
c
0.5
r x
0
-0.5 FIG.
3. Five
members
of the family
of PCC
interpolators.
(9c)
262
PARK
AND
SCHOWENGERDT
Equation (9a) illustrates clearly the linear dependence of this family of interpolators upon the parameter a. The function To(x), illustrated in Fig. 4, can be thought of as a smooth approximation to the triangular function associated with linear interpolation and r,(x), also shown in Fig. 4, is the smooth correction which must be added to r,,(x) to yield a function r(x) = To(x) + w,(x) whose slope at x = 1 is (Y. Equation (9) satisfies conditions (i) through (v) and, in particular, Eq. (7) becomes f
ro(x-n)=
(104
1
and fj t-,(x - n) = 0 n=-CC
for all x. Equations (5) and (9) define parametric reconstruction process, i.e.,
(10’4
cubic convolution
(PCC) as an image
where
h(x) = g,(x)*r0(x)= “=-ccl E b+hb - n>
(11’4
and
g,(x) = dx)*r,(x)
= it
dnhb
- 4.
Because the interpolation functions To(x) and rI(x) are spatially limited to two and four sample intervals, respectively, and because they are just piecewise cubic polynomials, PCC represents a computationally efficient algorithm. In addition, this algorithm gives the user the ability to control the visual appearance (or radiometric
FIG. 4. The two
componentsof
the family
of PCC interpolators.
PARAMETRIC
263
CUBIC CONVOLUTION
q,(x):
STANDARD CUBIC CONVOLUTION (a =-
FIG. 5. A one-dimensional image reconstructed from 31 samples for several values of a: 0, -0.5. - 1. Also shown is g,(x), the high frequency image component which adds local gradient information.
distribution) of the reconstructed image by changing the PCC parameter (Y. From Eq. (1 l), all PCC images for a given sampled image are linear combinations of the components g,(x) and g,(x). Thus an interactive search for the “best” reconstructed image can be implemented quite efficiently. Figure 5 illustrates PCC image reconstruction for a representative portion of a sampled (one-dimensional) image. The 31 image samples are interpolated by PCC for three values of (Y: 0, - 0.5, and - 1. The reconstructed image corresponding to (Y= 0 is just g,(x); the other two reconstructed images are formed from g,(x) by adding ag, (x). Of the three reconstructed images, the one for (r = - 0.5 is best; for (Y = - 1, edge overshoots occur at x = 23.5 and 26.5 and artificial stepping occurs for x between 5.0 and 11.0; for LY= 0, the slope of go(x) is zero at all image samples causing the appearance of artifacts, particularly in the range 5.0 < x < 11 .O. On the basis of Fig. 5 and the results of Ref. [8], we conjecture that for a large class of images, PCC with (Y= - 0.5 is superior to both standard cubic convolution ((Y = - 1) and to the smooth approximation to linear interpolation provided by (Y = 0. Later in this paper, an analysis in the frequency domain is presented which confirms this conjecture. 4. RECONSTRUCTION
FILTER
The interpretation of Eq. (3) in the frequency (Fourier) domain is that image reconstruction can be thought of as a process in which the sampled image spectrum is multiplied by a reconstruction filter i(v) where i(v)
= /”
-m
r(x)e-2nxvidx
(12)
is the Fourier transform of the interpolation function T(X). Since x is normalized in units of sample interval, the frequency coordinate v has units of cycles per sample
264
PARK
FIG.
6. Five
members
AND
SCHOWENGERDT
of the family
of PCC
reconstruction
filters
interval. In this coordinate system the Nyquist (or folding) frequency is 0.5. It can be shown that the Fourier transform of Eq. (9) is
i(v) = to(v) + “i,(V)
(134
where
+I(~)= -&*
[sinc*( y) - sinc(2v)]
(13’4
and
i,(v) = - 2
[3 sinc2(2v) - 2 sinc(2Y) - sinc(4v)l.
(13c)
(4’
Figure 6, which is the frequency domain counterpart of Fig. 3, illustrates i(v) for five values of (Y: 1, 0, - 0.5, - 1, and - 2. For all CX,i(v) is a low-pass filter. Figure 7, the frequency domain counterpart to Fig. 4, illustrates tO(v) and i,(v). Since all the functions in Eq. (13a) are even functions of v, only values of v 2 0 are illustrated.
FIG.
7. The
two components
of the family
of PCC
reconstruction
filters
PARAMETRIC
265
CUBIC CONVOLUTION
From Fig. 6 it can be seen that the value (Y = - 0.5 has particular significance; for a! < - 0.5, i(v) is a reconstruction filter which enhances low frequencies while for (Y > - 0.5, i(v) suppresses low frequencies. This effect can be examined in more detail by expanding i(v) in powers of (av) which yields 3(v) = 1 - $2,
+ l)(~v)*
+ A(1601 + l)(~v)~
+ ... .
From Eq. (14) i(0) = 1 and i(v) is either curved up or down at v = 0 depending upon whether (~CX+ 1) is negative or positive. The value (Y= - 0.5 is optimal in that it eliminates the (TV)* term in Eq. (14) and provides the best low frequency approximation to the ideal reconstruction filter, i.e., to the Fourier transform of Eq. (6) which is Iv( < 0.5 otherwise. Note, from Fig. 6, that the choice (Y= -0.5 also provides an excellent high frequency approximation to Eq. (15), i.e., for (Y = -0.5, i(v) (Eq. (13)) is effectively zero for all Iv] > 1.0. Thus the reconstruction filter corresponding to (Y= -0.5 is also particularly effective at suppressing aliasing. 5. SAMPLING
AND RECONSTRUCTION
BLUR
It is well known that for practical systems, image sampling and reconstruction inevitably produce some degradation in the reconstructed image [lo], i.e., g,(x) (Eq. (3)) is not an exact copy of g(x) (Eq. (1)). One measure of this degradation (mean square radiometric error) (16)
called sampling and reconstruction blur (SR blur) in reference [7], is particularly amenable to analysis. For example, it can be shown that: (i) if the image is band-limited, i.e., the image energy spectrum, ] g( v)] *, is ’ zero for all ]v( greater than a cut-off frequency vC; and (ii) if the image is sufficiently sampled, i.e., v, is less than or equal to the Nyquist frequency, 0.5 cycles per sample interval, then Eq. (16) reduces to E& = /“’ e2(v)lg(v)lz - “c
dv
where e*(v)
= II - i(v
+ C li(v - n)l*. nto
(17b)
If the band-limited image is not sufficiently sampled, then SR blur depends upon the sample-scene phase, i.e., the location of the sampling grid relative to the scene. However, Eq. (17) remains valid in the sense that it represents the average SR blur over the ensemble of all possible (random, equally likely) sample-scene phases.
266
PARK AND SCHOWENGERDT
For PCC image reconstruction, Eq. (17b) defines e*(v) in terms of the reconstruction filter, Eq. (13). As Eq. (17a) reveals, the cause of SR blur is the presence of significant image energy at frequencies where e*(v) is not zero. With reference to Fig. 1, just as the image subsystem can be interpreted as a presampling filter acting upon the scene energy spectrum ]f(v)l* to produce image blur, the reconstruction subsystem can be interpreted as a postsampling filter acting upon the image energy spectrum ]g( v)]* to produce SR blur. Note that e*(v) is not just ]I - i(v)]*; sampling causes the appearance of the side-band terms Cn tO ]i( v - n )I *. The importance of the parameter (Yis that it provides a means of shaping the e*(v) vs v curve and thereby controlling SR blur. 6. OPTIMAL
a: IMAGE-INDEPENDENT
CASE
Since the reconstruction filter (Eq. (13)) is real, Eq. (17b) can be simplified by omitting the magnitude notation and including i*(v) in the infinite series to yield e*(v) = 1 - 2?(v) +
f P(v - n). n=--M
(18)
Because the associated interpolation function, Eq. (9) is spatially limited to the interval Ix] 5 2, it can be shown [7] that the infinite series in Eq. (18) is equivalent to i*(v
f
- n) = (r * r)(O) + 2 i
n=-CC
(r * r)(n)cos(2mnv)
n=l
(19)
where (r * r)(n) represents r(x) * r(x) evaluated at x = n. The coefficients (r * r)(n) can be determined algebraically and combined with Eq. (19) to yield E
i*(v
- n) = 1 - J$sin*(nu)
- $$sin*(2nv)
n=-CC + gsin*(2*v)[l
+ 6sin*(nv)].
Since Eq. (20) is quadratic in (Y, e*(v) will be also. Note that e*(v) is an even function of v. For most digital images the preponderance of energy exists at low frequencies. This implies that in an image-independent sense the optimal choice of (Yis that one which yields the best low frequency approximation to the equation e*(v) = 0. This optimality criteria is reinforced by the observation that, for the ideal reconstruction filter (Eq. (15)), e*(v) is identically zero for Iv] < 0.5. If Eqs. (18) and (20) are expanded in powers of (TV) and combined, the resulting low frequency expansion is e*(v)
= -&(2a
+ l)‘(pv)*
+ &-(2,
+ 1)(14a! + 9)(7r~)~ + ... .
Equation (21) is a particularly satisfying result because it again emphasizes the (image-independent) optimality of (Y= - 0.5; if (2a + 1) = 0, then both the (TV)* and (a~)~ terms are eliminated and e*(v) = (constant)( rv)6
(22)
PARAMETRIC
CUBIC
CONVOLUTION
267
for low frequencies. (We have verified numerically that the constant in Eq. (22) is not zero.) In Ref. [7], cubic convolution (with (Y = - 1) was compared with the standard nearest neighbor and linear image reconstruction techniques. The results of this analysis may be used to derive low frequency expansions for nearest neighbor
e’(v)= &7v)2- &rv)4 +&+v)6 + . ..
(234
and for linear
f?‘(v)= &7v)” - $&rv)” + *.. . By comparing Eqs. (21), (22), (23a), and (23b), we see that for images dominated by low frequency energy, linear interpolation is superior to nearest neighbor and PCC is superior to linear but only if a = - 0.5. Figure 8 is a plot of e(y) (not e’(v)) vs v for four interpolators: nearest neighbor, linear, standard cubic convolution (a = - I), and PCC with (Y = -0.5. The curves for nearest neighbor, linear, and PCC with (Y = -0.5 are similar in shape but different in magnitude; the curve for standard cubic convolution is markedly different with a low frequency bump caused by the nonoptimality of (Y= - 1. Since
FIG.
convolution frequency
8. The function e(v) for four common interpolators, nearest neighbor, linear, standard cubic ((1 = - I), and PCC with a = - 0.5. The nonoptimality of a = - 1 is associated with the low bump in the curve. The Nyquist frequency is 0.5 cycles per sample interval.
268
PARK AND SCHOWENGERDT
I
10‘3 0
FIG.
I, 0.2
I 0.4
I I
I,
I
0.6
0.8
I
I I.0
9. The function e(v) for PCC with a = 1, 0, - 0.5, - I, and - 2.
e(v) for PCC with (Y = - 0.5 is less than either linear or nearest neighbor over the range Iv] I 0.5, we conclude that for any band-limited and sufficiently sampled image, parametric cubic convolution with a = - 0.5 yields less SR blur than either linear or nearest neighbor interpolation. The irony of this important theoretical result is that it is not also true for the standard value (Y = - 1, e.g., for this value, and for any image
with all its energy at frequencies below v = 0.12, linear interpolation would yield less SR blur. The insights provided by Eqs. (21), (22), (23a), and (23b) are consistent with the results of Ref. [8], based upon a local spatial domain analysis. In Ref. [8] it was proven that cubic convolution is best in that it is third-order convergent if (Y= - 0.5; linear interpolation is second-order convergent and thus is the intermediate choice; nearest neighbor is the worst because it is only first-order convergent. For (Y* - 0.5, cubic convolution is just first-order convergent. Figure 9 is a plot which illustrates e(v) for PCC corresponding to the five values of LY: 1, 0, - 0.5, - 1, and - 2. As this figure demonstrates, a parametric implementation of cubic convolution is desirable because (Ycan be used to shape e*(v) and thereby control (average) SR blur. 7. OPTIMAL
a: IMAGE-DEPENDENT
CASE
If the image energy spectrum is known or can be estimated, then a value of (Ycan be easily calculated which optimizes PCC to the frequency content of the image. To
PARAMETRIC
CUBIC
269
CONVOLUTION
establish this, Eqs. (13) (18), and (20) can be combined to yield e”(v) = e,(v) - 2ae,(v)
+ a2e2(v)
(244
where e,(v)
= 2 - 2P0(v) - gsin2(?iv)
e,(v) = i,(v)
(W
f gsin’(2nv)
(24~)
and e*(v) = &sin2(2nu)[l Figure 10 is a plot of ee(y), nonnegative for all v but e,(v) From Eqs. (17) and (24) it spectrum lg(v)12 the (average) /”
+ 6sin2(nv)].
(244
e,(y), and e2(v). As this figure illustrates, e2(v) is is not; below the Nyquist frequency e,(v) is negative. follows that for any band-limited image with energy SR blur is
e,( v)lg( v)12 dv - 2aj”’ e, (v)lg( v)12 dv + a”J” e2( v)lg( v)I’ dv. - Yc - “c - vc
(25)
Since Eq. (25) is quadratic in LYand the a2 coefficient is positive, the SR blur will be minimal in an image-dependent sense when
(y’(y
opt
1”’ e,(v)lh)12 dv = - “c
(26)
“’ e2( v)lg( v)12 dv . - “c
J
In other words, for an image with energy spectrum lg(v)12, the (average) SR blur -
e?(u)
= e,(v)
- zael(v)
+ cl%2(lJ) -\ 8’ \ \ \ el (u) ,I’ \ I \
0.3
\
0.2 '; 1, al
en(u) r;
0. I
c 0
\
\
/
I .o
/
-0.1
FIG. vertical
10. The three components axis for eo( v).)
of the function
e*(v)
for
PCC.
(Note
the different
scaling
of the
270
PARK
AND
SCHOWENGERDT
associated with parametric cubic convolution is least when OLsatisfies Eq. (26). Since e,(v) is negative below the Nyquist frequency, for most images a,rt will be negative; in particular, this must be true if the image is sufficiently sampled. However, if the image is undersampled and the image spectrum has significant energy above the Nyquist frequency, aopt may be nearly zero or even positive. To illustrate Eq. (26) for a realistic image energy spectrum, consider an image which is formed by scanning an edge with an ideal aperture. The resultant energy spectrum [7] is
li?w12 =
sinc2( sv) 4T2v2
(27)
where s is a dimensionless parameter which represents the width of the scanning aperture (IFOV) relative to the sample interval. Equivalently, s is the sample rate, i.e., the number of samples per IFOV. Equation (27) provides an important (one-dimensional) energy spectrum model for remote sensing applications since many images of interest are agricultural or urban and dominated by edges. For most applications the range 0.5 I s I 2.0 is realistic; at s = 1.0 consecutive IFOVs are just contiguous while s = 2.0 corresponds to a 50% overlap in consecutive IFOVs. Figure 11 illustrates (Y,-+~as a function of s. It is interesting to note that for this realistic range of sample rates, (Y,~*is negative, relatively constant, and less than the optimal image-independent value (Y= -0.5 but larger than the common value (Y= - 1.0. For example, aopt = -0.68, -0.70, and -0.61 at s = 1.0, 1.5, and 2.0, respectively. An average value of a,rt for this range of s is approximately - 3. For comparison purposes, Fig. 11 also illustrates (Y,~~for the image spectrum It(v
= sinc*(sv).
(28)
This spectrum corresponds to the image of a point source and so all scene energies ]f (v)] 2 are equal. For this extreme case, when the image is significantly undersampled (s < 0.9), (Y,,,~is actually positive. However, for typical sample rates (s 2: 1.5),
0.5 ‘, a0
0
-0.5
-1.0
FIG.
1 I. The
optimum
parameter
c
as a function
of s (samples/IFOV)
for two different
scenes.
PARAMETRIC
CUBIC CONVOLUTION
271
a opt is again negative and slightly less than (Y = - 0.5. Equations (27) and (28) define a range of energy spectra within which most remote sensing images can be expected to lie. 8. DISCUSSION
Cubic convolution is an image reconstruction algorithm which has enjoyed recent popularity, particularly in remote sensing applications where it is frequently necessary to use resampling to provide accurate radiometry in a coordinate system different from that of the original image. Unfortunately, the traditional implementation of cubic convolution has ignored its parametric nature and its associated ability to shape the energy spectrum of the reconstructed image. A parametric implementation of cubic convolution has been presented in this paper; the algorithm (PCC) is computationally efficient and the PCC parameter (Yprovides the ability to control the frequency characteristics of the reconstruction process. Standard (nonparametric) cubic convolution, defined by (Y= - 1, is widely used as a compromise approximation to the ideal sine interpolator; however, as this paper and [S] demonstrate, the (image-independent) reconstruction properties of this algorithm can be significantly improved by using instead (Y= -0.5. In particular, the choice (Y = -0.5 optimizes the low frequency suppression of sampling and reconstruction error and yields a nonparametric algorithm whose frequency characteristics are similar but superior to nearest neighbor and linear interpolation. Since PCC is parametric, it is natural to ask if (Ycan be optimized to the image. Specifically, although (Y= - 0.5 represents a good default value for general purpose implementation, if the image energy spectrum is known (or can be estimated), for what value of cwis the sampling and reconstruction error minimized? The general answer is provided by Eq. (26); in particular, for representative sampling rates and a scene dominated by edges, the optimal value of (Yis approximately - 3. An alternate approach to optimizing the choice of (Y is based upon visual appearance using a display generated with PCC and varying (Yinteractively. Such an interactive search can be implemented efficiently from Eq. [l l] because all PCC reconstructed images (corresponding to a given sampled image) are just the sum of two component images; one is a low-pass filtered version of the original image and the other is a high-pass filtered version that enhances high contrast features, i.e., edges, lines, etc., by an amount proportional to CL The reconstruction process would then be analogous to interactive high frequency enhancement [ 131. Finally, it should be mentioned that the optimization of (Y as discussed in the previous paragraph is based upon a global analysis and thus yields one value of LX which is used throughout the image. Many of the results of this paper could be extended to a local, adaptive reconstruction process in which (Ywould vary throughout the image. However, a locally adaptive implementation of PCC would be practical only if coupled with a very fast algorithm for estimating the local image frequency content; otherwise any increase in the image quality of the reconstructed image would be negated by a significant increase in computational expense.
ACKNOWLEDGMENT
Robert A. Schowengerdt was supported during this research by Air Force Office of Scientific Research Grant AFOSR-8 1-O170.
212
PARK
AND
SCHOWENGERDT
REFERENCES I. S. S. Rifman and D. M. McKinnon, Evaluation of digital correction techniques for ERTS images, TRW Report 20634-6003-m-02, July 1974. 2. F. C. Billingsley, Data preprocessing and processing, in Manual of Remote Sensing, 2nd ed., Chap. 16, American Society of Photogrammetry, to appear, 1983. 3. R. Bernstein. Digital image processing of earth observation sensor data, IBM J. Res. Deuel. 20, 1976. 40-57. 4. K. W. Simon, Digital image reconstruction and resampling of Landsat imagery, Symposium on Machine Processing of Remotely Sensed Data, Purdue University, Lafayette, Ind., June 1975, IEEE 75 CHlCKI9-O-C, pp. 3A-l-3A-11. 5. S. Shlien, Geometric correction, registration and resampling of Landsat imagery, Canad. J. Remote Sensing 5, 1979, 74-89. 6. P. F. Holkenbrink, Manual on Characteristics of La&at Computer-Compatible Tapes Produced bv the EROS Data Center Digital Image Processing System, U.S. Geological Survey EROS Data Center Report, 1978. 7. S. K. Park and R. A. Schowengerdt, Image sampling, reconstruction, and the effect of sample-scene phasing, J. Appl. Opt. 21, 1982. 8. R. G. Keys, Cubic convolution interpolation for digital image processing, IEEE Trans. Acoust. Speech Signal Process. ASSP-29, 198 1, 1153- 1160. 9. J. D. Gaskill, Linear Systems, Fourier Transforms, und Optics, Wiley, New York, 1978. IO. W. K. Pratt, DIgital Image Processing, Wiley, New York, 1978. I I. H. C. Andrews and B. R. Hunt, Digital Image Restoration, Prentice-Hall, Englewood Cliffs, N.J.. 1977. 12. C. E. Shannon, Communications in the presence of noise, Proc. IRE 37, 1949, 10-21. 13. J. -S. Lee, Digital image enhancement and noise filtering by use of local statistics. IEEE Trans. Pattern Anal. Mach. Intel. PAMI-2. 1980, 165-168.
VISION,
GRAPHICS,
AND
IMAGE
PROCESSING
23,
258-212
Image Reconstruction by Parametric Cubic Convolution STEPHEN K. PARK NASA,
Langley
Reseurch
Center,
Hampton,
Virginia
23465
AND ROBERT A. SCHOWENGERDT Office
of Arid
La&
Studies
and Electrical Tucson,
Engineering Department, Arizona 85721
University
of Arizona,
Received May 2 1. 1982; revised July 27, 1982 A parametric implementation of cubic convolution image reconstruction is presented which is generally superior to the standard algorithm and which can be optimized to the frequency content of the image. I. INTRODUCTION
Cubic convolution is a process originally developed for the reconstruction (i.e., interpolation, resampling) of Landsat digital images [ 11. Since its development, this process has come to be generally acknowledged as providing a good compromise between computational complexity and reconstruction accuracy [2]. Although it is recognized that cubic convolution image reconstruction is actually interpolation implemented by convolving the image samples with a smooth, spatially limited, one-parameter family of functions constructed by joining cubic polynomials in a splinelike fashion [3, 41, traditionally the parametric nature of cubic convolution has not been exploited and one member of this family has received virtually all the attention; it is the interpolation function originally proposed [l] and it remains the standard for comparison [5] and software implementation [6]. In this paper cubic convolution is formulated and analyzed parametrically, i.e., the reconstruction properties of the one-parameter family of cubic convolution interpolation functions are discussed and the image degradation (i.e., sampling and reconstruction blur [7]) associated with reasonable choices of this parameter is analyzed. It is demonstrated, with an analysis in the frequency domain, that in an image-independent sense there is an optimal value for this parameter and it is not the standard value commonly referenced in the literature. This important conclusion is consistent with a recently published result [8], which is based upon a complementary spatial domain analysis. In addition, it is demonstrated that in an image-dependent sense, cubic convolution can be tailored to any class of images characterized by a common energy spectrum, i.e., a different optimal parameter value exists for each different image energy spectrum. Thus, this paper establishes theoretically that a parametric implementation of cubic convolution image reconstruction is possible which is generally superior to the standard and which can be adapted to the frequency content of the image. 258 0734-189X/83
$3.00
Copyright 0 1983 by Academic Press, Inc. All rights of reproduction in any form reserved.
PARAMETRIC 2. IMAGE
CUBIC CONVOLUTION
259
RECONSTRUCTION
In a recent article [7], the authors analyzed the degradation caused by image sampling and interpolative reconstruction. The present paper uses the results of this analysis as well as the same notation and terminology; for simplicity the analysis is one-dimensional. A typical imaging, sampling and reconstruction system is illustrated in Fig. 1. The input scene (or object) is denoted f( x), the PSF of the imaging subsystem is denoted h(x), and the output of the imaging subsystem (prior to sampling) is the image g(x). Thus
g(x)=m*w
(1)
where * denotes convolution. As in [7], the spatial coordinate x is normalized in units of sample interval and the output of the sampling subsystem is the sampled image g,(x) given by g,(x)
= dx)comb(x)
(24
where comb(x)
=
F 8(x-n) n=-CC
(‘W
is the familiar sampling function [9]. Image reconstruction is commonly modeled as [lo]
g,(x) = g,(x)* 44
(3)
where the interpolation function r(x) is the reconstruction subsystem’s response to a sampled input which is one at x =0 and zero at x = n = rtl, t-2,... . The objective of image reconstruction is to choose T(X) so that the reconstructed image g,(x) accurately reproduces the output of the imaging subsystem g(x). Image reconstruction should not be confused with (linear) image restoration [I l] whose objective is to choose r(x) so that g,(x) accurately reproduces the input to the imaging subsystem f(x). Image reconstruction is an interpolation process in which the image samples (Eq. (2 a )) are reproduced exactly and thus, in general, r(0) = 1 r(n)
f(x) v SCENE
FIG.
IMAGING s!(x) SUBSYSTEM +h(x) IMAGE
= 0
n = +l,
+2, . . . .
SAMPLING g$) RECONSTRUCTION q,(x) a SUBSYSTEM SUBSYSTEH .comb(x) SAMPLED ?r(x) RECDNSTRUCTED - IMAGE I HAGF
1. Block diagram of a general imaging, sampling, and reconstruction system.
(4)
260
PARK
AND
SCHOWENGERDT
Equation (3) can be written equivalently as
d-4 = E s(n)rb - 4 n=-CC
(5)
which illustrates that, for a fixed x, the extent to which neighboring image samples contribute to the reconstructed image is determined by the spread of the interpolation function. For digital image reconstruction, computational considerations dictate that this spread should be small; typically, T(X) is identically zero for 1x1larger than two. It is well known [9-121 that if the image g(x) is band-limited and sufficiently sampled, then it can be reconstructed exactly by using the ideal interpolation function r(x)
sin( 7rx ) = sine(x) = 7rx .
In practice, exact digital reconstruction is impossible because the ideal interpolation function is not spatially limited and the sum in Eq. (5) must be truncated by using a truncated form of Eq. (6) as the interpolation function. However, truncating Eq. (6) to a practical width, say 1x1 I 2, introduces undesirable ripples in the spectrum (i.e., Fourier transform) of the resultant interpolator [5] because of the slope discontinuity at x = + 2. One method of suppressing these ripples is to construct the interpolation function as a smooth, spatially limited approximation to sine(x) on the interval - 2 5 x I 2 and eliminate the slope discontinuity at the end points; cubic convolution is an implementation of this method. 3. PARAMETRIC
CUBIC
CONVOLUTION
(PCC)
The interpolation function commonly associated with cubic convolution image reconstruction is the piecewise cubic polynomial illustrated in Fig. 2. By inspection this function (i) is even, (ii) is identically zero for all Ix] 2 2; (iii) is smooth, i.e., continuous with a continuous slope for all x, and (iv) satisfies Eq. (4). In addition, it can be verified that this function satisfies an important fifth property: (v) for UN x f
r(x-n)=
1.
n=-W Equation (7) is satisfied if and only if the image reconstruction exactly a constant (i.e., uniform) image.
FIG.
comparison,
2. The
interpolation function commonly associated with cubic the ideal interpolation function, sinc( x). is also illustrated.
process reproduces
convolution
((Y = - 1). car
PARAMETRIC
CUBIC
261
CONVOLUTION
It is known [3] that the function in Fig. 2 is just one member (a = - 1) of a one-parameter family of piecewise cubic polynomials given by (a + 2)1x13 - (a + 3)1x12 + 1
44 = (Ylx13 - 5a1x1* + thqxl-
1x1< 1 1 I Ix] < 2 otherwise.
4a
0
(8)
The parameter LYhas physical significance; it is the slope of T(X) at x = 1. The common value [3-61 is (r = - 1, which duplicates the slope of the ideal interpolator (Eq. (6)) at x = 1 (see Fig. 2); however, it was recently proven [8] that cx= -0.5 is superior to (Y= - 1 in that it provides an interpolator with better convergence properties. An intermediate choice, cx= -0.75, is also sometimes advocated; this choice results from forcing Eq. (8) to have second derivative continuity at x = 1. It can be shown directly that Eq. (8) satisfies conditions (i) through (v) for all (Y. Various members of this family of interpolators are illustrated in Fig. 3 correspondingtoa= l,O, - 0.5, - 1, and - 2. Since T(X) is an even function, only values for x 2 0 are illustrated. Note that for (Y = 0, T(X) is nonnegative and identically zero for 1x1 2 1. Equation (8) can be written equivalently as
r(x) = To(X)+ q(x) where
cw + Nlxl - l)*
Ix] < 1 otherwise
0
(9b)
and
dx)
=
14*(1x1- 1)
I4 < 1
(1x1 - l)(]x] - 2)* 0
1 I Ix] I 2 otherwise.
I.0
c
0.5
r x
0
-0.5 FIG.
3. Five
members
of the family
of PCC
interpolators.
(9c)
262
PARK
AND
SCHOWENGERDT
Equation (9a) illustrates clearly the linear dependence of this family of interpolators upon the parameter a. The function To(x), illustrated in Fig. 4, can be thought of as a smooth approximation to the triangular function associated with linear interpolation and r,(x), also shown in Fig. 4, is the smooth correction which must be added to r,,(x) to yield a function r(x) = To(x) + w,(x) whose slope at x = 1 is (Y. Equation (9) satisfies conditions (i) through (v) and, in particular, Eq. (7) becomes f
ro(x-n)=
(104
1
and fj t-,(x - n) = 0 n=-CC
for all x. Equations (5) and (9) define parametric reconstruction process, i.e.,
(10’4
cubic convolution
(PCC) as an image
where
h(x) = g,(x)*r0(x)= “=-ccl E b+hb - n>
(11’4
and
g,(x) = dx)*r,(x)
= it
dnhb
- 4.
Because the interpolation functions To(x) and rI(x) are spatially limited to two and four sample intervals, respectively, and because they are just piecewise cubic polynomials, PCC represents a computationally efficient algorithm. In addition, this algorithm gives the user the ability to control the visual appearance (or radiometric
FIG. 4. The two
componentsof
the family
of PCC interpolators.
PARAMETRIC
263
CUBIC CONVOLUTION
q,(x):
STANDARD CUBIC CONVOLUTION (a =-
FIG. 5. A one-dimensional image reconstructed from 31 samples for several values of a: 0, -0.5. - 1. Also shown is g,(x), the high frequency image component which adds local gradient information.
distribution) of the reconstructed image by changing the PCC parameter (Y. From Eq. (1 l), all PCC images for a given sampled image are linear combinations of the components g,(x) and g,(x). Thus an interactive search for the “best” reconstructed image can be implemented quite efficiently. Figure 5 illustrates PCC image reconstruction for a representative portion of a sampled (one-dimensional) image. The 31 image samples are interpolated by PCC for three values of (Y: 0, - 0.5, and - 1. The reconstructed image corresponding to (Y= 0 is just g,(x); the other two reconstructed images are formed from g,(x) by adding ag, (x). Of the three reconstructed images, the one for (r = - 0.5 is best; for (Y = - 1, edge overshoots occur at x = 23.5 and 26.5 and artificial stepping occurs for x between 5.0 and 11.0; for LY= 0, the slope of go(x) is zero at all image samples causing the appearance of artifacts, particularly in the range 5.0 < x < 11 .O. On the basis of Fig. 5 and the results of Ref. [8], we conjecture that for a large class of images, PCC with (Y= - 0.5 is superior to both standard cubic convolution ((Y = - 1) and to the smooth approximation to linear interpolation provided by (Y = 0. Later in this paper, an analysis in the frequency domain is presented which confirms this conjecture. 4. RECONSTRUCTION
FILTER
The interpretation of Eq. (3) in the frequency (Fourier) domain is that image reconstruction can be thought of as a process in which the sampled image spectrum is multiplied by a reconstruction filter i(v) where i(v)
= /”
-m
r(x)e-2nxvidx
(12)
is the Fourier transform of the interpolation function T(X). Since x is normalized in units of sample interval, the frequency coordinate v has units of cycles per sample
264
PARK
FIG.
6. Five
members
AND
SCHOWENGERDT
of the family
of PCC
reconstruction
filters
interval. In this coordinate system the Nyquist (or folding) frequency is 0.5. It can be shown that the Fourier transform of Eq. (9) is
i(v) = to(v) + “i,(V)
(134
where
+I(~)= -&*
[sinc*( y) - sinc(2v)]
(13’4
and
i,(v) = - 2
[3 sinc2(2v) - 2 sinc(2Y) - sinc(4v)l.
(13c)
(4’
Figure 6, which is the frequency domain counterpart of Fig. 3, illustrates i(v) for five values of (Y: 1, 0, - 0.5, - 1, and - 2. For all CX,i(v) is a low-pass filter. Figure 7, the frequency domain counterpart to Fig. 4, illustrates tO(v) and i,(v). Since all the functions in Eq. (13a) are even functions of v, only values of v 2 0 are illustrated.
FIG.
7. The
two components
of the family
of PCC
reconstruction
filters
PARAMETRIC
265
CUBIC CONVOLUTION
From Fig. 6 it can be seen that the value (Y = - 0.5 has particular significance; for a! < - 0.5, i(v) is a reconstruction filter which enhances low frequencies while for (Y > - 0.5, i(v) suppresses low frequencies. This effect can be examined in more detail by expanding i(v) in powers of (av) which yields 3(v) = 1 - $2,
+ l)(~v)*
+ A(1601 + l)(~v)~
+ ... .
From Eq. (14) i(0) = 1 and i(v) is either curved up or down at v = 0 depending upon whether (~CX+ 1) is negative or positive. The value (Y= - 0.5 is optimal in that it eliminates the (TV)* term in Eq. (14) and provides the best low frequency approximation to the ideal reconstruction filter, i.e., to the Fourier transform of Eq. (6) which is Iv( < 0.5 otherwise. Note, from Fig. 6, that the choice (Y= -0.5 also provides an excellent high frequency approximation to Eq. (15), i.e., for (Y = -0.5, i(v) (Eq. (13)) is effectively zero for all Iv] > 1.0. Thus the reconstruction filter corresponding to (Y= -0.5 is also particularly effective at suppressing aliasing. 5. SAMPLING
AND RECONSTRUCTION
BLUR
It is well known that for practical systems, image sampling and reconstruction inevitably produce some degradation in the reconstructed image [lo], i.e., g,(x) (Eq. (3)) is not an exact copy of g(x) (Eq. (1)). One measure of this degradation (mean square radiometric error) (16)
called sampling and reconstruction blur (SR blur) in reference [7], is particularly amenable to analysis. For example, it can be shown that: (i) if the image is band-limited, i.e., the image energy spectrum, ] g( v)] *, is ’ zero for all ]v( greater than a cut-off frequency vC; and (ii) if the image is sufficiently sampled, i.e., v, is less than or equal to the Nyquist frequency, 0.5 cycles per sample interval, then Eq. (16) reduces to E& = /“’ e2(v)lg(v)lz - “c
dv
where e*(v)
= II - i(v
+ C li(v - n)l*. nto
(17b)
If the band-limited image is not sufficiently sampled, then SR blur depends upon the sample-scene phase, i.e., the location of the sampling grid relative to the scene. However, Eq. (17) remains valid in the sense that it represents the average SR blur over the ensemble of all possible (random, equally likely) sample-scene phases.
266
PARK AND SCHOWENGERDT
For PCC image reconstruction, Eq. (17b) defines e*(v) in terms of the reconstruction filter, Eq. (13). As Eq. (17a) reveals, the cause of SR blur is the presence of significant image energy at frequencies where e*(v) is not zero. With reference to Fig. 1, just as the image subsystem can be interpreted as a presampling filter acting upon the scene energy spectrum ]f(v)l* to produce image blur, the reconstruction subsystem can be interpreted as a postsampling filter acting upon the image energy spectrum ]g( v)]* to produce SR blur. Note that e*(v) is not just ]I - i(v)]*; sampling causes the appearance of the side-band terms Cn tO ]i( v - n )I *. The importance of the parameter (Yis that it provides a means of shaping the e*(v) vs v curve and thereby controlling SR blur. 6. OPTIMAL
a: IMAGE-INDEPENDENT
CASE
Since the reconstruction filter (Eq. (13)) is real, Eq. (17b) can be simplified by omitting the magnitude notation and including i*(v) in the infinite series to yield e*(v) = 1 - 2?(v) +
f P(v - n). n=--M
(18)
Because the associated interpolation function, Eq. (9) is spatially limited to the interval Ix] 5 2, it can be shown [7] that the infinite series in Eq. (18) is equivalent to i*(v
f
- n) = (r * r)(O) + 2 i
n=-CC
(r * r)(n)cos(2mnv)
n=l
(19)
where (r * r)(n) represents r(x) * r(x) evaluated at x = n. The coefficients (r * r)(n) can be determined algebraically and combined with Eq. (19) to yield E
i*(v
- n) = 1 - J$sin*(nu)
- $$sin*(2nv)
n=-CC + gsin*(2*v)[l
+ 6sin*(nv)].
Since Eq. (20) is quadratic in (Y, e*(v) will be also. Note that e*(v) is an even function of v. For most digital images the preponderance of energy exists at low frequencies. This implies that in an image-independent sense the optimal choice of (Yis that one which yields the best low frequency approximation to the equation e*(v) = 0. This optimality criteria is reinforced by the observation that, for the ideal reconstruction filter (Eq. (15)), e*(v) is identically zero for Iv] < 0.5. If Eqs. (18) and (20) are expanded in powers of (TV) and combined, the resulting low frequency expansion is e*(v)
= -&(2a
+ l)‘(pv)*
+ &-(2,
+ 1)(14a! + 9)(7r~)~ + ... .
Equation (21) is a particularly satisfying result because it again emphasizes the (image-independent) optimality of (Y= - 0.5; if (2a + 1) = 0, then both the (TV)* and (a~)~ terms are eliminated and e*(v) = (constant)( rv)6
(22)
PARAMETRIC
CUBIC
CONVOLUTION
267
for low frequencies. (We have verified numerically that the constant in Eq. (22) is not zero.) In Ref. [7], cubic convolution (with (Y = - 1) was compared with the standard nearest neighbor and linear image reconstruction techniques. The results of this analysis may be used to derive low frequency expansions for nearest neighbor
e’(v)= &7v)2- &rv)4 +&+v)6 + . ..
(234
and for linear
f?‘(v)= &7v)” - $&rv)” + *.. . By comparing Eqs. (21), (22), (23a), and (23b), we see that for images dominated by low frequency energy, linear interpolation is superior to nearest neighbor and PCC is superior to linear but only if a = - 0.5. Figure 8 is a plot of e(y) (not e’(v)) vs v for four interpolators: nearest neighbor, linear, standard cubic convolution (a = - I), and PCC with (Y = -0.5. The curves for nearest neighbor, linear, and PCC with (Y = -0.5 are similar in shape but different in magnitude; the curve for standard cubic convolution is markedly different with a low frequency bump caused by the nonoptimality of (Y= - 1. Since
FIG.
convolution frequency
8. The function e(v) for four common interpolators, nearest neighbor, linear, standard cubic ((1 = - I), and PCC with a = - 0.5. The nonoptimality of a = - 1 is associated with the low bump in the curve. The Nyquist frequency is 0.5 cycles per sample interval.
268
PARK AND SCHOWENGERDT
I
10‘3 0
FIG.
I, 0.2
I 0.4
I I
I,
I
0.6
0.8
I
I I.0
9. The function e(v) for PCC with a = 1, 0, - 0.5, - I, and - 2.
e(v) for PCC with (Y = - 0.5 is less than either linear or nearest neighbor over the range Iv] I 0.5, we conclude that for any band-limited and sufficiently sampled image, parametric cubic convolution with a = - 0.5 yields less SR blur than either linear or nearest neighbor interpolation. The irony of this important theoretical result is that it is not also true for the standard value (Y = - 1, e.g., for this value, and for any image
with all its energy at frequencies below v = 0.12, linear interpolation would yield less SR blur. The insights provided by Eqs. (21), (22), (23a), and (23b) are consistent with the results of Ref. [8], based upon a local spatial domain analysis. In Ref. [8] it was proven that cubic convolution is best in that it is third-order convergent if (Y= - 0.5; linear interpolation is second-order convergent and thus is the intermediate choice; nearest neighbor is the worst because it is only first-order convergent. For (Y* - 0.5, cubic convolution is just first-order convergent. Figure 9 is a plot which illustrates e(v) for PCC corresponding to the five values of LY: 1, 0, - 0.5, - 1, and - 2. As this figure demonstrates, a parametric implementation of cubic convolution is desirable because (Ycan be used to shape e*(v) and thereby control (average) SR blur. 7. OPTIMAL
a: IMAGE-DEPENDENT
CASE
If the image energy spectrum is known or can be estimated, then a value of (Ycan be easily calculated which optimizes PCC to the frequency content of the image. To
PARAMETRIC
CUBIC
269
CONVOLUTION
establish this, Eqs. (13) (18), and (20) can be combined to yield e”(v) = e,(v) - 2ae,(v)
+ a2e2(v)
(244
where e,(v)
= 2 - 2P0(v) - gsin2(?iv)
e,(v) = i,(v)
(W
f gsin’(2nv)
(24~)
and e*(v) = &sin2(2nu)[l Figure 10 is a plot of ee(y), nonnegative for all v but e,(v) From Eqs. (17) and (24) it spectrum lg(v)12 the (average) /”
+ 6sin2(nv)].
(244
e,(y), and e2(v). As this figure illustrates, e2(v) is is not; below the Nyquist frequency e,(v) is negative. follows that for any band-limited image with energy SR blur is
e,( v)lg( v)12 dv - 2aj”’ e, (v)lg( v)12 dv + a”J” e2( v)lg( v)I’ dv. - Yc - “c - vc
(25)
Since Eq. (25) is quadratic in LYand the a2 coefficient is positive, the SR blur will be minimal in an image-dependent sense when
(y’(y
opt
1”’ e,(v)lh)12 dv = - “c
(26)
“’ e2( v)lg( v)12 dv . - “c
J
In other words, for an image with energy spectrum lg(v)12, the (average) SR blur -
e?(u)
= e,(v)
- zael(v)
+ cl%2(lJ) -\ 8’ \ \ \ el (u) ,I’ \ I \
0.3
\
0.2 '; 1, al
en(u) r;
0. I
c 0
\
\
/
I .o
/
-0.1
FIG. vertical
10. The three components axis for eo( v).)
of the function
e*(v)
for
PCC.
(Note
the different
scaling
of the
270
PARK
AND
SCHOWENGERDT
associated with parametric cubic convolution is least when OLsatisfies Eq. (26). Since e,(v) is negative below the Nyquist frequency, for most images a,rt will be negative; in particular, this must be true if the image is sufficiently sampled. However, if the image is undersampled and the image spectrum has significant energy above the Nyquist frequency, aopt may be nearly zero or even positive. To illustrate Eq. (26) for a realistic image energy spectrum, consider an image which is formed by scanning an edge with an ideal aperture. The resultant energy spectrum [7] is
li?w12 =
sinc2( sv) 4T2v2
(27)
where s is a dimensionless parameter which represents the width of the scanning aperture (IFOV) relative to the sample interval. Equivalently, s is the sample rate, i.e., the number of samples per IFOV. Equation (27) provides an important (one-dimensional) energy spectrum model for remote sensing applications since many images of interest are agricultural or urban and dominated by edges. For most applications the range 0.5 I s I 2.0 is realistic; at s = 1.0 consecutive IFOVs are just contiguous while s = 2.0 corresponds to a 50% overlap in consecutive IFOVs. Figure 11 illustrates (Y,-+~as a function of s. It is interesting to note that for this realistic range of sample rates, (Y,~*is negative, relatively constant, and less than the optimal image-independent value (Y= -0.5 but larger than the common value (Y= - 1.0. For example, aopt = -0.68, -0.70, and -0.61 at s = 1.0, 1.5, and 2.0, respectively. An average value of a,rt for this range of s is approximately - 3. For comparison purposes, Fig. 11 also illustrates (Y,~~for the image spectrum It(v
= sinc*(sv).
(28)
This spectrum corresponds to the image of a point source and so all scene energies ]f (v)] 2 are equal. For this extreme case, when the image is significantly undersampled (s < 0.9), (Y,,,~is actually positive. However, for typical sample rates (s 2: 1.5),
0.5 ‘, a0
0
-0.5
-1.0
FIG.
1 I. The
optimum
parameter
c
as a function
of s (samples/IFOV)
for two different
scenes.
PARAMETRIC
CUBIC CONVOLUTION
271
a opt is again negative and slightly less than (Y = - 0.5. Equations (27) and (28) define a range of energy spectra within which most remote sensing images can be expected to lie. 8. DISCUSSION
Cubic convolution is an image reconstruction algorithm which has enjoyed recent popularity, particularly in remote sensing applications where it is frequently necessary to use resampling to provide accurate radiometry in a coordinate system different from that of the original image. Unfortunately, the traditional implementation of cubic convolution has ignored its parametric nature and its associated ability to shape the energy spectrum of the reconstructed image. A parametric implementation of cubic convolution has been presented in this paper; the algorithm (PCC) is computationally efficient and the PCC parameter (Yprovides the ability to control the frequency characteristics of the reconstruction process. Standard (nonparametric) cubic convolution, defined by (Y= - 1, is widely used as a compromise approximation to the ideal sine interpolator; however, as this paper and [S] demonstrate, the (image-independent) reconstruction properties of this algorithm can be significantly improved by using instead (Y= -0.5. In particular, the choice (Y = -0.5 optimizes the low frequency suppression of sampling and reconstruction error and yields a nonparametric algorithm whose frequency characteristics are similar but superior to nearest neighbor and linear interpolation. Since PCC is parametric, it is natural to ask if (Ycan be optimized to the image. Specifically, although (Y= - 0.5 represents a good default value for general purpose implementation, if the image energy spectrum is known (or can be estimated), for what value of cwis the sampling and reconstruction error minimized? The general answer is provided by Eq. (26); in particular, for representative sampling rates and a scene dominated by edges, the optimal value of (Yis approximately - 3. An alternate approach to optimizing the choice of (Y is based upon visual appearance using a display generated with PCC and varying (Yinteractively. Such an interactive search can be implemented efficiently from Eq. [l l] because all PCC reconstructed images (corresponding to a given sampled image) are just the sum of two component images; one is a low-pass filtered version of the original image and the other is a high-pass filtered version that enhances high contrast features, i.e., edges, lines, etc., by an amount proportional to CL The reconstruction process would then be analogous to interactive high frequency enhancement [ 131. Finally, it should be mentioned that the optimization of (Y as discussed in the previous paragraph is based upon a global analysis and thus yields one value of LX which is used throughout the image. Many of the results of this paper could be extended to a local, adaptive reconstruction process in which (Ywould vary throughout the image. However, a locally adaptive implementation of PCC would be practical only if coupled with a very fast algorithm for estimating the local image frequency content; otherwise any increase in the image quality of the reconstructed image would be negated by a significant increase in computational expense.
ACKNOWLEDGMENT
Robert A. Schowengerdt was supported during this research by Air Force Office of Scientific Research Grant AFOSR-8 1-O170.
212
PARK
AND
SCHOWENGERDT
REFERENCES I. S. S. Rifman and D. M. McKinnon, Evaluation of digital correction techniques for ERTS images, TRW Report 20634-6003-m-02, July 1974. 2. F. C. Billingsley, Data preprocessing and processing, in Manual of Remote Sensing, 2nd ed., Chap. 16, American Society of Photogrammetry, to appear, 1983. 3. R. Bernstein. Digital image processing of earth observation sensor data, IBM J. Res. Deuel. 20, 1976. 40-57. 4. K. W. Simon, Digital image reconstruction and resampling of Landsat imagery, Symposium on Machine Processing of Remotely Sensed Data, Purdue University, Lafayette, Ind., June 1975, IEEE 75 CHlCKI9-O-C, pp. 3A-l-3A-11. 5. S. Shlien, Geometric correction, registration and resampling of Landsat imagery, Canad. J. Remote Sensing 5, 1979, 74-89. 6. P. F. Holkenbrink, Manual on Characteristics of La&at Computer-Compatible Tapes Produced bv the EROS Data Center Digital Image Processing System, U.S. Geological Survey EROS Data Center Report, 1978. 7. S. K. Park and R. A. Schowengerdt, Image sampling, reconstruction, and the effect of sample-scene phasing, J. Appl. Opt. 21, 1982. 8. R. G. Keys, Cubic convolution interpolation for digital image processing, IEEE Trans. Acoust. Speech Signal Process. ASSP-29, 198 1, 1153- 1160. 9. J. D. Gaskill, Linear Systems, Fourier Transforms, und Optics, Wiley, New York, 1978. IO. W. K. Pratt, DIgital Image Processing, Wiley, New York, 1978. I I. H. C. Andrews and B. R. Hunt, Digital Image Restoration, Prentice-Hall, Englewood Cliffs, N.J.. 1977. 12. C. E. Shannon, Communications in the presence of noise, Proc. IRE 37, 1949, 10-21. 13. J. -S. Lee, Digital image enhancement and noise filtering by use of local statistics. IEEE Trans. Pattern Anal. Mach. Intel. PAMI-2. 1980, 165-168.