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

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

0
76
lượt xem
5

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

Mô tả tài liệu

UNITARY TRANSFORMS Two-dimensional unitary transforms have found two major applications in image processing. Transforms have been utilized to extract features from images. For example, with the Fourier transform, the average value or dc term is proportional to the average image amplitude, and the high-frequency terms (ac term) give an indication of the amplitude and orientation of edges within an image.

Chủ đề:

Bình luận(0)

Lưu

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

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) 8 UNITARY TRANSFORMS Two-dimensional unitary transforms have found two major applications in image processing. Transforms have been utilized to extract features from images. For example, with the Fourier transform, the average value or dc term is proportional to the average image amplitude, and the high-frequency terms (ac term) give an indica- tion of the amplitude and orientation of edges within an image. Dimensionality reduction in computation is a second image processing application. Stated simply, those transform coefficients that are small may be excluded from processing opera- tions, such as filtering, without much loss in processing accuracy. Another applica- tion in the field of image coding is transform image coding, in which a bandwidth reduction is achieved by discarding or grossly quantizing low-magnitude transform coefficients. In this chapter we consider the properties of unitary transforms com- monly used in image processing. 8.1. GENERAL UNITARY TRANSFORMS A unitary transform is a specific type of linear transformation in which the basic lin- ear operation of Eq. 5.4-1 is exactly invertible and the operator kernel satisfies cer- tain orthogonality conditions (1,2). The forward unitary transform of the N 1 × N 2 image array F ( n1, n 2 ) results in a N1 × N 2 transformed image array as defined by N1 N2 F ( m 1, m 2 ) = ∑ ∑ F ( n 1, n 2 )A ( n 1, n 2 ; m 1, m 2 ) (8.1-1) n1 = 1 n2 = 1 185
2. 186 UNITARY TRANSFORMS where A ( n1, n 2 ; m1 , m 2 ) represents the forward transform kernel. A reverse or inverse transformation provides a mapping from the transform domain to the image space as given by N1 N2 F ( n 1, n 2 ) = ∑ ∑ F ( m 1, m 2 )B ( n 1, n 2 ; m 1, m2 ) (8.1-2) m1 = 1 m2 = 1 where B ( n1, n 2 ; m 1, m2 ) denotes the inverse transform kernel. The transformation is unitary if the following orthonormality conditions are met: ∑ ∑ A ( n1, n2 ; m1, m2 )A∗ ( j1, j2 ; m1, m2 ) = δ ( n1 – j1, n 2 – j2 ) (8.1-3a) m1 m2 ∑ ∑ B ( n1, n 2 ; m1, m2 )B∗ ( j1, j2 ; m1, m2 ) = δ ( n 1 – j 1, n 2 – j2 ) (8.1-3b) m1 m 2 ∑ ∑ A ( n1, n2 ; m 1, m 2 )A∗ ( n 1, n 2 ; k 1, k 2 ) = δ ( m 1 – k 1, m 2 – k 2 ) (8.1-3c) n1 n2 ∑ ∑ B ( n1, n 2 ; m1, m2 )B∗ ( n1, n2 ; k1, k2 ) = δ ( m 1 – k 1, m 2 – k 2 ) (8.1-3c) n1 n2 The transformation is said to be separable if its kernels can be written in the form A ( n1, n 2 ; m 1, m 2 ) = A C ( n 1, m 1 )AR ( n 2, m 2 ) (8.1-4a) B ( n1, n 2 ; m 1, m 2 ) = B C ( n 1, m 1 )BR ( n 2, m 2 ) (8.1-4b) where the kernel subscripts indicate row and column one-dimensional transform operations. A separable two-dimensional unitary transform can be computed in two steps. First, a one-dimensional transform is taken along each column of the image, yielding N1 P ( m 1, n 2 ) = ∑ F ( n 1, n 2 )A C ( n 1, m 1 ) (8.1-5) n1= 1 Next, a second one-dimensional unitary transform is taken along each row of P ( m1, n 2 ), giving N2 F ( m 1, m 2 ) = ∑ P ( m1, n 2 )A R ( n 2, m 2 ) (8.1-6) n2 = 1
3. GENERAL UNITARY TRANSFORMS 187 Unitary transforms can conveniently be expressed in vector-space form (3). Let F and f denote the matrix and vector representations of an image array, and let F and f be the matrix and vector forms of the transformed image. Then, the two-dimen- sional unitary transform written in vector form is given by f = Af (8.1-7) where A is the forward transformation matrix. The reverse transform is f f = Bf (8.1-8) where B represents the inverse transformation matrix. It is obvious then that –1 B = A (8.1-9) For a unitary transformation, the matrix inverse is given by –1 T A = A∗ (8.1-10) and A is said to be a unitary matrix. A real unitary matrix is called an orthogonal matrix. For such a matrix, –1 T A = A (8.1-11) If the transform kernels are separable such that A = AC ⊗ AR (8.1-12) where A R and A C are row and column unitary transform matrices, then the trans- formed image matrix can be obtained from the image matrix by T F = A C FA R (8.1-13a) The inverse transformation is given by T F = BC F BR (8.1-13b)
4. 188 UNITARY TRANSFORMS –1 –1 where B C = A C and BR = A R . Separable unitary transforms can also be expressed in a hybrid series–vector space form as a sum of vector outer products. Let a C ( n 1 ) and a R ( n 2 ) represent rows n1 and n2 of the unitary matrices AR and AR, respectively. Then, it is easily verified that N1 N2 T F = ∑ ∑ F ( n 1, n 2 )a C ( n 1 )a R ( n 2 ) (8.1-14a) n1 = 1 n2 = 1 Similarly, N1 N2 T F = ∑ ∑ F ( m 1, m 2 )b C ( m 1 )b R ( m 2 ) (8.1-14b) m1 = 1 m2 = 1 where b C ( m 1 ) and b R ( m 2 ) denote rows m1 and m2 of the unitary matrices BC and BR, respectively. The vector outer products of Eq. 8.1-14 form a series of matrices, called basis matrices, that provide matrix decompositions of the image matrix F or its unitary transformation F. There are several ways in which a unitary transformation may be viewed. An image transformation can be interpreted as a decomposition of the image data into a generalized two-dimensional spectrum (4). Each spectral component in the trans- form domain corresponds to the amount of energy of the spectral function within the original image. In this context, the concept of frequency may now be generalized to include transformations by functions other than sine and cosine waveforms. This type of generalized spectral analysis is useful in the investigation of specific decom- positions that are best suited for particular classes of images. Another way to visual- ize an image transformation is to consider the transformation as a multidimensional rotation of coordinates. One of the major properties of a unitary transformation is that measure is preserved. For example, the mean-square difference between two images is equal to the mean-square difference between the unitary transforms of the images. A third approach to the visualization of image transformation is to consider Eq. 8.1-2 as a means of synthesizing an image with a set of two-dimensional mathe- matical functions B ( n1, n 2 ; m 1, m 2 ) for a fixed transform domain coordinate ( m 1, m2 ) . In this interpretation, the kernel B ( n 1, n 2 ; m 1, m 2 ) is called a two-dimen- sional basis function and the transform coefficient F ( m1, m 2 ) is the amplitude of the basis function required in the synthesis of the image. In the remainder of this chapter, to simplify the analysis of two-dimensional uni- tary transforms, all image arrays are considered square of dimension N. Further- more, when expressing transformation operations in series form, as in Eqs. 8.1-1 and 8.1-2, the indices are renumbered and renamed. Thus the input image array is denoted by F(j, k) for j, k = 0, 1, 2,..., N - 1, and the transformed image array is rep- resented by F(u, v) for u, v = 0, 1, 2,..., N - 1. With these definitions, the forward uni- tary transform becomes
5. FOURIER TRANSFORM 189 N–1 N–1 F ( u, v ) = ∑ ∑ F ( j, k )A ( j, k ; u, v ) (8.1-15a) j=0 k=0 and the inverse transform is N–1 N–1 F ( j, k ) = ∑ ∑ F ( u, v )B ( j, k ; u, v ) (8.1-15b) u=0 v=0 8.2. FOURIER TRANSFORM The discrete two-dimensional Fourier transform of an image array is defined in series form as (5–10) N–1 N–1 1  – 2πi  F ( u, v ) = --- N - ∑ ∑ F ( j, k ) exp  ----------- ( uj + vk )   N  (8.2-1a) j=0 k=0 where i = – 1 , and the discrete inverse transform is given by N–1 N–1 1  2πi  F ( j, k ) = --- N - ∑ ∑ F ( u, v ) exp  -------- ( uj + vk )   N  (8.2-1b) u=0 v=0 The indices (u, v) are called the spatial frequencies of the transformation in analogy with the continuous Fourier transform. It should be noted that Eq. 8.2-1 is not uni- versally accepted by all authors; some prefer to place all scaling constants in the inverse transform equation, while still others employ a reversal in the sign of the kernels. Because the transform kernels are separable and symmetric, the two dimensional transforms can be computed as sequential row and column one-dimensional trans- forms. The basis functions of the transform are complex exponentials that may be decomposed into sine and cosine components. The resulting Fourier transform pairs then become  – 2πi   2π   2π  A ( j, k ; u, v ) = exp  ----------- ( uj + vk )  = cos  ----- ( uj + vk )  – i sin  ----- ( uj + vk )  - - (8.2-2a)  N  N  N   2πi   2π   2π  B ( j, k ; u, v ) = exp  ------- ( uj + vk )  = cos  ----- ( uj + vk )  + i sin  ----- ( uj + vk )  - - - (8.2-2b)  N   N   N  Figure 8.2-1 shows plots of the sine and cosine components of the one-dimensional Fourier basis functions for N = 16. It should be observed that the basis functions are a rough approximation to continuous sinusoids only for low frequencies; in fact, the
6. 190 UNITARY TRANSFORMS FIGURE 8.2-1 Fourier transform basis functions, N = 16. highest-frequency basis function is a square wave. Also, there are obvious redun- dancies between the sine and cosine components. The Fourier transform plane possesses many interesting structural properties. The spectral component at the origin of the Fourier domain N–1 N–1 1 F ( 0, 0 ) = --- N - ∑ ∑ F ( j, k ) (8.2-3) j=0 k=0 is equal to N times the spatial average of the image plane. Making the substitutions u = u + mN , v = v + nN in Eq. 8.2-1, where m and n are constants, results in
7. FOURIER TRANSFORM 191 FIGURE 8.2-2. Periodic image and Fourier transform arrays. N–1 N–1 1  – 2πi  F ( u + mN, v + nN ) = --- N - ∑ ∑ F ( j, k ) exp  ----------- ( uj + vk )  exp { –2πi ( mj + nk ) }  N  j=0 k=0 (8.2-4) For all integer values of m and n, the second exponential term of Eq. 8.2-5 assumes a value of unity, and the transform domain is found to be periodic. Thus, as shown in Figure 8.2-2a, F ( u + mN, v + nN ) = F ( u, v ) (8.2-5) for m, n = 0, ± 1, ± 2, … . The two-dimensional Fourier transform of an image is essentially a Fourier series representation of a two-dimensional field. For the Fourier series representation to be valid, the field must be periodic. Thus, as shown in Figure 8.2-2b, the original image must be considered to be periodic horizontally and vertically. The right side of the image therefore abuts the left side, and the top and bottom of the image are adjacent. Spatial frequencies along the coordinate axes of the transform plane arise from these transitions. If the image array represents a luminance field, F ( j, k ) will be a real positive function. However, its Fourier transform will, in general, be complex. Because the 2 transform domain contains 2N components, the real and imaginary, or phase and magnitude components, of each coefficient, it might be thought that the Fourier transformation causes an increase in dimensionality. This, however, is not the case because F ( u, v ) exhibits a property of conjugate symmetry. From Eq. 8.2-4, with m and n set to integer values, conjugation yields
8. 192 UNITARY TRANSFORMS FIGURE 8.2-3. Fourier transform frequency domain. N–1 N–1 1  – 2πi  F * ( u + mN, v + nN ) = --- N - ∑ ∑ F ( j, k ) exp  ----------- ( uj + vk )   N  (8.2-6) j=0 k=0 By the substitution u = – u and v = – v it can be shown that F ( u, v ) = F * ( – u + mN, – v + nN ) (8.2-7) for n = 0, ± 1, ±2, … . As a result of the conjugate symmetry property, almost one- half of the transform domain samples are redundant; that is, they can be generated from other transform samples. Figure 8.2-3 shows the transform plane with a set of redundant components crosshatched. It is possible, of course, to choose the left half- plane samples rather than the upper plane samples as the nonredundant set. Figure 8.2-4 shows a monochrome test image and various versions of its Fourier transform, as computed by Eq. 8.2-1a, where the test image has been scaled over unit range 0.0 ≤ F ( j, k ) ≤ 1.0. Because the dynamic range of transform components is much larger than the exposure range of photographic film, it is necessary to com- press the coefficient values to produce a useful display. Amplitude compression to a unit range display array D ( u, v ) can be obtained by clipping large-magnitude values according to the relation
9. FOURIER TRANSFORM 193 (a) Original (b) Clipped magnitude, nonordered (c) Log magnitude, nonordered (d) Log magnitude, ordered FIGURE 8.2-4. Fourier transform of the smpte_girl_luma image.  1.0 if F ( u, v ) ≥ c F max (8.2-8a)  D ( u, v ) =  F ( u, v )  -------------------- c F max -  if F ( u, v ) < c F max (8.2-8b) where 0.0 < c ≤ 1.0 is the clipping factor and F max is the maximum coefficient magnitude. Another form of amplitude compression is to take the logarithm of each component as given by log { a + b F ( u, v ) } D ( u, v ) = ------------------------------------------------ - (8.2-9) log { a + b F max }
10. 194 UNITARY TRANSFORMS where a and b are scaling constants. Figure 8.2-4b is a clipped magnitude display of the magnitude of the Fourier transform coefficients. Figure 8.2-4c is a logarithmic display for a = 1.0 and b = 100.0. In mathematical operations with continuous signals, the origin of the transform domain is usually at its geometric center. Similarly, the Fraunhofer diffraction pat- tern of a photographic transparency of transmittance F ( x, y ) produced by a coher- ent optical system has its zero-frequency term at the center of its display. A computer-generated two-dimensional discrete Fourier transform with its origin at its center can be produced by a simple reordering of its transform coefficients. Alterna- tively, the quadrants of the Fourier transform, as computed by Eq. 8.2-la, can be j+k reordered automatically by multiplying the image function by the factor ( – 1 ) prior to the Fourier transformation. The proof of this assertion follows from Eq. 8.2-4 with the substitution m = n = 1 . Then, by the identity -- 2 - j+k exp { iπ ( j + k ) } = ( – 1 ) (8.2-10) Eq. 8.2-5 can be expressed as N–1 N–1 1 j+k  – 2πi  F ( u + N ⁄ 2, v + N ⁄ 2 ) = --- N - ∑ ∑ F ( j, k ) ( – 1 ) exp  ----------- ( uj + vk )   N  j=0 k=0 (8.2-11) Figure 8.2-4d contains a log magnitude display of the reordered Fourier compo- nents. The conjugate symmetry in the Fourier domain is readily apparent from the photograph. The Fourier transform written in series form in Eq. 8.2-1 may be redefined in vector-space form as f = Af (8.2-12a) T f = A∗ f (8.2-12b) where f and f are vectors obtained by column scanning the matrices F and F, respectively. The transformation matrix A can be written in direct product form as A = AC ⊗ A R (8.2-13)
11. COSINE, SINE, AND HARTLEY TRANSFORMS 195 where 0 0 0 0 W W W … W 0 1 2 N–1 W W W … W AR = AC = 0 2 4 2(N – 1) (8.2-14) W W W … W … … 2 0 · · (N – 1) W … W with W = exp { – 2πi ⁄ N }. As a result of the direct product decomposition of A, the image matrix and transformed image matrix are related by F = A C FA R (8.2-15a) F = A C∗ F A R∗ (8.2-15b) The properties of the Fourier transform previously proved in series form obviously hold in the matrix formulation. One of the major contributions to the field of image processing was the discovery (5) of an efficient computational algorithm for the discrete Fourier transform (DFT). Brute-force computation of the discrete Fourier transform of a one-dimensional 2 sequence of N values requires on the order of N complex multiply and add opera- tions. A fast Fourier transform (FFT) requires on the order of N log N operations. For large images the computational savings are substantial. The original FFT algo- rithms were limited to images whose dimensions are a power of 2 (e.g., 9 N = 2 = 512 ). Modern algorithms exist for less restrictive image dimensions. Although the Fourier transform possesses many desirable analytic properties, it has a major drawback: Complex, rather than real number computations are necessary. Also, for image coding it does not provide as efficient image energy compaction as other transforms. 8.3. COSINE, SINE, AND HARTLEY TRANSFORMS The cosine, sine, and Hartley transforms are unitary transforms that utilize sinusoidal basis functions, as does the Fourier transform. The cosine and sine transforms are not simply the cosine and sine parts of the Fourier transform. In fact, the cosine and sine parts of the Fourier transform, individually, are not orthogonal functions. The Hartley transform jointly utilizes sine and cosine basis functions, but its coefficients are real numbers, as contrasted with the Fourier transform whose coefficients are, in general, complex numbers.
12. 196 UNITARY TRANSFORMS 8.3.1. Cosine Transform The cosine transform, discovered by Ahmed et al. (12), has found wide application in transform image coding. In fact, it is the foundation of the JPEG standard (13) for still image coding and the MPEG standard for the coding of moving images (14). The forward cosine transform is defined as (12) N–1 N–1 2 π  π  ∑ ∑ F ( j, k ) cos  --- [ u ( j + --- ) ]  cos  --- [ v ( k + --- ) ]  1 1 F ( u, v ) = --- C ( u )C ( v ) - - - N N 2  N 2  j=0 k=0 (8.3-1a) N–1 N–1 2 π  π  ∑ ∑ C ( u )C ( v )F ( u, v ) cos  --- [ u ( j + --- ) ]  cos  --- [ v ( k + --- ) ]  1 1 F ( j, k ) = --- - - - N N 2  N 2  j=0 k=0 (8.3-1b) –1 ⁄ 2 where C ( 0 ) = ( 2 ) and C ( w ) = 1 for w = 1, 2,..., N – 1. It has been observed that the basis functions of the cosine transform are actually a class of discrete Che- byshev polynomials (12). Figure 8.3-1 is a plot of the cosine transform basis functions for N = 16. A photo- graph of the cosine transform of the test image of Figure 8.2-4a is shown in Figure 8.3-2a. The origin is placed in the upper left corner of the picture, consistent with matrix notation. It should be observed that as with the Fourier transform, the image energy tends to concentrate toward the lower spatial frequencies. The cosine transform of a N × N image can be computed by reflecting the image about its edges to obtain a 2N × 2N array, taking the FFT of the array and then extracting the real parts of the Fourier transform (15). Algorithms also exist for the direct computation of each row or column of Eq. 8.3-1 with on the order of N log N real arithmetic operations (12,16). 8.3.2. Sine Transform The sine transform, introduced by Jain (17), as a fast algorithmic substitute for the Karhunen–Loeve transform of a Markov process is defined in one-dimensional form by the basis functions 2  ( j + 1 ) ( u + 1 )π  A ( u, j ) = ------------ sin  -------------------------------------  - - (8.3-2) N+1  N+1  for u, j = 0, 1, 2,..., N – 1. Consider the tridiagonal matrix
13. COSINE, SINE, AND HARTLEY TRANSFORMS 197 FIGURE 8.3-1. Cosine transform basis functions, N = 16. · 1 –α 0 … 0 · · –α 1 –α · · · · T = · · · · (8.3-3) · · · –α 1 –α 0 … 0 –α 1 2 where α = ρ ⁄ ( 1 + ρ ) and 0.0 ≤ ρ ≤ 1.0 is the adjacent element correlation of a Markov process covariance matrix. It can be shown (18) that the basis functions of
14. 198 UNITARY TRANSFORMS (a) Cosine (b) Sine (c) Hartley FIGURE 8.3-2. Cosine, sine, and Hartley transforms of the smpte_girl_luma image, log magnitude displays Eq. 8.3-2, inserted as the elements of a unitary matrix A, diagonalize the matrix T in the sense that T ATA = D (8.3-4) Matrix D is a diagonal matrix composed of the terms 2 1–ρ D ( k, k ) = ------------------------------------------------------------------------ (8.3-5) 2 1 – 2ρ cos { kπ ⁄ ( N + 1 ) } + ρ for k = 1, 2,..., N. Jain (17) has shown that the cosine and sine transforms are interre- lated in that they diagonalize a family of tridiagonal matrices.
15. COSINE, SINE, AND HARTLEY TRANSFORMS 199 FIGURE 8.3-3. Sine transform basis functions, N = 15. The two-dimensional sine transform is defined as N–1 N–1 2  ( j + 1 ) ( u + 1 )π   ( k + 1 ) ( v + 1 )π  F ( u, v ) = ------------ N+1 - ∑ ∑ F ( j, k ) sin  -------------------------------------   N+1 -  sin  -------------------------------------  (8.3-6)  N+1  j=0 k=0 Its inverse is of identical form. Sine transform basis functions are plotted in Figure 8.3-3 for N = 15. Figure 8.3-2b is a photograph of the sine transform of the test image. The sine transform can also be computed directly from Eq. 8.3-10, or efficiently with a Fourier trans- form algorithm (17).
16. 200 UNITARY TRANSFORMS 8.3.3. Hartley Transform Bracewell (19,20) has proposed a discrete real-valued unitary transform, called the Hartley transform, as a substitute for the Fourier transform in many filtering appli- cations. The name derives from the continuous integral version introduced by Hart- ley in 1942 (21). The discrete two-dimensional Hartley transform is defined by the transform pair N–1 N–1 1  2π  F ( u, v ) = --- N - ∑ ∑ F ( j, k ) cas  ------ ( uj + vk )   N  (8.3-7a) j=0 k=0 N–1 N–1 1  2π  F ( j, k ) = --- N - ∑ ∑ F ( u, v ) cas  ----- ( uj + vk )   N -  (8.3-7b) u=0 v=0 where casθ ≡ cos θ + sin θ . The structural similarity between the Fourier and Hartley transforms becomes evident when comparing Eq. 8.3-7 and Eq. 8.2-2. It can be readily shown (17) that the cas θ function is an orthogonal function. Also, the Hartley transform possesses equivalent but not mathematically identical structural properties of the discrete Fourier transform (20). Figure 8.3-2c is a photo- graph of the Hartley transform of the test image. The Hartley transform can be computed efficiently by a FFT-like algorithm (20). The choice between the Fourier and Hartley transforms for a given application is usually based on computational efficiency. In some computing structures, the Hart- ley transform may be more efficiently computed, while in other computing environ- ments, the Fourier transform may be computationally superior. 8.4. HADAMARD, HAAR, AND DAUBECHIES TRANSFORMS The Hadamard, Haar, and Daubechies transforms are related members of a family of nonsinusoidal transforms. 8.4.1. Hadamard Transform The Hadamard transform (22,23) is based on the Hadamard matrix (24), which is a square array of plus and minus 1s whose rows and columns are orthogonal. A nor- malized N × N Hadamard matrix satisfies the relation T HH = I (8.4-1) The smallest orthonormal Hadamard matrix is the 2 × 2 Hadamard matrix given by