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

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

0
47
lượt xem
6

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

Mô tả tài liệu

LINEAR PROCESSING TECHNIQUES Most discrete image processing computational algorithms are linear in nature; an output image array is produced by a weighted linear combination of elements of an input array. The popularity of linear operations stems from the relative simplicity of spatial linear processing as opposed to spatial nonlinear processing. However, for image processing operations, conventional linear processing is often computationally infeasible without efficient computational algorithms because of the large image arrays....

Chủ đề:

Bình luận(0)

Lưu

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

1. 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) 9 LINEAR PROCESSING TECHNIQUES Most discrete image processing computational algorithms are linear in nature; an output image array is produced by a weighted linear combination of elements of an input array. The popularity of linear operations stems from the relative simplicity of spatial linear processing as opposed to spatial nonlinear processing. However, for image processing operations, conventional linear processing is often computation- ally infeasible without efficient computational algorithms because of the large image arrays. This chapter considers indirect computational techniques that permit more efficient linear processing than by conventional methods. 9.1. TRANSFORM DOMAIN PROCESSING Two-dimensional linear transformations have been defined in Section 5.4 in series form as N1 N2 P ( m 1, m 2 ) = ∑ ∑ F ( n 1, n 2 )T ( n 1, n 2 ; m 1, m 2 ) (9.1-1) n1 = 1 n2 = 1 and defined in vector form as p = Tf (9.1-2) It will now be demonstrated that such linear transformations can often be computed more efficiently by an indirect computational procedure utilizing two-dimensional unitary transforms than by the direct computation indicated by Eq. 9.1-1 or 9.1-2. 213
2. 214 LINEAR PROCESSING TECHNIQUES FIGURE 9.1-1. Direct processing and generalized linear filtering; series formulation. Figure 9.1-1 is a block diagram of the indirect computation technique called gen- eralized linear filtering (1). In the process, the input array F ( n1, n 2 ) undergoes a two-dimensional unitary transformation, resulting in an array of transform coeffi- cients F ( u 1, u 2 ) . Next, a linear combination of these coefficients is taken according to the general relation M1 M2 ˜ F ( w 1, w 2 ) = ∑ ∑ F ( u 1, u 2 )T ( u 1, u 2 ; w 1, w 2 ) (9.1-3) u1 = 1 u2 = 1 where T ( u 1, u 2 ; w 1, w 2 ) represents the linear filtering transformation function. Finally, an inverse unitary transformation is performed to reconstruct the processed array P ( m1, m 2 ) . If this computational procedure is to be more efficient than direct computation by Eq. 9.1-1, it is necessary that fast computational algorithms exist for the unitary transformation, and also the kernel T ( u 1, u 2 ; w 1, w 2 ) must be reasonably sparse; that is, it must contain many zero elements. The generalized linear filtering process can also be defined in terms of vector- space computations as shown in Figure 9.1-2. For notational simplicity, let N1 = N2 = N and M1 = M2 = M. Then the generalized linear filtering process can be described by the equations f = [ A 2 ]f (9.1-4a) N ˜ = Tf f f (9.1-4b) –1 p = [A 2 ] ˜ f (9.1-4c) M
3. TRANSFORM DOMAIN PROCESSING 215 FIGURE 9.1-2. Direct processing and generalized linear filtering; vector formulation. 2 2 2 2 where A 2 is a N × N unitary transform matrix, T is a M × N linear filtering N 2 2 transform operation, and A 2 is a M × M unitary transform matrix. From M Eq. 9.1-4, the input and output vectors are related by –1 p = [A 2 ] T [ A 2 ]f (9.1-5) M N Therefore, equating Eqs. 9.1-2 and 9.1-5 yields the relations between T and T given by –1 T = [A 2 ] T [A 2] (9.1-6a) M N –1 T = [A 2 ]T [ A 2 ] (9.1-6b) M N 2 2 If direct processing is employed, computation by Eq. 9.1-2 requires k P ( M N ) oper- ations, where 0 ≤ k P ≤ 1 is a measure of the sparseness of T. With the generalized linear filtering technique, the number of operations required for a given operator are: 4 Forward transform: N by direct transformation 2 2N log 2 N by fast transformation 2 2 Filter multiplication: kT M N 4 Inverse transform: M by direct transformation 2 2M log 2 M by fast transformation
4. 216 LINEAR PROCESSING TECHNIQUES where 0 ≤ k T ≤ 1 is a measure of the sparseness of T. If k T = 1 and direct unitary transform computation is performed, it is obvious that the generalized linear filter- ing concept is not as efficient as direct computation. However, if fast transform algorithms, similar in structure to the fast Fourier transform, are employed, general- ized linear filtering will be more efficient than direct processing if the sparseness index satisfies the inequality 2 2 k T < k P – ------ log 2 N – ----- log 2 M - - (9.1-7) 2 2 M N In many applications, T will be sufficiently sparse such that the inequality will be satisfied. In fact, unitary transformation tends to decorrelate the elements of T caus- ing T to be sparse. Also, it is often possible to render the filter matrix sparse by setting small-magnitude elements to zero without seriously affecting computational accuracy (1). In subsequent sections, the structure of superposition and convolution operators is analyzed to determine the feasibility of generalized linear filtering in these appli- cations. 9.2. TRANSFORM DOMAIN SUPERPOSITION The superposition operations discussed in Chapter 7 can often be performed more efficiently by transform domain processing rather than by direct processing. Figure 9.2-1a and b illustrate block diagrams of the computational steps involved in direct finite area or sampled image superposition. In Figure 9.2-1d and e, an alternative form of processing is illustrated in which a unitary transformation operation is per- formed on the data vector f before multiplication by a finite area filter matrix D or sampled image filter matrix B. An inverse transform reconstructs the output vector. From Figure 9.2-1, for finite-area superposition, because q = Df (9.2-1a) and –1 q = [A 2 ] D [ A 2 ]f (9.2-1b) M N then clearly the finite-area filter matrix may be expressed as –1 D = [A 2 ]D [ A 2 ] (9.2-2a) M N
5. TRANSFORM DOMAIN SUPERPOSITION 217 FIGURE 9.2-1. Data and transform domain superposition.
6. 218 LINEAR PROCESSING TECHNIQUES Similarly, –1 B = [A 2 ]B [ A 2 ] (9.2-2b) M N If direct finite-area superposition is performed, the required number of 2 2 computational operations is approximately N L , where L is the dimension of the impulse response matrix. In this case, the sparseness index of D is L 2 k D =  ---  N - (9.2-3a) 2 2 Direct sampled image superposition requires on the order of M L operations, and the corresponding sparseness index of B is L 2 k B =  ----  - (9.2-3b) M Figure 9.2-1f is a block diagram of a system for performing circulant superposition by transform domain processing. In this case, the input vector kE is the extended data vector, obtained by embedding the input image array F ( n1, n 2 ) in the left cor- ner of a J × J array of zeros and then column scanning the resultant matrix. Follow- ing the same reasoning as above, it is seen that –1 k E = Cf E = [ A 2 ] C [ A 2 ]f (9.2-4a) J J and hence, –1 C = [ A 2 ]C [ A 2 ] (9.2-4b) J J As noted in Chapter 7, the equivalent output vector for either finite-area or sampled image superposition can be obtained by an element selection operation of kE. For finite-area superposition, (M) (M) q = [ S1 J ⊗ S1 J ]k E (9.2-5a) and for sampled image superposition (M) (M) g = [ S2 J ⊗ S2 J ]k E (9.2-5b)
7. TRANSFORM DOMAIN SUPERPOSITION 219 Also, the matrix form of the output for finite-area superposition is related to the extended image matrix KE by (M) (M) T Q = [ S1 J ]K E [ S1 J ] (9.2-6a) For sampled image superposition, (M) (M) T G = [ S2 J ]K E [ S2 J ] (9.2-6b) The number of computational operations required to obtain kE by transform domain processing is given by the previous analysis for M = N = J. 4 Direct transformation 3J 2 2 Fast transformation: J + 4J log 2 J 2 If C is sparse, many of the J filter multiplication operations can be avoided. From the discussion above, it can be seen that the secret to computationally effi- cient superposition is to select a transformation that possesses a fast computational algorithm that results in a relatively sparse transform domain superposition filter matrix. As an example, consider finite-area convolution performed by Fourier domain processing (2,3). Referring to Figure 9.2-1, let A 2 = AK ⊗ AK (9.2-7) K where 1 - ( x – 1) (y – 1 )   AK = ------- W with W ≡ exp  – 2πi  ----------- K  K  (K) 2 for x, y = 1, 2,..., K. Also, let h E denote the K × 1 vector representation of the extended spatially invariant impulse response array of Eq. 7.3-2 for J = K. The Fou- (K) rier transform of h E is denoted as (K) (K ) hE = [ A 2 ]h E (9.2-8) K 2 2 These transform components are then inserted as the diagonal elements of a K × K matrix ( K) ( K) (K) 2 H = diag [ h E ( 1 ), …, h E ( K ) ] (9.2-9)
8. 220 LINEAR PROCESSING TECHNIQUES Then, it can be shown, after considerable manipulation, that the Fourier transform domain superposition matrices for finite area and sampled image convolution can be written as (4) (M) D = H [ PD ⊗ PD ] (9.2-10) for N = M – L + 1 and (N) B = [ PB ⊗ PB ] H (9.2-11) where N = M + L + 1 and –(u – 1 ) (L – 1 ) 1 1 – WM P D ( u, v ) = -------- --------------------------------------------------------------- - - (9.2-12a) M 1 – W M – ( u – 1 ) – W N –( v – 1 ) –( v – 1 ) ( L – 1 ) 1 1 – WN PB ( u, v ) = ------- --------------------------------------------------------------- - - (9.2-12b) N 1 – W M –( u – 1 ) – W N– ( v – 1 ) Thus the transform domain convolution operators each consist of a scalar weighting (K ) matrix H and an interpolation matrix ( P ⊗ P ) that performs the dimensionality con- 2 2 version between the N - element input vector and the M - element output vector. Generally, the interpolation matrix is relatively sparse, and therefore, transform domain superposition is quite efficient. Now, consider circulant area convolution in the transform domain. Following the previous analysis it is found (4) that the circulant area convolution filter matrix reduces to a scalar operator (J ) C = JH (9.2-13) Thus, as indicated in Eqs. 9.2-10 to 9.2-13, the Fourier domain convolution filter matrices can be expressed in a compact closed form for analysis or operational stor- age. No closed-form expressions have been found for other unitary transforms. Fourier domain convolution is computationally efficient because the convolution operator C is a circulant matrix, and the corresponding filter matrix C is of diagonal form. Actually, as can be seen from Eq. 9.1-6, the Fourier transform basis vectors are eigenvectors of C (5). This result does not hold true for superposition in general, nor for convolution using other unitary transforms. However, in many instances, the filter matrices D, B, and C are relatively sparse, and computational savings can often be achieved by transform domain processing.
9. FAST FOURIER TRANSFORM CONVOLUTION 221 Signal Fourier Hadamard (a) Finite length convolution (b) Sampled data convolution (c) Circulant convolution FIGURE 9.2-2. One-dimensional Fourier and Hadamard domain convolution matrices. Figure 9.2-2 shows the Fourier and Hadamard domain filter matrices for the three forms of convolution for a one-dimensional input vector and a Gaussian-shaped impulse response (6). As expected, the transform domain representations are much more sparse than the data domain representations. Also, the Fourier domain circulant convolution filter is seen to be of diagonal form. Figure 9.2-3 illustrates the structure of the three convolution matrices for two-dimensional convolution (4). 9.3. FAST FOURIER TRANSFORM CONVOLUTION As noted previously, the equivalent output vector for either finite-area or sampled image convolution can be obtained by an element selection operation on the extended output vector kE for circulant convolution or its matrix counterpart KE.
10. 222 LINEAR PROCESSING TECHNIQUES Spatial domain Fourier domain (a) Finite-area convolution (b) Sampled image convolution (c) Circulant convolution FIGURE 9.2-3. Two-dimensional Fourier domain convolution matrices. This result, combined with Eq. 9.2-13, leads to a particularly efficient means of con- volution computation indicated by the following steps: 1. Embed the impulse response matrix in the upper left corner of an all-zero J × J matrix, J ≥ M for finite-area convolution or J ≥ N for sampled infinite-area convolution, and take the two-dimensional Fourier trans- form of the extended impulse response matrix, giving
11. FAST FOURIER TRANSFORM CONVOLUTION 223 HE = AJ HE AJ (9.3-1) 2. Embed the input data array in the upper left corner of an all-zero J × J matrix, and take the two-dimensional Fourier transform of the extended input data matrix to obtain FE = A J FE A J (9.3-2) 3. Perform the scalar multiplication K E ( m, n ) = JH E ( m, n )F E ( m, n ) (9.3-3) where 1 ≤ m, n ≤ J . 4. Take the inverse Fourier transform –1 –1 KE = [ A 2 ] HE [ A 2 ] (9.3-4) J J 5. Extract the desired output matrix (M) (M) T Q = [ S1 J ]K E [ S1 J ] (9.3-5a) or (M) (M) T G = [ S2 J ]K E [ S2 J ] (9.3-5b) It is important that the size of the extended arrays in steps 1 and 2 be chosen large enough to satisfy the inequalities indicated. If the computational steps are performed with J = N, the resulting output array, shown in Figure 9.3-1, will contain erroneous terms in a boundary region of width L – 1 elements, on the top and left-hand side of the output field. This is the wraparound error associated with incorrect use of the Fourier domain convolution method. In addition, for finite area (D-type) convolu- tion, the bottom and right-hand-side strip of output elements will be missing. If the computation is performed with J = M, the output array will be completely filled with the correct terms for D-type convolution. To force J = M for B-type convolution, it is necessary to truncate the bottom and right-hand side of the input array. As a conse- quence, the top and left-hand-side elements of the output array are erroneous.
12. 224 LINEAR PROCESSING TECHNIQUES FIGURE 9.3-1. Wraparound error effects. Figure 9.3-2 illustrates the Fourier transform convolution process with proper zero padding. The example in Figure 9.3-3 shows the effect of no zero padding. In both examples, the image has been filtered using a 11 × 11 uniform impulse response array. The source image of Figure 9.3-3 is 512 × 512 pixels. The source image of Figure 9.3-2 is 502 × 502 pixels. It has been obtained by truncating the bot- tom 10 rows and right 10 columns of the source image of Figure 9.3-3. Figure 9.3-4 shows computer printouts of the upper left corner of the processed images. Figure 9.3-4a is the result of finite-area convolution. The same output is realized in Figure 9.3-4b for proper zero padding. Figure 9.3-4c shows the wraparound error effect for no zero padding. In many signal processing applications, the same impulse response operator is used on different data, and hence step 1 of the computational algorithm need not be repeated. The filter matrix HE may be either stored functionally or indirectly as a computational algorithm. Using a fast Fourier transform algorithm, the forward and 2 inverse transforms require on the order of 2J log 2 J operations each. The scalar 2 2 multiplication requires J operations, in general, for a total of J ( 1 + 4 log 2 J ) opera- tions. For an N × N input array, an M × M output array, and an L × L impulse 2 2 response array, finite-area convolution requires N L operations, and sampled 2 2 image convolution requires M L operations. If the dimension of the impulse response L is sufficiently large with respect to the dimension of the input array N, Fourier domain convolution will be more efficient than direct convolution, perhaps by an order of magnitude or more. Figure 9.3-5 is a plot of L versus N for equality
13. FAST FOURIER TRANSFORM CONVOLUTION 225 (a) HE (b) E (c) FE (d ) E (e) KE (f ) E FIGURE 9.3-2. Fourier transform convolution of the candy_502_luma image with proper zero padding, clipped magnitude displays of Fourier images.
14. 226 LINEAR PROCESSING TECHNIQUES (a ) H E (b ) E (c ) F E (d ) E (e ) k E (f ) E FIGURE 9.3-3. Fourier transform convolution of the candy_512_luma image with improper zero padding, clipped magnitude displays of Fourier images.
15. FAST FOURIER TRANSFORM CONVOLUTION 227 0.001 0.002 0.003 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.013 0.013 0.013 0.013 0.013 0.002 0.005 0.007 0.009 0.011 0.014 0.016 0.018 0.021 0.023 0.025 0.025 0.026 0.026 0.026 0.003 0.007 0.010 0.014 0.017 0.020 0.024 0.027 0.031 0.034 0.038 0.038 0.038 0.039 0.039 0.005 0.009 0.014 0.018 0.023 0.027 0.032 0.036 0.041 0.046 0.050 0.051 0.051 0.051 0.051 0.006 0.011 0.017 0.023 0.028 0.034 0.040 0.045 0.051 0.057 0.063 0.063 0.063 0.064 0.064 0.007 0.014 0.020 0.027 0.034 0.041 0.048 0.054 0.061 0.068 0.075 0.076 0.076 0.076 0.076 0.008 0.016 0.024 0.032 0.040 0.048 0.056 0.064 0.072 0.080 0.088 0.088 0.088 0.088 0.088 0.009 0.018 0.027 0.036 0.045 0.054 0.064 0.073 0.082 0.091 0.100 0.100 0.100 0.100 0.101 0.010 0.020 0.031 0.041 0.051 0.061 0.071 0.081 0.092 0.102 0.112 0.112 0.112 0.113 0.113 0.011 0.023 0.034 0.045 0.056 0.068 0.079 0.090 0.102 0.113 0.124 0.124 0.125 0.125 0.125 0.012 0.025 0.037 0.050 0.062 0.074 0.087 0.099 0.112 0.124 0.136 0.137 0.137 0.137 0.137 0.012 0.025 0.037 0.049 0.062 0.074 0.086 0.099 0.111 0.124 0.136 0.136 0.136 0.136 0.136 0.012 0.025 0.037 0.049 0.062 0.074 0.086 0.099 0.111 0.123 0.135 0.135 0.135 0.135 0.135 0.012 0.025 0.037 0.049 0.061 0.074 0.086 0.098 0.110 0.123 0.135 0.135 0.135 0.135 0.134 0.012 0.025 0.037 0.049 0.061 0.074 0.086 0.098 0.110 0.122 0.134 0.134 0.134 0.134 0.134 (a) Finite-area convolution 0.001 0.002 0.003 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.013 0.013 0.013 0.013 0.013 0.002 0.005 0.007 0.009 0.011 0.014 0.016 0.018 0.021 0.023 0.025 0.025 0.026 0.026 0.026 0.003 0.007 0.010 0.014 0.017 0.020 0.024 0.027 0.031 0.034 0.038 0.038 0.038 0.039 0.039 0.005 0.009 0.014 0.018 0.023 0.027 0.032 0.036 0.041 0.046 0.050 0.051 0.051 0.051 0.051 0.006 0.011 0.017 0.023 0.028 0.034 0.040 0.045 0.051 0.057 0.063 0.063 0.063 0.064 0.064 0.007 0.014 0.020 0.027 0.034 0.041 0.048 0.054 0.061 0.068 0.075 0.076 0.076 0.076 0.076 0.008 0.016 0.024 0.032 0.040 0.048 0.056 0.064 0.072 0.080 0.088 0.088 0.088 0.088 0.088 0.009 0.018 0.027 0.036 0.045 0.054 0.064 0.073 0.082 0.091 0.100 0.100 0.100 0.100 0.101 0.010 0.020 0.031 0.041 0.051 0.061 0.071 0.081 0.092 0.102 0.112 0.112 0.112 0.113 0.113 0.011 0.023 0.034 0.045 0.056 0.068 0.079 0.090 0.102 0.113 0.124 0.124 0.125 0.125 0.125 0.012 0.025 0.037 0.050 0.062 0.074 0.087 0.099 0.112 0.124 0.136 0.137 0.137 0.137 0.137 0.012 0.025 0.037 0.049 0.062 0.074 0.086 0.099 0.111 0.124 0.136 0.136 0.136 0.136 0.136 0.012 0.025 0.037 0.049 0.062 0.074 0.086 0.099 0.111 0.123 0.135 0.135 0.135 0.135 0.135 0.012 0.025 0.037 0.049 0.061 0.074 0.086 0.098 0.110 0.123 0.135 0.135 0.135 0.135 0.134 0.012 0.025 0.037 0.049 0.061 0.074 0.086 0.098 0.110 0.122 0.134 0.134 0.134 0.134 0.134 (b) Fourier transform convolution with proper zero padding 0.771 0.700 0.626 0.552 0.479 0.407 0.334 0.260 0.187 0.113 0.040 0.036 0.034 0.033 0.034 0.721 0.655 0.587 0.519 0.452 0.385 0.319 0.252 0.185 0.118 0.050 0.047 0.044 0.044 0.045 0.673 0.612 0.550 0.488 4.426 0.365 0.304 0.243 0.182 0.122 0.061 0.057 0.055 0.055 0.055 0.624 0.569 0.513 0.456 0.399 0.344 0.288 0.234 0.180 0.125 0.071 0.067 0.065 0.065 0.065 0.578 0.528 0.477 0.426 0.374 0.324 0.274 0.225 0.177 0.129 0.081 0.078 0.076 0.075 0.075 0.532 0.488 0.442 0.396 0.350 0.305 0.260 0.217 0.174 0.133 0.091 0.088 0.086 0.085 0.086 0.486 0.448 0.407 0.367 0.326 0.286 0.246 0.208 0.172 0.136 0.101 0.098 0.096 0.096 0.096 0.438 0.405 0.371 0.336 0.301 0.266 0.232 0.200 0.169 0.139 0.110 0108 0.107 0.106 0.106 0.387 0.361 0.333 0.304 0.275 0.246 0.218 0.191 0.166 0.142 0.119 0.118 0.117 0.116 0.116 0.334 0.313 0.292 0.270 0.247 0.225 0.203 0.182 0.163 0.145 0.128 0.127 0.127 0.127 0.127 0.278 0.264 0.249 0.233 0.218 0.202 0.186 0.172 0.159 0.148 0.136 0.137 0.137 0.137 0.137 0.273 0.260 0.246 0.231 0.216 0.200 0.185 0.171 0.158 0.147 0.136 0.136 0.136 0.136 0.136 0.266 0.254 0.241 0.228 0.213 0.198 0.183 0.169 0.157 0.146 0.135 0.135 0.135 0.135 0.135 0.257 0.246 0.234 0.222 0.209 0.195 0.181 0.168 0.156 0.145 0.135 0.135 0.135 0.135 0.134 0.247 0.237 0.227 0.215 0.204 0.192 0.179 0.166 0.155 0.144 0.134 0.134 0.134 0.134 0.134 (c) Fourier transform convolution without zero padding FIGURE 9.3-4. Wraparound error for Fourier transform convolution, upper left corner of processed image. between direct and Fourier domain finite area convolution. The jaggedness of the plot, in this example, arises from discrete changes in J (64, 128, 256,...) as N increases. Fourier domain processing is more computationally efficient than direct process- ing for image convolution if the impulse response is sufficiently large. However, if the image to be processed is large, the relative computational advantage of Fourier domain processing diminishes. Also, there are attendant problems of computational
16. 228 LINEAR PROCESSING TECHNIQUES FIGURE 9.3-5. Comparison of direct and Fourier domain processing for finite-area convo- lution. accuracy with large Fourier transforms. Both difficulties can be alleviated by a block-mode filtering technique in which a large image is separately processed in adjacent overlapped blocks (2, 7–9). Figure 9.3-6a illustrates the extraction of a NB × NB pixel block from the upper left corner of a large image array. After convolution with a L × L impulse response, the resulting M B × M B pixel block is placed in the upper left corner of an output FIGURE 9.3-6. Geometric arrangement of blocks for block-mode filtering.
17. FOURIER TRANSFORM FILTERING 229 data array as indicated in Figure 9.3-6a. Next, a second block of N B × N B pixels is extracted from the input array to produce a second block of M B × M B output pixels that will lie adjacent to the first block. As indicated in Figure 9.3-6b, this second input block must be overlapped by (L – 1) pixels in order to generate an adjacent output block. The computational process then proceeds until all input blocks are filled along the first row. If a partial input block remains along the row, zero-value elements can be added to complete the block. Next, an input block, overlapped by (L –1) pixels with the first row blocks, is extracted to produce the first block of the second output row. The algorithm continues in this fashion until all output points are computed. A total of 2 2 O F = N + 2N log 2 N (9.3-6) operations is required for Fourier domain convolution over the full size image array. With block-mode filtering with N B × N B input pixel blocks, the required number of operations is 2 2 2 O B = R ( N B + 2NB log 2 N ) (9.3-7) where R represents the largest integer value of the ratio N ⁄ ( N B + L – 1 ). Hunt (9) has determined the optimum block size as a function of the original image size and impulse response size. 9.4. FOURIER TRANSFORM FILTERING The discrete Fourier transform convolution processing algorithm of Section 9.3 is often utilized for computer simulation of continuous Fourier domain filtering. In this section we consider discrete Fourier transform filter design techniques. 9.4.1. Transfer Function Generation The first step in the discrete Fourier transform filtering process is generation of the discrete domain transfer function. For simplicity, the following discussion is limited to one-dimensional signals. The extension to two dimensions is straightforward. Consider a one-dimensional continuous signal f C ( x ) of wide extent which is bandlimited such that its Fourier transform f C ( ω ) is zero for ω greater than a cut- off frequency ω 0. This signal is to be convolved with a continuous impulse function h C ( x ) whose transfer function h C ( ω ) is also bandlimited to ω 0 . From Chapter 1 it is known that the convolution can be performed either in the spatial domain by the operation
18. 230 LINEAR PROCESSING TECHNIQUES ∞ gC ( x ) = ∫–∞ fC ( α )hC ( x – α ) dα (9.4-1a) or in the continuous Fourier domain by 1 ∞ g C ( x ) = ----- 2π - ∫–∞ fC ( ω )hC ( ω ) exp { iωx } dω (9.4-1b) Chapter 7 has presented techniques for the discretization of the convolution inte- gral of Eq. 9.4-1. In this process, the continuous impulse response function h C ( x ) must be truncated by spatial multiplication of a window function y(x) to produce the windowed impulse response b C ( x ) = h C ( x )y ( x ) (9.4-2) where y(x) = 0 for x > T . The window function is designed to smooth the truncation effect. The resulting convolution integral is then approximated as x+T gC ( x ) = ∫x – T fC ( α )b C ( x – α ) dα (9.4-3) Next, the output signal g C ( x ) is sampled over 2J + 1 points at a resolution ∆ = π ⁄ ω 0, and the continuous integration is replaced by a quadrature summation at the same resolution ∆ , yielding the discrete representation j+K g C ( j∆ ) = ∑ f C ( k∆ )b C [ ( j – k )∆ ] (9.4-4) k=j–K where K is the nearest integer value of the ratio T ⁄ ∆. Computation of Eq. 9.4-4 by discrete Fourier transform processing requires formation of the discrete domain transfer function b D ( u ) . If the continuous domain impulse response function h C ( x ) is known analytically, the samples of the windowed impulse response function are inserted as the first L = 2K + 1 elements of a J-element sequence and the remaining J – L elements are set to zero. Thus, let b D ( p ) = b C ( – K ), …, b C ( 0 ), …, b C ( K ) , 0, …, 0 (9.4-5)              L terms where 0 ≤ p ≤ P – 1. The terms of b D ( p ) can be extracted from the continuous impulse response function h C ( x ) and the window function by the sampling operation
19. FOURIER TRANSFORM FILTERING 231 b D ( p ) = y ( x )h C ( x )δ ( x – p∆ ) (9.4-6) The next step in the discrete Fourier transform convolution algorithm is to perform a discrete Fourier transform of b D ( p ) over P points to obtain P–1 1  – 2πipu  b D ( u ) = ------- P ∑ b D ( p ) exp  -----------------   P -  (9.4-7) p=1 where 0 ≤ u ≤ P – 1 . If the continuous domain transfer function hC ( ω ) is known analytically, then b D ( u ) can be obtained directly. It can be shown that   b D ( u ) = ---------------- exp  – iπ ( L – 1 )  h C  2πu  1 - ------------------------- - --------- - (9.4-8a) 2 P  P∆  4 Pπ   bD( P – u ) = b * ( u ) D (9.4-8b) for u = 0, 1,..., P/2, where bC ( ω ) = h C ( ω ) * y ( ω ) (9.4-8c) and y ( ω ) is the continuous domain Fourier transform of the window function y(x). If h C ( ω ) and y ( ω ) are known analytically, then, in principle, h C ( ω ) can be obtained by analytically performing the convolution operation of Eq. 9.4-8c and evaluating the resulting continuous function at points 2πu ⁄ P∆. In practice, the analytic convo- lution is often difficult to perform, especially in two dimensions. An alternative is to perform an analytic inverse Fourier transformation of the transfer function h C ( ω ) to obtain its continuous domain impulse response h C ( x ) and then form b D ( u ) from the steps of Eqs. 9.4-5 to 9.4-7. Still another alternative is to form b D ( u ) from h C ( ω ) according to Eqs. 9.4-8a and 9.4-8b, take its discrete inverse Fourier transform, win- dow the resulting sequence, and then form b D ( u ) from Eq. 9.4-7. 9.4.2. Windowing Functions The windowing operation performed explicitly in the spatial domain according to Eq. 9.4-6 or implicitly in the Fourier domain by Eq. 9.4-8 is absolutely imperative if the wraparound error effect described in Section 9.3 is to be avoided. A common mistake in image filtering is to set the values of the discrete impulse response func- tion arbitrarily equal to samples of the continuous impulse response function. The corresponding extended discrete impulse response function will generally possess nonzero elements in each of its J elements. That is, the length L of the discrete
20. 232 LINEAR PROCESSING TECHNIQUES impulse response embedded in the extended vector of Eq. 9.4-5 will implicitly be set equal to J. Therefore, all elements of the output filtering operation will be subject to wraparound error. A variety of window functions have been proposed for discrete linear filtering (10–12). Several of the most common are listed in Table 9.4-1 and sketched in Figure 9.4-1. Figure 9.4-2 shows plots of the transfer functions of these window functions. The window transfer functions consist of a main lobe and sidelobes whose peaks decrease in magnitude with increasing frequency. Examination of the structure of Eq. 9.4-8 indicates that the main lobe causes a loss in frequency response over the signal passband from 0 to ω 0 , while the sidelobes are responsible for an aliasing error because the windowed impulse response function b C ( ω ) is not bandlimited. A tapered window function reduces the magnitude of the sidelobes and consequently attenuates the aliasing error, but the main lobe becomes wider, causing the signal frequency response within the passband to be reduced. A design trade-off must be made between these complementary sources of error. Both sources of degradation can be reduced by increasing the truncation length of the windowed impulse response, but this strategy will either result in a shorter length output sequence or an increased number of computational operations. TABLE 9.4-1. Window Functions a Function Definition Rectangular w(n) = 1 0≤n≤L–1 2n -  ----------- 0≤n–L–1 ----------- - L–1 2 Barlett (triangular) w(n) =   2 – ----------- 2 - L–1 ----------- ≤ n ≤ L – 1 -  L–1 2 1  2πn  Hanning w(n) = --  1 – cos  -----------  - - 0≤n≤L–1 2  L – 1  2πn-   Hamming w(n) = 0.54 - 0.46 cos  -----------  0≤n≤L–1  L – 1  2πn   4πn  Blackman w(n) = 0.42 – 0.5 cos  -----------  + 0.08 cos  -----------  - - 0≤n≤L–1  L–1 L – 1   2 2 1 ⁄ 2 I0  ωa [ ( ( L – 1 ) ⁄ 2 ) – [ n – ( ( L – 1 ) ⁄ 2 ) ] ]  Kaiser   ----------------------------------------------------------------------------------------------------------------- 0 ≤ n ≤ L – 1 - I0 { ωa [ ( L – 1 ) ⁄ 2 ] } a I 0 { · } is the modified zeroth-order Bessel function of the first kind and ω a is a design parameter.