Xử lý hình ảnh kỹ thuật số P7

Chia sẻ: Do Xon Xon | Ngày: | Loại File: PDF | Số trang:23

Thêm vào BST

Báo xấu

75
lượt xem 7
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

SUPERPOSITION AND CONVOLUTION In Chapter 1, superposition and convolution operations were derived for continuous two-dimensional image fields. This chapter provides a derivation of these operations for discrete two-dimensional images. Three types of superposition and convolution operators are defined: finite area, sampled image, and circulant area. The finite-area operator is a linear filtering process performed on a discrete image data array.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Xử lý hình ảnh kỹ thuật số P7

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) 7 SUPERPOSITION AND CONVOLUTION In Chapter 1, superposition and convolution operations were derived for continuous two-dimensional image fields. This chapter provides a derivation of these operations for discrete two-dimensional images. Three types of superposition and convolution operators are defined: finite area, sampled image, and circulant area. The finite-area operator is a linear filtering process performed on a discrete image data array. The sampled image operator is a discrete model of a continuous two-dimensional image filtering process. The circulant area operator provides a basis for a computationally efficient means of performing either finite-area or sampled image superposition and convolution. 7.1. FINITE-AREA SUPERPOSITION AND CONVOLUTION Mathematical expressions for finite-area superposition and convolution are devel- oped below for both series and vector-space formulations. 7.1.1. Finite-Area Superposition and Convolution: Series Formulation Let F ( n1, n 2 ) denote an image array for n1, n2 = 1, 2,..., N. For notational simplicity, all arrays in this chapter are assumed square. In correspondence with Eq. 1.2-6, the image array can be represented at some point ( m 1 , m 2 ) as a sum of amplitude weighted Dirac delta functions by the discrete sifting summation F ( m 1, m 2 ) = ∑∑ F ( n 1, n 2 )δ ( m 1 – n 1 + 1, m 2 – n 2 + 1 ) (7.1-1) n1 n2 161
162 SUPERPOSITION AND CONVOLUTION The term 1 if m1 = n 1 and m 2 = n 2 (7.1-2a)  δ ( m 1 – n 1 + 1, m 2 – n 2 + 1 ) =   0 otherwise (7.1-2b) is a discrete delta function. Now consider a spatial linear operator O { · } that pro- duces an output image array Q ( m 1, m 2 ) = O { F ( m 1, m 2 ) } (7.1-3) by a linear spatial combination of pixels within a neighborhood of ( m 1, m 2 ) . From the sifting summation of Eq. 7.1-1,   Q ( m 1, m 2 ) = O  ∑ ∑ F ( n 1, n 2 )δ ( m1 – n 1 + 1, m 2 – n 2 + 1 ) (7.1-4a)  n1 n2  or Q ( m 1, m 2 ) = ∑∑ F ( n 1, n 2 )O { δ ( m 1 – n 1 + 1, m 2 – n 2 + 1 ) } (7.1-4b) n1 n2 recognizing that O { · } is a linear operator and that F ( n 1, n 2 ) in the summation of Eq. 7.1-4a is a constant in the sense that it does not depend on ( m 1, m 2 ) . The term O { δ ( t 1, t 2 ) } for ti = m i – n i + 1 is the response at output coordinate ( m 1, m 2 ) to a unit amplitude input at coordinate ( n 1, n 2 ) . It is called the impulse response function array of the linear operator and is written as δ ( m 1 – n 1 + 1, m 2 – n 2 + 1 ; m 1, m 2 ) = O { δ ( t1, t2 ) } for 1 ≤ t1, t2 ≤ L (7.1-5) and is zero otherwise. For notational simplicity, the impulse response array is con- sidered to be square. In Eq. 7.1-5 it is assumed that the impulse response array is of limited spatial extent. This means that an output image pixel is influenced by input image pixels only within some finite area L × L neighborhood of the corresponding output image pixel. The output coordinates ( m 1, m 2 ) in Eq. 7.1-5 following the semicolon indicate that in the general case, called finite area superposition, the impulse response array can change form for each point ( m 1, m 2 ) in the processed array Q ( m 1, m 2 ). Follow- ing this nomenclature, the finite area superposition operation is defined as
FINITE-AREA SUPERPOSITION AND CONVOLUTION 163 FIGURE 7.1-1. Relationships between input data, output data, and impulse response arrays for finite-area superposition; upper left corner justified array definition. Q ( m 1, m 2 ) = ∑∑ F ( n 1, n 2 )H ( m 1 – n 1 + 1, m 2 – n 2 + 1 ; m 1, m 2 ) (7.1-6) n 1 n2 The limits of the summation are MAX { 1, m i – L + 1 } ≤ n i ≤ MIN { N, m i } (7.1-7) where MAX { a, b } and MIN { a, b } denote the maximum and minimum of the argu- ments, respectively. Examination of the indices of the impulse response array at its extreme positions indicates that M = N + L - 1, and hence the processed output array Q is of larger dimension than the input array F. Figure 7.1-1 illustrates the geometry of finite-area superposition. If the impulse response array H is spatially invariant, the superposition operation reduces to the convolution operation. Q ( m 1, m 2 ) = ∑∑ F ( n 1, n 2 )H ( m 1 – n 1 + 1, m2 – n 2 + 1 ) (7.1-8) n1 n2 Figure 7.1-2 presents a graphical example of convolution with a 3 × 3 impulse response array. Equation 7.1-6 expresses the finite-area superposition operation in left-justified form in which the input and output arrays are aligned at their upper left corners. It is often notationally convenient to utilize a definition in which the output array is cen- tered with respect to the input array. This definition of centered superposition is given by
164 SUPERPOSITION AND CONVOLUTION FIGURE 7.1-2. Graphical example of finite-area convolution with a 3 × 3 impulse response array; upper left corner justified array definition. Q c ( j 1, j2 ) = ∑∑ F ( n 1, n 2 )H ( j 1 – n 1 + L c, j2 – n 2 + L c ; j1, j 2 ) (7.1-9) n 1 n2 where – ( L – 3 ) ⁄ 2 ≤ ji ≤ N + ( L – 1 ) ⁄ 2 and L c = ( L + 1 ) ⁄ 2 . The limits of the summa- tion are MAX { 1, j i – ( L – 1 ) ⁄ 2 } ≤ n i ≤ MIN { N, j i + ( L – 1 ) ⁄ 2 } (7.1-10) Figure 7.1-3 shows the spatial relationships between the arrays F, H, and Qc for cen- tered superposition with a 5 × 5 impulse response array. In digital computers and digital image processors, it is often convenient to restrict the input and output arrays to be of the same dimension. For such systems, Eq. 7.1-9 needs only to be evaluated over the range 1 ≤ j i ≤ N . When the impulse response
FINITE-AREA SUPERPOSITION AND CONVOLUTION 165 FIGURE 7.1-3. Relationships between input data, output data, and impulse response arrays for finite-area superposition; centered array definition. array is located on the border of the input array, the product computation of Eq. 7.1-9 does not involve all of the elements of the impulse response array. This situa- tion is illustrated in Figure 7.1-3, where the impulse response array is in the upper left corner of the input array. The input array pixels “missing” from the computation are shown crosshatched in Figure 7.1-3. Several methods have been proposed to deal with this border effect. One method is to perform the computation of all of the impulse response elements as if the missing pixels are of some constant value. If the constant value is zero, the result is called centered, zero padded superposition. A variant of this method is to regard the missing pixels to be mirror images of the input array pixels, as indicated in the lower left corner of Figure 7.1-3. In this case the centered, reflected boundary superposition definition becomes Q c ( j 1, j2 ) = ∑∑ F ( n′ , n′ )H ( j 1 – n 1 + L c , j 2 – n 2 + L c ; j 1, j 2 ) 1 2 (7.1-11) n 1 n2 where the summation limits are ji – ( L – 1 ) ⁄ 2 ≤ ni ≤ ji + ( L – 1 ) ⁄ 2 (7.1-12)
166 SUPERPOSITION AND CONVOLUTION and  2 – ni for n i ≤ 0 (7.1-13a)    n' i =  n i for 1 ≤ n i ≤ N (7.1-13b)    2N – n  i for n i > N (7.1-13c) In many implementations, the superposition computation is limited to the range ( L + 1 ) ⁄ 2 ≤ j i ≤ N – ( L – 1 ) ⁄ 2 , and the border elements of the N × N array Qc are set to zero. In effect, the superposition operation is computed only when the impulse response array is fully embedded within the confines of the input array. This region is described by the dashed lines in Figure 7.1-3. This form of superposition is called centered, zero boundary superposition. If the impulse response array H is spatially invariant, the centered definition for convolution becomes Q c ( j 1, j 2 ) = ∑∑ F ( n 1, n 2 )H ( j 1 – n 1 + L c, j2 – n 2 + L c ) (7.1-14) n1 n2 The 3 × 3 impulse response array, which is called a small generating kernel (SGK), is fundamental to many image processing algorithms (1). When the SGK is totally embedded within the input data array, the general term of the centered convolution operation can be expressed explicitly as Q c ( j1, j 2 ) = H ( 3, 3 )F ( j1 – 1, j 2 – 1 ) + H ( 3, 2 )F ( j 1 – 1, j2 ) + H ( 3, 1 )F ( j 1 – 1, j 2 + 1 ) + H ( 2, 3 )F ( j 1, j2 – 1 ) + H ( 2, 2 )F ( j1, j 2 ) + H ( 2, 1 )F ( j 1, j 2 + 1 ) + H ( 1, 3 )F ( j 1 + 1, j 2 – 1 ) + H ( 1, 2 )F ( j 1 + 1, j2 ) + H ( 1, 1 )F ( j 1 + 1, j 2 + 1 ) (7.1-15) for 2 ≤ j i ≤ N – 1 . In Chapter 9 it will be shown that convolution with arbitrary-size impulse response arrays can be achieved by sequential convolutions with SGKs. The four different forms of superposition and convolution are each useful in var- ious image processing applications. The upper left corner–justified definition is appropriate for computing the correlation function between two images. The cen- tered, zero padded and centered, reflected boundary definitions are generally employed for image enhancement filtering. Finally, the centered, zero boundary def- inition is used for the computation of spatial derivatives in edge detection. In this application, the derivatives are not meaningful in the border region.
FINITE-AREA SUPERPOSITION AND CONVOLUTION 167 0.040 0.080 0.120 0.160 0.200 0.200 0.200 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.080 0.160 0.240 0.320 0.400 0.400 0.400 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.120 0.240 0.360 0.480 0.600 0.600 0.600 0.000 0.000 1.000 1.000 1.000 1.000 1.000 0.160 0.320 0.480 0.640 0.800 0.800 0.800 0.000 0.000 1.000 1.000 1.000 1.000 1.000 0.200 0.400 0.600 0.800 1.000 1.000 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 0.200 0.400 0.600 0.800 1.000 1.000 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 0.200 0.400 0.600 0.800 1.000 1.000 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 (a) Upper left corner justified (b) Centered, zero boundary 0.360 0.480 0.600 0.600 0.600 0.600 0.600 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.480 0.640 0.800 0.800 0.800 0.800 0.800 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.600 0.800 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.600 0.800 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.600 0.800 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.600 0.800 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.600 0.800 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 (c) Centered, zero padded (d) Centered, reflected FIGURE 7.1-4 Finite-area convolution boundary conditions, upper left corner of convolved image. Figure 7.1-4 shows computer printouts of pixels in the upper left corner of a convolved image for the four types of convolution boundary conditions. In this example, the source image is constant of maximum value 1.0. The convolution impulse response array is a 5 × 5 uniform array. 7.1.2. Finite-Area Superposition and Convolution: Vector-Space Formulation 2 If the arrays F and Q of Eq. 7.1-6 are represented in vector form by the N × 1 vec- 2 tor f and the M × 1 vector q, respectively, the finite-area superposition operation can be written as (2) q = Df (7.1-16) 2 2 where D is a M × N matrix containing the elements of the impulse response. It is convenient to partition the superposition operator matrix D into submatrices of dimension M × N . Observing the summation limits of Eq. 7.1-7, it is seen that D 1, 1 0 … 0 D 2, 1 D 2, 2 … 0 … … D = D L, 1 D L, 2 D M – L + 1, N (7.1-17) 0 D L + 1, 1 … … … … 0 … 0 D M, N
168 SUPERPOSITION AND CONVOLUTION 11 0 0 0 21 11 0 0 31 21 0 0 0 31 0 0 12 0 11 0 22 12 21 11 32 22 31 21 11 12 13 0 32 0 31 H = 21 22 23 D= 13 0 12 0 31 32 33 23 13 22 12 33 23 32 22 0 33 0 32 0 0 13 0 0 0 23 13 0 0 33 23 0 0 0 33 (a) (b) FIGURE 7.1-5 Finite-area convolution operators: (a) general impulse array, M = 4, N = 2, L = 3; (b) Gaussian-shaped impulse array, M = 16, N = 8, L = 9. The general nonzero term of D is then given by Dm ( m 1, n 1 ) = H ( m 1 – n 1 + 1, m 2 – n 2 + 1 ; m1, m 2 ) (7.1-18) 2 ,n 2 Thus, it is observed that D is highly structured and quite sparse, with the center band of submatrices containing stripes of zero-valued elements. If the impulse response is position invariant, the structure of D does not depend explicitly on the output array coordinate ( m 1, m2 ) . Also, Dm = Dm + 1 ,n 2 + 1 (7.1-19) 2 ,n 2 2 As a result, the columns of D are shifted versions of the first column. Under these conditions, the finite-area superposition operator is known as the finite-area convo- lution operator. Figure 7.1-5a contains a notational example of the finite-area con- volution operator for a 2 × 2 (N = 2) input data array, a 4 × 4 (M = 4) output data array, and a 3 × 3 (L = 3) impulse response array. The integer pairs (i, j) at each ele- ment of D represent the element (i, j) of H ( i, j ). The basic structure of D can be seen more clearly in the larger matrix depicted in Figure 7.l-5b. In this example, M = 16,
FINITE-AREA SUPERPOSITION AND CONVOLUTION 169 N = 8, L = 9, and the impulse response has a symmetrical Gaussian shape. Note that D is a 256 × 64 matrix in this example. Following the same technique as that leading to Eq. 5.4-7, the matrix form of the superposition operator may be written as M N T Q = ∑ ∑ Dm ,n Fvn um (7.1-20) m =1 n =1 If the impulse response is spatially invariant and is of separable form such that T H = h C hR (7.1-21) where hR and hC are column vectors representing row and column impulse responses, respectively, then D = DC ⊗ DR (7.1-22) The matrices D R and DC are M × N matrices of the form hR ( 1 ) 0 … 0 hR ( 2 ) hR ( 1 ) … hR ( 3 ) hR ( 2 ) … 0 hR ( 1 ) (7.1-23) … DR = hR ( L ) … … … 0 … 0 … 0 hR ( L ) The two-dimensional convolution operation may then be computed by sequential row and column one-dimensional convolutions. Thus T Q = D C FD R (7.1-24) In vector form, the general finite-area superposition or convolution operator requires 2 2 N L operations if the zero-valued multiplications of D are avoided. The separable operator of Eq. 7.1-24 can be computed with only NL ( M + N ) operations.
170 SUPERPOSITION AND CONVOLUTION 7.2. SAMPLED IMAGE SUPERPOSITION AND CONVOLUTION Many applications in image processing require a discretization of the superposition integral relating the input and output continuous fields of a linear system. For exam- ple, image blurring by an optical system, sampling with a finite-area aperture or imaging through atmospheric turbulence, may be modeled by the superposition inte- gral equation ˜ ∞ ∞ ˜ ˜ G ( x, y ) = ∫–∞ ∫–∞ F ( α, β )J ( x, y ; α, β ) dα dβ (7.2-1a) ˜ ˜ where F ( x ,y ) and G ( x, y ) denote the input and output fields of a linear system, ˜ respectively, and the kernel J ( x, y ; α , β ) represents the impulse response of the linear system model. In this chapter, a tilde over a variable indicates that the spatial indices of the variable are bipolar; that is, they range from negative to positive spatial limits. In this formulation, the impulse response may change form as a function of its four indices: the input and output coordinates. If the linear system is space invariant, the output image field may be described by the convolution integral ˜ ∞ ∞ ˜ ˜ G ( x, y ) = ∫–∞ ∫–∞ F ( α, β )J ( x – α, y – β ) dα dβ (7.2-1b) For discrete processing, physical image sampling will be performed on the output image field. Numerical representation of the integral must also be performed in order to relate the physical samples of the output field to points on the input field. Numerical representation of a superposition or convolution integral is an impor- tant topic because improper representations may lead to gross modeling errors or numerical instability in an image processing application. Also, selection of a numer- ical representation algorithm usually has a significant impact on digital processing computational requirements. As a first step in the discretization of the superposition integral, the output image field is physically sampled by a ( 2J + 1 ) × ( 2J + 1 ) array of Dirac pulses at a resolu- tion ∆S to obtain an array whose general term is ˜ ˜ G ( j 1 ∆S, j 2 ∆S) = G ( x, y )δ ( x – j1 ∆S, y – j2 ∆S) (7.2-2) where – J ≤ ji ≤ J . Equal horizontal and vertical spacing of sample pulses is assumed for notational simplicity. The effect of finite area sample pulses can easily be incor- ˜ porated by replacing the impulse response with J ( x, y ; α, β ) * P ( – x ,– y ), where P ( – x ,– y ) represents the pulse shape of the sampling pulse. The delta function may be brought under the integral sign of the superposition integral of Eq. 7.2-la to give ˜ ∞ ∞ ˜ ˜ G ( j 1 ∆S, j 2 ∆ S) = ∫–∞ ∫–∞ F ( α, β )J ( j1 ∆S, j 2 ∆S ; α, β ) dα dβ (7.2-3)
SAMPLED IMAGE SUPERPOSITION AND CONVOLUTION 171 It should be noted that the physical sampling is performed on the observed image spatial variables (x, y); physical sampling does not affect the dummy variables of integration ( α, β ) . Next, the impulse response must be truncated to some spatial bounds. Thus, let ˜ J ( x, y ; α, β ) = 0 (7.2-4) for x > T and y > T . Then, j 1 ∆S + T j 2 ∆S + T ˜ ˜ ˜ G ( j 1 ∆S, j1 ∆S) = ∫j ∆S – T ∫ j ∆S – T 1 2 F ( α, β ) J ( j 1 ∆S, j 2 ∆S ; α, β ) dα dβ (7.2-5) Truncation of the impulse response is equivalent to multiplying the impulse response by a window function V(x, y), which is unity for x < T and y < T and zero elsewhere. By the Fourier convolution theorem, the Fourier spectrum of G(x, y) is equivalently convolved with the Fourier transform of V(x, y), which is a two- dimensional sinc function. This distortion of the Fourier spectrum of G(x, y) results in the introduction of high-spatial-frequency artifacts (a Gibbs phenomenon) at spa- tial frequency multiples of 2π ⁄ T . Truncation distortion can be reduced by using a shaped window, such as the Bartlett, Blackman, Hamming, or Hanning windows (3), which smooth the sharp cutoff effects of a rectangular window. This step is especially important for image restoration modeling because ill-conditioning of the superposition operator may lead to severe amplification of the truncation artifacts. In the next step of the discrete representation, the continuous ideal image array ˜ F ( α, β ) is represented by mesh points on a rectangular grid of resolution ∆I and dimension ( 2K + 1 ) × ( 2K + 1 ). This is not a physical sampling process, but merely an abstract numerical representation whose general term is described by ˜ ˜ F ( k 1 ∆I, k2 ∆I) = F ( α, β )δ ( α – k 1 ∆I, α – k 2 ∆I) (7.2-6) where K iL ≤ k i ≤ K iU , with K iU and K iL denoting the upper and lower index limits. If the ultimate objective is to estimate the continuous ideal image field by pro- cessing the physical observation samples, the mesh spacing ∆I should be fine enough to satisfy the Nyquist criterion for the ideal image. That is, if the spectrum of the ideal image is bandlimited and the limits are known, the mesh spacing should be set at the corresponding Nyquist spacing. Ideally, this will permit perfect interpola- ˜ ˜ tion of the estimated points F ( k 1 ∆I, k 2 ∆I) to reconstruct F ( x, y ). The continuous integration of Eq. 7.2-5 can now be approximated by a discrete summation by employing a quadrature integration formula (4). The physical image samples may then be expressed as K 1U K 2U ˜ ˜ ˜ ˜ G ( j 1 ∆ S, j 2 ∆S) = ∑ ∑ F ( k 1 ∆ I, k 2 ∆ I)W ( k 1, k 2 )J ( j 1 ∆ S, j 2 ∆ S ; k 1 ∆I, k 2 ∆I ) k 1 = K 1L k 2 = K 2L (7.2-7)
172 SUPERPOSITION AND CONVOLUTION ˜ where W ( k 1 , k 2 ) is a weighting coefficient for the particular quadrature formula employed. Usually, a rectangular quadrature formula is used, and the weighting coefficients are unity. In any case, it is notationally convenient to lump the weight- ing coefficient and the impulse response function together so that ˜ ˜ ˜ H ( j1 ∆S, j 2 ∆S; k 1 ∆I, k 2 ∆I) = W ( k 1, k 2 )J ( j 1 ∆S, j 2 ∆S ; k 1 ∆I, k 2 ∆I) (7.2-8) Then, K 1U K 2U ˜ ˜ ˜ G ( j 1 ∆S, j 2 ∆S) = ∑ ∑ F ( k 1 ∆I, k 2 ∆I )H ( j 1 ∆S, j 2 ∆S ; k 1 ∆I, k 2 ∆I ) (7.2-9) k 1 = K 1L k 2 = K 2L ˜ Again, it should be noted that H is not spatially discretized; the function is simply evaluated at its appropriate spatial argument. The limits of summation of Eq. 7.2-9 are ∆S T ∆S T K iL = j i ------ – ----- - - K iU = ji ------ + ----- - - (7.2-10) ∆I ∆I N ∆I ∆I N where [ · ] N denotes the nearest integer value of the argument. Figure 7.2-1 provides an example relating actual physical sample values ˜ ˜ G ( j1 ∆ S, j2 ∆S) to mesh points F ( k 1 ∆I, k 2 ∆I ) on the ideal image field. In this exam- ple, the mesh spacing is twice as large as the physical sample spacing. In the figure, FIGURE 7.2-1. Relationship of physical image samples to mesh points on an ideal image field for numerical representation of a superposition integral.
SAMPLED IMAGE SUPERPOSITION AND CONVOLUTION 173 FIGURE 7.2-2. Relationship between regions of physical samples and mesh points for numerical representation of a superposition integral. the values of the impulse response function that are utilized in the summation of Eq. 7.2-9 are represented as dots. An important observation should be made about the discrete model of Eq. 7.2-9 for a sampled superposition integral; the physical area of the ideal image field ˜ F ( x, y ) containing mesh points contributing to physical image samples is larger ˜ than the sample image G ( j 1 ∆S, j2 ∆S) regardless of the relative number of physical samples and mesh points. The dimensions of the two image fields, as shown in Figure 7.2-2, are related by J ∆S + T = K ∆I (7.2-11) to within an accuracy of one sample spacing. At this point in the discussion, a discrete and finite model for the sampled super- ˜ position integral has been obtained in which the physical samples G ( j1 ∆S, j2 ∆S) are related to points on an ideal image field F 1˜ ( k ∆I, k ∆I) by a discrete mathemati- 2 cal superposition operation. This discrete superposition is an approximation to con- tinuous superposition because of the truncation of the impulse response function ˜ J ( x, y ; α, β ) and quadrature integration. The truncation approximation can, of course, be made arbitrarily small by extending the bounds of definition of the impulse response, but at the expense of large dimensionality. Also, the quadrature integration approximation can be improved by use of complicated formulas of quadrature, but again the price paid is computational complexity. It should be noted, however, that discrete superposition is a perfect approximation to continuous super- position if the spatial functions of Eq. 7.2-1 are all bandlimited and the physical
174 SUPERPOSITION AND CONVOLUTION sampling and numerical representation periods are selected to be the corresponding Nyquist period (5). It is often convenient to reformulate Eq. 7.2-9 into vector-space form. Toward ˜ ˜ this end, the arrays G and F are reindexed to M × M and N × N arrays, respectively, such that all indices are positive. Let ˜ F ( n 1 ∆I, n2 ∆I ) = F ( k 1 ∆I, k2 ∆I ) (7.2-12a) where n i = k i + K + 1 and let G ( m 1 ∆S, m 2 ∆S) = G ( j 1 ∆S, j2 ∆S ) (7.2-12b) where m i = j i + J + 1 . Also, let the impulse response be redefined such that ˜ H ( m1 ∆S, m 2 ∆S ; n 1 ∆ I, n 2 ∆I ) = H ( j1 ∆S, j2 ∆S ; k 1 ∆I, k 2 ∆I ) (7.2-12c) Figure 7.2-3 illustrates the geometrical relationship between these functions. The discrete superposition relationship of Eq. 7.2-9 for the shifted arrays becomes N 1U N 2U G ( m 1 ∆S, m 2 ∆S) = ∑ ∑ F ( n 1 ∆I, n 2 ∆I ) H ( m1 ∆ S, m 2 ∆S ; n 1 ∆I, n 2 ∆I ) n 1 = N 1L n 2 = N 2L (7.2-13) for ( 1 ≤ m i ≤ M ) where N iL = m i ∆S ------ - N iU = m i ∆S + 2T ------ ----- - - ∆I N ∆I ∆I N Following the techniques outlined in Chapter 5, the vectors g and f may be formed by column scanning the matrices G and F to obtain g = Bf (7.2-14) 2 2 where B is a M × N matrix of the form B 1, 1 B1, 2 … B ( 1, L ) 0… 0 … … … 0 B 2, 2 B = (7.2-15) 0 … 0 … 0 B M, N – L + 1 … B M, N
SAMPLED IMAGE SUPERPOSITION AND CONVOLUTION 175 FIGURE 7.2-3. Sampled image arrays. The general term of B is defined as Bm 2, n 2 ( m 1, n 1 ) = H ( m1 ∆S, m 2 ∆S ; n 1 ∆I, n 2 ∆I ) (7.2-16) for 1 ≤ m i ≤ M and m i ≤ n i ≤ m i + L – 1 where L = [ 2T ⁄ ∆I ] N represents the nearest odd integer dimension of the impulse response in resolution units ∆I . For descrip- tional simplicity, B is called the blur matrix of the superposition integral. If the impulse response function is translation invariant such that H ( m 1 ∆S, m 2 ∆S ; n 1 ∆I, n 2 ∆I ) = H ( m 1 ∆S – n 1 ∆I, m2 ∆S – n 2 ∆I ) (7.2-17) then the discrete superposition operation of Eq. 7.2-13 becomes a discrete convolu- tion operation of the form N 1U N 2U G ( m 1 ∆S, m 2 ∆S ) = ∑ ∑ F ( n 1 ∆I, n 2 ∆I )H ( m1 ∆S – n 1 ∆I, m 2 ∆S – n 2 ∆I ) n 1 = N 1L n 2 = N 2L (7.2-18) If the physical sample and quadrature mesh spacings are equal, the general term of the blur matrix assumes the form Bm 2, n 2 ( m 1, n 1 ) = H ( m 1 – n 1 + L, m 2 – n 2 + L ) (7.2-19)
176 SUPERPOSITION AND CONVOLUTION 11 12 13 H = 21 22 23 31 32 33 33 23 13 0 32 22 12 0 31 21 11 0 0 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 0 0 0 0 B= 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 0 0 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 (a) (b) FIGURE 7.2-4. Sampled infinite area convolution operators: (a) General impulse array, M = 2, N = 4, L = 3; (b) Gaussian-shaped impulse array, M = 8, N = 16, L = 9. In Eq. 7.2-19, the mesh spacing variable ∆I is understood. In addition, Bm = Bm (7.2-20) 2, n2 2 + 1, n 2 + 1 Consequently, the rows of B are shifted versions of the first row. The operator B then becomes a sampled infinite area convolution operator, and the series form rep- resentation of Eq. 7.2-19 reduces to m1 + L – 1 m 2 + L – 1 G ( m 1 ∆S, m 2 ∆S ) = ∑ ∑ F ( n 1, n 2 )H ( m 1 – n 1 + L, m2 – n 2 + L ) (7.2-21) n1 = m1 n2 = m2 where the sampling spacing is understood. Figure 7.2-4a is a notational example of the sampled image convolution operator for a 4 × 4 (N = 4) data array, a 2 × 2 (M = 2) filtered data array, and a 3 × 3 (L = 3) impulse response array. An extension to larger dimension is shown in Figure 7.2-4b for M = 8, N = 16, L = 9 and a Gaussian-shaped impulse response. When the impulse response is spatially invariant and orthogonally separable, B = BC ⊗ BR (7.2-22) where B R and B C are M × N matrices of the form
CIRCULANT SUPERPOSITION AND CONVOLUTION 177 hR ( L ) hR ( L – 1 ) … hR ( 1 ) 0 … 0 0 hR ( L ) … BR = (7.2-23) 0 … 0 … 0 hR ( L ) … hR ( 1 ) The two-dimensional convolution operation then reduces to sequential row and col- umn convolutions of the matrix form of the image array. Thus T G = B C FB R (7.2-24) 2 2 The superposition or convolution operator expressed in vector form requires M L operations if the zero multiplications of B are avoided. A separable convolution operator can be computed in matrix form with only ML ( M + N ) operations. 7.3. CIRCULANT SUPERPOSITION AND CONVOLUTION In circulant superposition (2), the input data, the processed output, and the impulse response arrays are all assumed spatially periodic over a common period. To unify the presentation, these arrays will be defined in terms of the spatially limited arrays considered previously. First, let the N × N data array F ( n1 ,n 2 ) be embedded in the upper left corner of a J × J array ( J > N ) of zeros, giving  F ( n 1, n 2 ) for 1 ≤ n i ≤ N (7.3-1a)   FE ( n 1, n 2 ) =     0 for N + 1 ≤ n i ≤ J (7.3-1b) In a similar manner, an extended impulse response array is created by embedding the spatially limited impulse array in a J × J matrix of zeros. Thus, let H(l , l ; m , m ) for 1 ≤ li ≤ L (7.3-2a)  1 2 1 2  H E ( l 1, l 2 ; m 1, m 2 ) =    0 for L + 1 ≤ l i ≤ J (7.3-2b) 
178 SUPERPOSITION AND CONVOLUTION Periodic arrays F E ( n 1, n 2 ) and H E ( l 1, l 2 ; m 1, m 2 ) are now formed by replicating the extended arrays over the spatial period J. Then, the circulant superposition of these functions is defined as J J KE ( m 1, m 2 ) = ∑ ∑ F E ( n 1, n 2 )H E ( m 1 – n 1 + 1, m 2 – n 2 + 1 ; m 1, m 2) (7.3-3) n1 = 1 n2 = 1 Similarity of this equation with Eq. 7.1-6 describing finite-area superposition is evi- dent. In fact, if J is chosen such that J = N + L – 1, the terms F E ( n 1, n 2 ) = F ( n 1, n 2 ) for 1 ≤ ni ≤ N . The similarity of the circulant superposition operation and the sam- pled image superposition operation should also be noted. These relations become clearer in the vector-space representation of the circulant superposition operation. 2 Let the arrays FE and KE be expressed in vector form as the J × 1 vectors fE and kE, respectively. Then, the circulant superposition operator can be written as k E = CfE (7.3-4) 2 2 where C is a J × J matrix containing elements of the array HE. The circulant superposition operator can then be conveniently expressed in terms of J × J subma- trices Cmn as given by C 1, 1 0 0 … 0 C 1, J – L + 2 … C 1, J C 2, 1 C 2, 2 0 … 0 0 … … · 0 C L – 1, J … C 2, 1 C L, 2 0 … 0 C = (7.3-5) 0 C L + 1, 2 … … … 0 … 0 … 0 C J, J – L + 1 C J, J – L + 2 … C J, J where Cm n 2 ( m 1, n 1 ) = H E ( k 1, k 2 ; m 1, m 2 ) (7.3-6) 2,
CIRCULANT SUPERPOSITION AND CONVOLUTION 179 11 12 13 H = 21 22 23 31 32 33 11 0 31 21 0 0 0 0 13 0 33 23 12 0 32 22 21 11 0 31 0 0 0 0 23 13 0 33 22 12 0 32 31 21 11 0 0 0 0 0 33 23 13 0 32 22 12 0 0 31 21 11 0 0 0 0 0 33 23 13 0 32 22 12 12 0 32 22 11 0 31 21 0 0 0 0 13 0 33 23 22 12 0 32 21 11 0 31 0 0 0 0 23 13 0 33 32 22 12 0 31 21 11 0 0 0 0 0 33 23 13 0 0 32 22 12 0 31 21 11 0 0 0 0 0 33 23 13 C= 13 0 33 23 12 0 32 22 11 0 31 21 0 0 0 0 23 13 0 33 22 12 0 32 21 11 0 31 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 0 0 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 0 0 0 0 0 0 0 0 13 0 33 23 12 0 32 22 11 0 31 21 0 0 0 0 23 13 0 33 22 12 0 32 21 11 0 31 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 0 0 0 0 0 0 33 23 13 0 32 22 12 0 31 21 11 (a) (b) FIGURE 7.3-1. Circulant convolution operators: (a) General impulse array, J = 4, L = 3; (b) Gaussian-shaped impulse array, J = 16, L = 9.
180 SUPERPOSITION AND CONVOLUTION for 1 ≤ n i ≤ J and 1 ≤ m i ≤ J with k i = ( m i – n i + 1 ) modulo J and HE(0, 0) = 0. It should be noted that each row and column of C contains L nonzero submatrices. If the impulse response array is spatially invariant, then Cm = Cm (7.3-7) 2, n2 2 + 1, n 2 + 1 and the submatrices of the rows (columns) can be obtained by a circular shift of the first row (column). Figure 7.3-la illustrates the circulant convolution operator for 16 × 16 (J = 4) data and filtered data arrays and for a 3 × 3 (L = 3) impulse response array. In Figure 7.3-lb, the operator is shown for J = 16 and L = 9 with a Gaussian- shaped impulse response. Finally, when the impulse response is spatially invariant and orthogonally separa- ble, C = CC ⊗ CR (7.3-8) where C R and CC are J × J matrices of the form hR ( 1 ) 0 … 0 h R ( L ) … hR ( 3 ) hR ( 2 ) hR ( 2 ) hR ( 1 ) … 0 0 hR ( 3 ) … … … … hR ( L – 1 ) … 0 hR ( L ) (7.3-9) CR = hR ( L ) hR ( L – 1 ) 0 0 hR ( L ) … 0 … 0 … 0 hR ( L ) … … h (2) h (1) R R Two-dimensional circulant convolution may then be computed as T KE = C C FE C R (7.3-10) 7.4. SUPERPOSITION AND CONVOLUTION OPERATOR RELATIONSHIPS The elements of the finite-area superposition operator D and the elements of the sampled image superposition operator B can be extracted from circulant superposi- tion operator C by use of selection matrices defined as (2)