Introduction to Two-Dimensional/Image Processing

MASSACHUSETTS INSTITUTE OF TECHNOLOGY DEPARTMENT OF MECHANICAL ENGINEERING 2.161 Signal Processing – Continuous and Discrete Introduction to Two-Dimensional/Image Processing 1 1 Introduction In this set of notes we extend the concepts of one-dimensional signal processing to two dimensions. Typical fields that involve the topics introduced here are imaging, including optics, and image processing. The most important concept to grasp is that we are no longer dealing with functions of time, and the primary domain is the spatial domain. Therefore, instead of dealing with f (t) we deal with f (x), where x is a spatial distance, or in the two-dimensional situation f (x, y). The methods developed here use Fourier techniques to transform the spatial representa­ tion f (x, y) to a frequency domain space F (j u, j v), where u, and v are spatial frequencies with units of radians/unit distance. As in temporal processing the difficult operation of convolution is replaced by simple multiplication in the frequency domain. In addition, an understanding of Fourier methods gives qualitative insights to image processing techniques such as “sharpening” and “blurring”. In these notes we begin by thoroughly examining one-dimensional Fourier methods and properties to gain insights on the Fourier transform and frequency domain processing, and then generalize to two-dimensional imaging. 2 Two-dimensional Imaging We now address the more usual situation of processing a two-dimensional object function f (x, y) through a linear image processing system, producing a two-dimensional image g(x, y). The sifting property of the Dirac delta function may be extended to two dimensions: g (x ,y ) = f(x ,y ) Ä f(x ,y ) h (x ,y ) L in e a r Im a g e P ro c e s s o r h (x ,y ) In p u t Im a g e O u tp u t Im a g e Figure 1: A two-dimensional linear image processing system. � � ∞ −∞ 1 f (x, y)δ(x − a)δ(y − b)dxdy = f (a, b) D. Rowell, November 20, 2008 1 (1) so that we can write the input f (x, y) as a superposition of two-dimensional impulses δ 2 (x, y) = δ(x)δ(y): � � ∞ O(σ, ρ)δ(x − σ)δ(y − ρ)dσdρ (2) f (x, y) = −∞ Assume that the action of the image processor can be expressed as a two-dimensional linear system operator L � 2 {} that maps the input image f (x, y) to the output image g(x, y), that is (3) g(x, y) = L � 2 {f (x, y)} so that we can write the image as �� � g(x, y) = L �2 � ∞ −∞ f (σ, ρ)δ(x − σ)δ(y − ρ)dσdρ (4) As in one-dimensional signal processing, we can reverse the order of the operator and inte­ gration � � ∞ g(x, y) = f (σ, ρ)�L2 {δ(x − σ)δ(y − ρ)} dσdρ (5) −∞ For a spatially invariant system, we define the two-dimensional point-spread function h(x, y) as the response to a two-dimensional delta function δ 2 (x, y) = δ(x)δ(y) at the origin of the object space: h(x, y) = L � 2 {δ(x)δ(y)} (6) Equation (5) can be written � � ∞ g(x, y) = −∞ f (σ, ρ)h(x − σ, y − ρ)dσdρ (7) which is a two-dimensional convolution between the point-spread function h(x, y) and the object space f (x, y). As in the one-dimensional case, we conclude that the point-spread function completely characterizes the imaging system, in the sense that the imaging fidelity is completely deter­ mined by h(x, y). Example 1 Determine the image formed by convolving a square uniform white image of size 20×20 by a linear image processing filter with point-spread function that is cir­ cularly symmetric with a radius of 5 units and height of unity. Solution: The resulting images shown in Fig. 2 were computed using the fol­ lowing Matlab script: 2 % Matlab script to demonstrate 2-D convolution % First create an array containing the point-spread function psf = zeros(40,40); for i=1:40 for j=1:40 r = sqrt((i-20)^2 + (j-20)^2); if (r < 5) psf(i,j)= 1; end; end; end; % % Now create the object object=zeros(40,40); for i=10:29 for j=10:29 object(i,j)=1; end; end; % % Set up a simple colormap so that all plots are in black and white bw=[0 0 0]; colormap(bw); % % Make 3-D plots of the two input functions mesh(psf); title(’Point-spread function’); pause; mesh(object); title(’Rectangular object’); pause; % % Use the conv2 function to do the 2-D convolution image=conv2(psf,object,’same’); mesh(image); title(’Resulting image’); pause; 3 The Two-Dimensional Fourier Transform In imaging processing work, representing two-dimensional functions, such as images, in the Fourier domain is of great importance. Let f (x, y) be a function of two independent variables x and y; the Fourier transform of f (x, y) is � ∞� ∞ f (x, y)e−j(ux+vy) dxdy (8) F (u, v) = −∞ −∞ 3 1 0.8 0.6 0.4 0.2 0 40 40 30 30 20 20 10 10 0 0 1 0.8 0.6 0.4 0.2 0 40 40 30 30 20 20 10 10 0 0 70 60 50 40 30 20 10 0 40 40 30 30 20 20 10 10 0 0 Figure 2: Two-dimensional convolution: (a) A circular point-spread function, (r = 5), (b) a uniform rectangular 20 ×20 object, and (c) the image formed by convolution. 4 where u and v are spatial frequencies with units of radians per unit length. By splitting the exponential into two halves it can be seen that the transform can be expressed as two one-dimensional transforms; first with respect to x and then with respect to y, � � ∞ �� ∞ −jux F (u, v) = f (x, y)e dx e−jvy dy. (9) −∞ −∞ In general F (u, v) is a complex-valued function of u and v. Example 2 Find the two-dimensional Fourier transform of the two-dimensional impulse f (x, y) = δ 2 (x, y) = δ(x)δ(y). Solution: Separating the variables as in Eq. 9, � � ∞ �� ∞ −jux δ(x)e dx δ(y)e−jvy dy Δ(u, v) = −∞ −∞ and by the sifting property of the delta function, Δ(u, v) = 1 Example 3 Find the Fourier transform of the function f (x, y) = rect(x, y). Solution: The function has a value of unity within the square defined by −0.5 > x < 0.5, y > −0.5 < 0.5. The transform is therefore � 1/2 � 1/2 F (u, v) = −1/2 = e−j(ux+vy) dxdy −1/2 sin u sin v = sinc(u, v) u v as shown in Fig. 3. The two-dimensional inverse Fourier transform is � ∞� ∞ 1 F (u, v)ej(ux+vy) dudv f (x, y) = 2 (2π) −∞ −∞ 5 (10) (a) (b) Figure 3: The two-dimensional rect function and its Fourier Transform. 6 3.1 Properties of the Two-dimensional Fourier Transform The following is list of some of the important properties of the transform. In all cases it is assumed F f (x, y) ⇐⇒ F (u, v). (1) Linearity: F af (x, y) + bg(x, y) ⇐⇒ aF (u, v) + bG(u, v) which follows directly from the linearity of the integration operation. (2) Shift Property: F f (x − α, y − β) ⇐⇒ F (u, v)e−j(uα+vβ) which can be shown by a change in the variables of integration. (3) Similarity (scaling): 1 F f (αx, βy) ⇐⇒ |αβ| F � u v , α β � . (4) Rotation by an Angle: we can write the transform in polar coordinates as F f (r, θ) ⇐⇒ F (ω, φ) If the spatial domain function, f , is rotated by an angle α then the frequency domain function, F , is also rotated by the same angle: F f (r, θ + α) ⇐⇒ F (ω, φ + α) (5) Convolution: The following result is analogous to the one dimensional case: � ∞� ∞ F f (σ, ρ)g(x − σ, y − ρ)dσdρ ⇐⇒ F (u, v)G(u, v) −∞ −∞ so that convolution in the spatial domain is equivalent to multiplication in the frequency domain, or F f (x, y) ⊗ g(x, y) ⇐⇒ F (u, v)G(u, v) The dual holds: 1 f (x, y)g(x, y) ⇐⇒ (2π)2 F � −∞ (6) Parseval’s Theorem: � ∞� ∞ f (x, y)g(x, y)dxdy = −∞ � ∞ −∞ ∞ −∞ 1 (2π)2 F (σ, ρ)G(u − σ, v − ρ)dσdρ � ∞ � ∞ F (u, v)G(u, v)dudv −∞ −∞ where f denotes the complex conjugate. When f (x, y) = g(x, y), we have Rayleigh’s theorem: � ∞� ∞ � ∞� ∞ 1 2 |f (x, y)| dxdy = |F (u, v)|2 dudv 2 (2π) −∞ −∞ −∞ −∞ which is a statement of the conservation of “energy” between the two domains. 7 3.2 Frequency domain description of system response By analogy with the one-dimensional case, if a two-dimensional imaging system has a pointspread function h(x, y) and its response to an object function f (x, y) is g(x, y), or g(x, y) = f (x, y) ⊗ h(x, y), then in the frequency domain G(u, v) = F (u, v)H(u, v) (11) where H(u, v) = F {h(x, y)}, and is the two-dimensional transfer function. As in the onedimensional case, the transfer function completely characterizes a linear imaging system. c o n v o lu tio n in p u t im a g e Im a g e F ilte r f(x ,y ) 2 -D F T h (x ,y ) F ( ju , jv ) 2 -D F T g (x ,y ) H ( ju , jv ) 2 -D o u tp u t im a g e IF T G ( ju , jv ) = F ( ju , jv ) H ( ju , jv ) m u ltip lic a tio n Figure 4: Two dimensional image filtering. 3.3 Frequency Domain Operations as a Filtering Process The multiplication of a spectral representation F (u, v) with a transfer function H(u, v) may be considered as a frequency domain filtering operation in which regions of the input spec­ trum are emphasized or attenuated by the characteristics of the transfer function: Low-pass filter: The low frequency region of the spectrum (adjacent to the origin) is em­ phasized, and the high frequency components are attenuated. This is generally a “blurring” operation that removes the fine detail in the image (and reduces noise). High-pass filter The high frequency region of the spectrum (distant from the origin) is emphasized and the low frequency region attenuated. This results in a sharpening of the image in which discontinuities (edges) are enhanced, but noise in the image is also enhanced. Band-Pass and Band-Reject Filters In these filters the transfer function is chosen to accentuate or attenuate an annular region of the input spectrum. 8 Example 4 Write a Matlab script to demonstrate the operations involved in two-dimensional image generation. Use the Fourier methods to convolve a circularly symmetric cylindrical (r = 5) point-spread function with a rectangular object. Solution: The script is shown below. % Matlab script to demonstrate 2-D convolution % Set up a simple colormap so that all plots are in black and white bw=[0 0 0]; colormap(bw); % % First create an array containing the circular point-spread % function and compute its spectrum % psf = zeros(64,64); for i=1:64 for j=1:64 r = sqrt((i-32)^2 + (j-32)^2); if (r < 5) psf(i,j)= 1; end; end; end; fftpsf=fftshift(fft2(psf)); % % Make 3-D plots of the PSF and its spectrum % figure(1); mesh(psf); title(’Point-spread function’); figure(2); mesh(abs(fftpsf)); title(’Magnitude of FFT of PSF: Modulation Transfer Function’); pause; % % Now create the object % object=zeros(64,64); for i=22:41 for j=22:41 object(i,j)=1; end; end; fftobj=fftshift(fft2(object)); 9 % % Plot these results % figure(1); mesh(object); title(’Rectangular object’); figure(2); mesh(abs(fftobj)); title(’Magnitude of FFT of object’); pause; % % Use the ".*" operator to do the element by element % frequency domain multiplication % fftimg = fftpsf.*fftobj; image=fftshift(ifft2(fftimg)); % % Plot the results % mesh(abs(fftimg)); title(’Product of FFTs of PSF and object’); figure(1); mesh(abs(image)); title(’Resulting image’); The three-dimensional mesh plots produced by this script are shown in Figs. 5 and 6. In Fig. 5 the spatial domain functions of the psf, object and resulting image are shown, in Fig. 6 the equivalent spectral plots are shown. 10 Point−spread function Rectangular object 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 80 0 80 70 60 70 60 60 60 50 40 50 40 40 40 30 20 30 20 20 0 20 10 0 0 10 0 Resulting image 70 60 50 40 30 20 10 0 80 70 60 60 50 40 40 30 20 20 0 10 0 Figure 5: Frequency domain based convolution:(a) the point-spread function, (b) the object, and (c) the resulting image. 11 Magnitude of FFT of PSF: Modulation Transfer Function Magnitude of FFT of object 70 400 60 350 50 300 250 40 200 30 150 20 100 10 50 0 80 0 80 70 60 70 60 60 60 50 40 50 40 40 40 30 20 30 20 20 0 20 10 0 0 10 0 Product of FFTs of PSF and object 4 x 10 3 2.5 2 1.5 1 0.5 0 80 70 60 60 50 40 40 30 20 20 0 10 0 Figure 6: Frequency domain based convolution: the magnitude of the Fourier transforms of (a) the point-spread function, (b) the object, and (c) the resulting image. 12