THE EFFECTS OF VARYING FREQUENCY RESPONSE IN THE LINEAR STAGE OF A HOMOMORPHIC IMAGE PROCESSOR by Brian R. Bell B.S., University of Colorado, 1982 A thesis submitted to the Faculty of the Gradv,ate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of Electrical Engineering 1988
This thesis for the Master of Science degree by Brian R. Bell has been approved for the Department of Electrical Engineering by oug A. Ross Date
1ll Bell, Brian R. (M.S., Electrical Engineering) The Effects of Varying Frequency Response in the Linear Stage of a Homomorphic Image Processor Thesis directed by Associate Professor Douglas A. Ross Background on two-dimensional signal processing, including separable and non-separable systems, and homomorphic filtering is provided, leading to an explanation of the use of a homomorphic image processor. This image processor consists of non-linear input and output stages surrounding a linear filter as the middle stage. The model chosen for the image determines the non-linear stages, reducing the design of any homomorphic image processor for this image model to the design of the linear 2D filter. Examples of the effects of varying the frequency response of separable and non-separable filters used as the linear stage of the processor are exhibited and discussed. The form and content of this abstract are approved. I recommend its publication. A. Ross
CONTENTS CHAPTER I. INTRODUCTION..................... 1 II. BACKGROUND....................... 3 2D Dig ita lSi g n a Pro c e s sin g . 3 2D Discrete Fourier Transform.. 4 Types of 2D Filters............ 5 Introduction to Homomorphic Filtering.................... 7 Ge n era liz e d Sup e r p 0 sit ion. . . 8 Homomorphic Image Processing... 10 Linear Phase FIR Filters....... 13 IMPLEMENTAT I ON AND COMPARISON OF FILTERS ......... 17 Characteristic Systems Imp eme n tat ion. . . . . . . . 1 7 Input Characteristic System.. 17 Output Characteristic System. 19 Linear Filter Stage Implementation............... 19 Separable Filters............ 19 Symmetric Filters............ 21 Comparison of Computing Requirements............... 23 Ch a r act e r i sit i c s 0 f Raw Ima g e . 2 4
Comp a r i son 0 f F i I t ere d Irna g e s . 25 S ymne t ric F i t e r s Variable Amplitude......... 25 S ynme t ric F i t e r s Variable Breakpoint........ 29 Symmetric vs. Separable e r s, Eq u a Amplification and Attenuation................ 32 Symmetric vs. Separable Filters, Unequal Amplification and Attenuation................ 34 Fix e d Poi n t Imp erne n tat ion. . 3 4 Application of Edge De t e c tor s . . . . . . . . . 38 IV. DISCUSSION AND CONCLUSIONS....... 41 Avenues for Further Study...... 44 BIBLIOGRAPHY. . . . . . . . . . . . . . . 46 APPENDIX A. PARALLEL CCMPUTING FOR 2D SIGNAL PROCESSING.. . . . . 47 v
vi FIGURES Figure 1. Homomorphic system, canonic representation . . . . . 9 2. Homomorphic image processor. . . . . . . . . . .. 11 3. Frequency sampling structure filter. . . . . . . . . .. 15 4. FSS frequency response. . . . . . . . . . . . . .. 16 5. Separable filter frequency response. . . . . . . . . .. 20 6. Symmetric filter frequency response . . . . . . . . .. 22 7. DFT of unfiltered image. . . . . . . . . . . . . .. 24
IMAGES Image Pl. Original photograph................... 18 P2. Unfiltered digital image.............. 26 P 3 S ymne t ric f i I t e r, +3 dB, -3 dB........ 2 7 P 4 S ymne t ric f i I t e r, +6 dB, -6 dB........ 2 7 P 5 S ymne t ric f i I t e r, + 1 2 dB, -1 2 dB...... 2 8 P6. Syrrmetric filter, +18 dB, -18 dB...... 28 P 7 S ymne t r i, c f i I t e r, b rea k poi n tat k =0 . 3 0 P 8 S ymne t ric f i I t e r, b rea k poi n tat k=2... 30 P9. Symnetric filter, breakpoint at k=4... 31 PI0. Symnetric filter, breakpoint at k=9... 31 Pll. Symnet r ic f i Iter, +12 dB, -12 dB...... 33 P12. Separable filter, +12 dB, -12 dB...... 33 P13. Symnetric filter, +18 dB, -6 dB....... 35 P14. Separable filter, +18 dB, -6 dB....... 35 PIS. Separable filter, floating point...... 36 P16. Separable filter, fixed point......... 36 P17. fixed point, Increased gaIn...................... 37 P18. Roberts edge detector on raw Image.... 39 P19. Roberts edge detector on enhanc ed image...................... 39
P20. P21. Sobel edge detector on raw image ..... Sobel edge detector on enhanc ed image ..................... viii 40 40
ACKNOWLEDGEMENTS I wish to thank Wyndham Hannaway for his generous technical assistance, the Computational Math Group at the University of Colorado at Denver for allowing me the use of their data processing facilities, and Dr. Doug Ross for his help and encouragement.
CHAPTER I INTRODUCTION The term refers to any process by which an image is transformed to another image, usually including some change to the image which either makes it more pleasing to look at or enhances or emphasizes some information in the image which is obscure in the original. There are many different image processing techniques, including histogram equalization, colorization, edge detection, blurring, and others. For general information on image processing techniques, see  and . Image processing applications include the enhancement of military reconnaissance and LANDSAT photographs. More recently, the importance of image processing in medical applications has grown tremendously. This thesis is concerned with contrast enhancement and dynamic range compression of images. Contrast enhancement allows the viewing of details in an image that may be obscure or even invisible in the original image. Dynamic range compression causes the difference between the light est and darkest areas of an image to be reduced, allowing more economical storage and transmission of the image over a channel. The method of image enhancement chosen for this work was a This processing technique allows for the simultaneous operations of contrast enhancement and dynamic range reduc-
2 tion. The processor is composed of two non-linear stages and one linear 2-dimensional (2D) filter. will be shown that the model of the image that is chosen fIXes the non-linear stages, reducing the design of homomorphic image processor for this problem to the design of the linear stage. The effects of varying the characteristics of the linear stage will be investigated. The next chapter provides background information on 2D signal processing and homomorphic filtering. is followed by a discussion of the implementation of the image processor and examples of the effects of varying the linear stage. Although the image used for this work was a black and white photograph, the results are easily extended to color image processing.
CHAPTER II BACKGROUND Signal processing generally is con cerned with one-dimensional (ID) signals and systems. An example of a ID signal would be any type of voltage signal over time. this case the single dimensional axis is time. Most ID signals encountered in practice have time as the dimension. image processing, however, the signals of interest are two dimensional (2D). Both of the dimensional axes are spatial, as opposed to temporal, and are orthogonal. Since digital image processing is the subject of this thesis, we will be concerned only with discrete images. A image here is defined as a 2D signal composed of a fmite number of pixels arranged in a rectangular Cartesian grid. Each pixel is assigned a number corresponding to the light intensity at that point in the image. A image is a discrete image for which the values of the pixels are quantized, usually to integer values. Just as a discrete ID signal can be fully represented in the fre quency domain as well as the time domain, a discrete image can be fully described in the frequency domain as well as the spatial domain. the case of images, however, the frequency domain representation must be 2D, similar to the spatial representation.
4 Systems that process ID signals can be considered as ID systems. Similarly, systems that process images must be considered to be 2D systems. will be shown that certain 2D systems can be decomposed into a cascade of ID systems; however, the resulting processing is truly 2D. other words, the processing performed cannot be achieved by a single ID system. discrete Fourier The temporal and frequency domain representations of discrete ID signals are linked through the ID Discrete Fourier Transform (DFT). Similarly, the frequency representation of a discrete image is obtained through the 2D DFT. a discrete 2D image is defined as m, where m and the 2D DFT is defined as E E The result is the 2D frequency representation of the spatial image m, Similarly, the original image can be obtained from the frequency representation by performing the 2D inverse DFT: = -E E (2) Thus the spatial image can be thought of as a weighted sum of 2D complex exponen tials. The calculation of the 2D DFT can be simplified from the equa tions shown above. Note that the equation for can be written as = E E
5 As pointed out by Rabiner and Gold , the term within the brackets is a set of ID DFT's obtained by allowing to vary from 0 to Simplifying this expression further, we can write = m=O where is the term in brackets above. This is again a series of ID DFT's, obtained by allowing to vary from 0 to (4) Examining the above expressions, it is seen that the calculation of the 2D DFT can be performed in the following manner. First, each row of pixels is replaced by its ID DFT. Then, each column of the resulting grid is replaced by its ID DFT. Fast Fourier Transform (FFT) techniques for ID DFT calculations are well known. Since the 2D DFT has been reduced to a series of ID DFT calculations, FFT's can be used to speed the calculations. All of the 2D D FT's used in this thesis were performed using a standard inplace decimation-in-time FFT program, originally written for ID applications. The convolution theorem for ID linear time invariant (L TI) systems can be extended to 2D signals and systems. A 2D LTI system will have an impulse response IT is the input to the system, the output m, will be 00 00 mo=-oo The direct calculation of this convolution can be extremely time-consuming.
6 Therefore, filtering problems are often transformed into the frequency domain before being solved. Analogous to the ID case, the DFT of m, can be obtained through the relation = where and are the 2D DFT's of and ( m, respectively. A 2D system (or signal) can be classified as either or A separable system is one for which the impulse response can be expressed as a product of ID responses: Note that this also implies separability of the DFT of the system, A non-separable system is one for which no such factorization is possible. (7) A system is one in which the impulse response is a function of distance from the origin, Le., = This type of system is not ,. direction dependent" -the impulse response at any point in the 2D plane is dependent only on distance from the origin. is easily seen that a circularly symmetric system is not separable. A separable system can be applied to an incoming 2D signal as a cascade of two ID systems. the 2D convolution is used as the starting point, this can be written as
7 = E mo=-oo no=-oo mol m-mo. IT the m variable denotes the column and the variable denotes the row, this equation states that the filtering operation can fIrst be performed along each column using the ID impulse response ( and then a similar opera tion can be performed for the rows, using the impulse response ( m) A simple case was used for the separable filters performed for this thesis, namely that For this case, a single ID filter network was used, fIrst being applied to the rows, and then applied to the columns of the resulting 2D grid. The separable filters used in this thesis were executed in the time domain using a filter network, discussed below. The circularly symmetric filters, however, were performed in the frequency domain, due to the diffIculties of designing 2D filter networks. This was done simply by multiplying the 2D DFT's of the system frequency response and incoming signal, and then performing the inverse 2D DFT on the result. homomorphic Homomorphic filtering is a non-linear filtering technique developed in the mid-1960's by Oppenheim The technique is based on the consideration of non-linear systems that obey a generalized principle of superposition. The principle of generalized superposition will fIrst be discussed in the next section, followed by a description of the use of homomorphic filtering in image processing.
8 Generalized superposition. The principle of superposition states that, for a linear system transformation any two input signals and and any scalar constant c, and c = c the following discussion, let denote a rule for combining input signals (addition, multiplication, etc.), 4> denote a rule for combining system outputs, : denote rule for combining a scalar with an input signal, and, denote a rule for combining system outputs with a scalar. Then a system transformation Q that obeys the generalized super position principle is one for which (12) and = (13) order to determine whether a system obeys this principle, the input operation and output operation 4> must be defmed. other words, a system may obey this principle with a certain and 4>, but with no other operations than these. The systems of interest here may be shown in as in Figure 1. The first system in the cascade, 8, has the property that (14)
and (15) c The effect of this system is to transform the two signals combined with the rule () to a conventional linear sum of signals, and transform the scalar combination rule : to a conventional scalar multiplication. Figure 1 Canonic representation of a homomorphic system. Adapted from A. V. Oppenheim and R. W. Schafer, (Englewood Cliffs, New Jersey: Prentice-Hall, 1975), p.482. The second system in the cascade is simply a linear system, performing the operation (16) and = (17) c
10 The third system, transforms the addition and scalar multipli cation operations to the 4> and : operations, respectively. and Yl 4> = (18) (19) The system is fIxed by the operations and :, thus it is referred to as the Similarly, is 4>. For a given class of signals combined with non-linear operations 4>, and I, the input system and output system are always the same. Therefore, within this class, once the input and output systems have been determined, the design of the homomorphic system reduces to the design of the linear fIlter. Homomorphic Image Processing. Discrete images will be modelled here as a product of two functions: the component and the component The homomorphic image processor to be used is shown in Figure 2. The input operation is simply multiplication. Therefore we can defme the characteristic system for the input as = . Then the operation of the fIrst stage is = (20) +In[!i(m,n)]
Thus the product of the input signals has been transformed to a sum of signals. Here it is assumed that the signals are both strictly positive and real, thus allowing use of the real logarithm. (In practice, an image is a pattern of light intensity, always being equal to or greater than zero. The difficulties around the pixels that equal zero, corresponding to total darkness, will discussed in another section.) Figure 2 Homomorphic image processor. Adapted from A. V. Oppenheim and R. W. Schafer, (Englewood Cliffs, New Jersey: Prentice-Hall, 1975), p.482. The output of the linear stage [1 is (21) The last stage must transform the sum back into a product, so the [1 function is chosen for The operation of this stage is = exp[ 1 exp[ 1
12 Note that the output of the system is guaranteed to be positive, even though it is likely that the output of the linear stage will sometimes be negative. Thus the output image is in accordance with the definition of an image as a pattern of light intensity that must be always greater than or equal to zero. With the characteristic systems defined, the linear stage must now be considered. Returning to the model of the image, the illuminative com ponent is the pattern of light and dark areas in the image, while the reflective component is the pattern of objects and textures that reflect the light contained in the illuminative component. Generally, the illumination across an image does not vary rapidly-most of the information contained in the illumination is contained in the low frequency regions of its 2D DFT. There are exceptions, of course. Sharp-edged shadows and, for instance, a doorway open to dark room will contain high frequency com ponents. The reflective component, however, contains information about the details and textures of the scene. This information will reside mostly in the higher frequency regions of the 2D DFT of the image. Thus, after the transformation from a product to a linear sum is performed by the first stage of the processor, the reflective and illuminative components of the image may be separately processed in the linear stage, due to the fact that they occupy different frequency regions. There will be some overlap, but hopefully the amount that each component spills over into the other's frequency region is negligible. The objective of the filtering process here is to enhance the contrast of the image, while simultaneously compressing the dynamic range.
13 The contrast enhancement will be achieved through amplifying the reflective component. The dynamic range of an image (difference between the highest and lowest valued pixels) is usually large, and due to the illuminative com ponent. Thus the dynamic range compression will be achieved through attenuating the illuminative component. Linear Phase FIR Filters. is desired to introduce a minimum amount of distortion into the image during the filtering operation. One source of distortion is differing group delays through a filter for different fre quency components of the incoming signal. The group delay may be made constant across the frequency axis by using a finite impulse response (FIR) digital filter with the property that ( = where is the length of the impulse response. This condition implies linear phase, which corresponds to constant group delay across frequency. All FIR filter networks used in this work had odd; therefore, only this case will be con sidered here. To demonstrate the linear phase property, the z-transform of the impulse response is written as (N-l)/2 N-l E E (23) n=O -1)/2 = E E n=O n=O [ = ] 2 Then evaluating for on the unit circle = -J = e 2 -3)/2 [[ 1 II 1 1 } cos .(24)
14 The term within braces is real; therefore the phase is /2, implying a group delay of /2 samples for any frequency. The objective of this thesis is to show the effects of varying the fre quency response of the linear stage of a homomorphic image processor. Therefore, it is desirable to use a linear stage whose response can be easily shaped. The (FSS) digital filter satisfies this requirement . also has the capability of displaying the linear phase discussed above. Figure 3 shows the structure of the FSS filter used for the separable filter operations in this work. is a cascade of a FIR with a bank of IIR (infmite impulse response) structures. All of the FSS filters used in this thesis had an impulse response length of =33. Figure 4 shows a typical frequency response of this filter, showing clearly the linear phase characteris tic. The straight lines in the upper plot are a linear interpolation of the fre quency samples, while the other line is the actual frequency response.
15 M, C, z' -r Z N Z-N -r 2 Z M2 C2 Z' _r2 Z' C'6 Z' 2 Z' YM k r Z C 33 r Figure 3 FSS network structure used for separable filters.
z C CJ z c c 15 10 5 0 -5 -10 -15 -20 4 2 0 : 0 4 5 7 10 15 16.5 o 10 15 16.5 Figure 4 -Typical frequency response of a linear-phase FSS filter. 16
CHAPTER IMPLEMENTATION AND COMPARISON OF FILTERS Characteristic Systems Implementation Input characteristic system. The photograph shown in Figure PI used the subject image. digitized to 8-bit resolution, with a pixel grid of 480 rows by 512 columns. the raw digital image, a value of zero indicated total darkness while 255 indicated the brightest possible pixel. The input characteristic system in the homomorphic image proces sor, the [1 function, is fixed for all the ftlters to be used-therefore it only performed once, providing the "'normalized log image'" input for all ftlters. The actual operation done follows. Single 8-bit pixel values were read from the raw image. IT a value equal to zero, it immedi ately changed to equal one. This is how the difficulty of taking the logarithm of a possibly zero value overcome. Then the logarithm of this integer value computed a floating point number. After all the pixel logarithm values were computed and stored in an array (the "'log imageN), the average of these values (the DC component) computed and then subtracted from each pixel log value. Then the pixel values were scaled such that the largest magnitude pixel would just fit into a 12-bit integer without overflow. These 12-bit values were then stored 16-bit integers, to form the normalized log image. The result that the normalized log image
19 consisted of integer values that closely approximated the floating point values, and had a DC component of zero. (The values were quantized to 12 bits instead of 16 to help avoid overflow in performing fixed-point arithmetic filters.) After the linear filter stage, the inverse operation had to be performed to obtain a digital image that could be restored to visual form. The scaling constant used in the input operation had been saved, as well as the DC component that had bee removed. The integer values in the filtered log image were inverse scaled into floating point numbers, then the DC component was added back in. Finally, the exponent function was applied, and the result was rounded and converted to an 8-bit integer value. a pixel value was too large or small to fit into the 8-bit quantity, it was set to either the maximum or minimum value, accordingly. All of the separable filters were implemented using the FSS digital filter discussed above. Filtering was first applied to the rows of the image, and then to the columns. The frequency response of the resulting 2D filter is a product of two identical ID filters, used in different dimensions. For example, if the frequency response of the FSS structure is as shown in Figure 4, the response of the 2D filter is (25) A graphic representation of the 2D response is shown in Figure 5 (containing 1 full quadrant in the frequency plane, 27!"x27!"). The origin is at the top of
20 the figure, with the x-axis going down to the left, and the y-axis going down to the right. Figure 5 -2D frequency response of separable filter based on ID filter of Figure 4. Note that, in this case, the amplification of the high frequencies in the ID filter is equal to the attenuation of the low frequencies. This results in a gain of 1 along the regions close to the center portion of the axes. The attenuation at the origin is twice the attenuation of the ID filter at the ori gin, and similarly the amplification in the center is twice the amplification of the ID filter at its center-point. The filter network was always applied to the rows in the left-toright direction. was found that image discontinuities from the right edge to the left edge contained high-frequency components that were amplified as part of the reflective component of the image. This, of course, is undesirable. would be nice to apply the filter first in one direction, then in the other, continuing to alternate on every row. The discontinuities would not
21 be expected to occur in adjacent pixels on an edge, and could thus be avoided. This, however, would introduce an additional indexing problem into the programming. Instead, every other row reversed before filtering began, and then reversed again after filtering concluded. The left-to right direction on every row maintained. The same technique applies to the columns. Symmetric filters. Unlike the separable filters, the symmetric filters do not lend themselves easily to execution with a filter network. These filters were executed in the frequency domain. Although an arbitrary frequency response could be chosen for fre quency domain filtering, there were two constraints imposed on these filters. First, the output of the filtering operation constrained to be real. This constraint could be met through making the magnitude response an even function the 2D plane, and the phase response an odd function. order to simplify things, the phase response chosen to be zero across the entire plane. Thus the phase linear, an odd function, and simple to implement. The second constraint that the magnitude response should have no discontinuities in the fIrst derivative; in other words, that it should be a smooth function. is known that sharp cutoffs in the magnitude response can cause "ringing" in the output image . This condition met through deriving the 2D magnitude response from the 1D magnitude response of the FSS structure used in the separable filters. The fIrst quadrant of the 2D response obtained by rotating the 1D response 90 degrees around the origin. The other quadrants were then calculated, based on the condition that the entire response must be an even function. The
22 regions in the center of the plane that were not covered by the rotation were simply set to the gain at center point of the ID response, thus forming a smooth plateau in the center. Figure 6 shows a graphic representation of the magnitude response of a symmetric filter derived from the frequency response shown in Figure 4. Figure 6 Magnitude frequency response of the circularly symmetric filter derived from the ID frequency response shown in Figure 4. Note that in this case, the amplification and attenuation are one-half of those in the separable filter. More importantly, the regions around the axes are now amplified, as opposed to the gain of 1 displayed in the separable filter. The execution of the symmetric filters, as mentioned above, was carried out in the frequency domain. First the 2D DFT of the input image was performed, then multiplied with the magnitude response of the filter. An inverse 2D DFT was then executed, resulting in the filtered image.
23 Comparison requirements. order to make a rough comparison of the computing requirements for the two filter methods, some approximations will be made. First the symmetric filter will be con sidered. The performance of a radix-2 =512 FFT requires approximately 9216 real multiplications and 13,824 real additions. Thus the 2D FFT calcu lations for obtaining the 2D DFT of the input image requires about 9.4 mil lion real multiplications and 14 million real additions. The filtering operation itself requires some 524,000 real multiplications. Finally, the inverse DFT operation will require the same amount of computation as the original DFT. The approximate total computing required is thus 19 million multiplications and 28 million additions. Inspection of the FSS network used in the separable filters leads to the conclusion that for each pixel that enters the network, there are 69 mul tiplications and 67 additions performed. Note that the number of multiplica tions could be reduced to 52; the extra ones are present due to ease of pro gramming for many different filters. Using the actual number of multiplica tions performed, the entire filtering operation requires 34 million multiplica tions and 33 million additions. Even if the number of multiplications were reduced, there would still be 25 million mulip lie at ions required. Thus the symmetric filter has a distinct advantage as far as total arithmetic computing power required. This comparison is admittedly a rough estimate and ignores issues such as the data scrambling required for FFT execution and any effects on performance of a limited amount of system memory. It should also be noted that the replacement of the FSS structure used for the separable filter with an appropriate IIR structure would probably reduce the computing requirements significantly.
24 Characteristics 0/ Raw Image Figure 7 shows a graphic representation of the 2D DFT of the raw (unfiltered) image, with the origin and axes in the same positions as the filter responses shown above. Some interesting features of this spectrum can Figure 7 -2D magnitude DFT of unfiltered image. be discerned. First, there is a preponderance of power in the x and y frequencies, compared to frequencies more in the center of the 2D plane. Inspection of the original photograph reveals why this is so. Many of the sharp boundaries between light and dark regions in the photograph are aligned in the x and y directions. These boundaries involve many different frequencies along the axes. The lines in the bricks also contribute to these frequencies.
25 Secondly, the power contained in the low frequencies is large, and falls off rapidly with increase in frequency. This is in accordance with the assumptions made in modelling the image, where the dynamic range (and thus power) was assumed to be much greater in the illuminative component of the image, occupying the lower frequency range. Third, there is very little information in the center region of the plane. This region corresponds to middle frequencies in an approximately 45 degree direction. Inspection of the original photograph reveals that this is feasible, since there are no distinguishable features, except perhaps for the woman's back, that are oriented in this direction. Figure P2 shows an unfiltered version of the original photograph obtained by reproducing the digital image on a CRT and taking a photo graph of it. Unfortunately, a great deal of information was lost in the digiti zation process. Much of this is due to low resolution and non-linearity of the photographic paper in the dark regions. Better performance could possibly have been achieved by digitizing from the negative of the photo, instead of a print. Comparison of Filtered Images filters, variable amplitude. Figures P3 through P6 show the results of filtering the image with four different circularly sym metric filters. All of these filters were based on a ID frequency response shaped like that shown in Figure 5, except that the lower breakpoint in the frequency response was placed at = o. The gains at the DC level and high frequencies were varied in the four different filters. The filter used for Fig ure P3 had attenuation and amplification of 3 dB, Figure P4 used 6 dB,
Figure P2 Unfiltered image after digitization. Figures P5 and P6 used 12 dB and 18 dB, respectively. The image in Figure P3, filtered plus and minus 3 dB, seems sharper than the original digital image (Figure P2). The chair to the right is much more evident, and the shadows on the driveway are clearly visible, unlike the unfiltered digital image. The reduction in dynamic range com bined with the enhanced contrast, although both are slight, has made this image more pleasing to look at, and has increased the information content evident to the observer. Also, the license plate on the car to the right in the garage has become somewhat more obvious. Figure P4 is the result of using a plus and minus 6 dB circularly symmetric filter. The reduction in dynamic range has now become quite obvious, and the shape of the rear of the car in the garage is starting to be
27 Figure P3 Symmetric filter, +3 dB, -3 dB. Figure P4 Symmetric filter, +6 dB, -6 dB.
28 Figure PS Symmetric filter, +12 dB, -12 dB. Figure P6 Symmetric filter, +18 dB, -18 dB.
29 defmed. This shape is entirely invisible in the unfiltered digital image. The shapes of the shadows on the driveway are now more defmed than they were in the previous image, especially on the left side of the image. Figure PS, the gain and attenuation of the filter has been dou bled again, to plus and minus 12 dB. The reduction in dynamic range has become dramatic, so much so that it is annoying. Note the bright outlines of the panels above the garage and the slight haloes around the woman and the young man to the left. These effects are from ringing of the filter in the spa tial domain. This is a different type of ringing than that discussed above. This is similar to the ringing that is observed when a step input is applied to a ID filter that has high gain in the high frequency regions. The shape of the car has now become readily evident, and objects in the back of the garage are just starting to appear. Note the details in the lower right of the garage that are invisible in the unfiltered digital image. Figure P6 is the result of turning the filter up to plus and minus 18 dB. The dynamic range has now been compressed to the point where the filter exhibits behaviour similar to an edge detector. A good deal of informa tion has been lost, due to this extreme compression. For instance, it is now hard to tell the woman's blouse is lighter than her slacks. However, due to the contrast enhancement, the objects in the back and to the right of the garage are more evident than before, and the rear of the car now stands out. The ringing has now increased to the point where there are multiple haloes around the woman and the young man. filters, variable breakpoint. Figures P7 through PI0 show the results of using a symmetrical filter with gain and attenuation of 12 dB, while varying the lower breakpoint of the frequency response. The
30 Figure P7 Symmetric filter, breakpoint at = o. Figure P8 Breakpoint at 2 (211"/33)
31 Figure pg -Breakpoint at 4 (211"/33) Figure PlO -Breakpoint at 9 (211"/33)
32 breakpoint for Figure P7 is 0, for P8 it is 2, and for Figures P9 and PI0 it is 4 and 9, respectively, where the frequency is = (271"/33) The transition width for these filters is 2 frequency samples, meaning that the upper breakpoint is 3 frequency samples above the lower, and the zero crossing occurs approximately 1.5 samples above the lower breakpoint. These images demonstrate that the separation between the illuminative and reflective components of this image is at a rather low fre quency. the fIrst of these images, there is signifIcant amplifIcation of detail. As the breakpoint is moved higher, however, the amount of the reflective component that falls within the frequency range of amplifIcation becomes rapidly less. By the time the breakpoint is at the frequency sample 4, a great deal of the reflective component is being attenuated rather than amplifIed. When the breakpoint is at sample 9, almost all of the power in the image is being attenuated. A good indicator of this effect is the chair to the right and its different appearance in the four pictures. An interesting observation concerning the ringing can be observed in these images. As the breakpoint is moved toward higher frequencies, the frequency of the ringing increases. This is especially obvious around the line dividing the upper portion of the driveway. This would be expected, because the dominant ringing frequency should be the lowest frequency component of the causative edge that is strongly amplifIed. Therefore, as the amplifIcation band of the filter is moved up, so does the dominant ringing frequency. Symmetric Figure PH shows the result of using a symmetric filter, while Figure P12 shows an image obtained with a separable filter. Both filters used 12 dB attenuation at low frequencies and 12 dB amplifIcation at high
33 Figure P11 Symmetric filter, +12 dB, -12 dB. Figure P12 Separable filter, +12 dB, -12 dB.
34 frequencies. Note that the separable filter, however, presents a gain of 1 for frequencies in the neighborhood of the greater part of the axes, as described above. As would be expected, the symmetric filter displays greater contrast enhancement, due to the large part of the image energy that is con tained along the x and y axis frequencies. The separable filter does not effect these frequencies. Note that the halo around the woman in Figure Pll is present along her entire back and leg down to the driveway. the separable filter result, however, the halo is only present along her back, where the frequencies involved are at approximately 45 degrees. This is again due to the fact that the separable filter does not amplify the high fre quencies in the x direction-thus the absence of overshoot along her leg. The same phenomenon is responsible for the absence of the bright haloes around the panels above the garage in the separable filter result. Symmetrie separable, unequal attenuation and ampli/iea Figures P13 and P14 show the results of again using symmetric and separable filters, however this time the high-frequency amplification is 18 dB while the low frequency attenuation is only 6 dB. This results in a separable filter that does amplify high frequency components along the axes, by 12 dB. Again, the symmetric filter is superior in contrast enhancement, due to the reasons given above. Note that the graininess in the lower right of the image seems to be more random in the output from the symmetric filter. This is because this filter treats all frequencies, regardless of orienta tion, in the same manner, while the separable filter is orientation-dependent. Fixed-point implementation. Figure P15 and P16 show the result of using two separable filters, identical except in that the lower image
35 Figure P13 Symmetric filter, +18 dB, -6 dB. Figure P14 Separable filter, +18 dB, -6 dB.
36 Figure PIS Separable filter, floating point. Figure PI6 Separable filter, fixed point.
37 was produced using fIxed-point arithemetic, while the upper used floating point calculations. The amplification and attenuation for each was 6 dB. The fIxed point arithmetic was performed using 16 bits of precision with 12 bits to right of the binary point. No provision was made for truncating over flows to maximum or minimum values. Even with the low gains used in these fIlters, the fIxed-point arith metic used has noticeably degraded the fIxed-point result. It is likely that this degradation is due to overflows occurring, rather than the limited preci sion used. Figure P17 shows a similar fIxed-point result, only this time the amplifIcation has been set to 18 dB. Overflows have Figure P17 Separable fIlter, fIxed point arithmetic, increased gain. degraded this image so severely as to make it nearly useless. Note that no scaling was performed within the fIltering process to
38 prevent overflows. may be possible to achieve better fIXed-point results through the use of such scaling, especially within the IIR structures con tained in the FSS network. Another improvement that might be made is to truncate overflows to the maximum or minimum 16-bit values. edge Figures P18 and PIg show the results of applying a Roberts edge detector flrst to the raw digital image, and secondly, to the output of the homomorphic system using a symmetric fllter with 12 dB amplification and 6 dB attenuation. This shows clearly how the edges in the image, composed of high-frequency components, have been amplifled. The car to the right, invisible in the edge-detected raw image, shows clearly in the edge-detected enhanced image. The bricks at each edge of the image indicate the same phenomenon. Figures P20 and P21 show the result of applying a 5x5 Sobel edge detection operator to the upper half of the same two images. Many edges become apparent within the garage after the flltering operation. Although some of these edges may be simply noise introduced by the high-gain fllter, careful inspection of the original photograph reveals that some of these edges do indicate items present in the photograph and entirely invisible in the raw or flltered images. For instance, there is a shelf in the back of the garage that runs just under the upper edge of the garage entrance. This shelf is barely visible in the original photograph, and is not evident in any of the fll tered images. The application of the Sobel edge detector to the flltered image, however, indicates the existence of this shelf.
39 Figure P18 -Roberts edge detector, raw image. Figure P19 -Roberts edge detector, enhanced image.
40 Figure P20 Sobel edge detector, raw image. Figure P21 Sobel edge detector, enhanced image.
CHAPTER IV DISCUSSION AND CONCLUSIONS The first conclusion that may be drawn from this work is that, for high-quality image enhancement, better techniques must be used in the conversion from photograph to digital image. The resolution used, 480 by 512 pixels, caused a loss of information that may have been important if this work was intended for something like gleaning more information from military reconnaisance photos. For instance, although the license plate on the car is easily read in the original photograph, it is less clear and sometimes unreadable in the mtered images, due to low resolution. A tradeoff is required, however, because any increase in resolution is accompanied by an increase in the amount of computing necessary to perform the image enhancement. Along these same lines, the digitization from a print of the photo graph instead of from a negative caused a loss of information in the dark regions of the image. This is due to the non-linearities displayed by photo graphic paper in the dark regions. Although an expensive process, the tech nology exists to achieve high-resolution digitization directly from 35 mm negatives. This could have improved the results of the mtered images signi ficantly.
42 The modelling of images as a product of two functions that basi cally occupy different frequency regions is useful in developing the homomorphic image processor. This processor can perform simultaneous dynamic range reduction and contrast enhancement, improving the appear ance of an image and even making visible features which were unseen in the original image. Once the input and output stages of the image processor have been determined, the design of any image processor is reduced to that of determining the linear 2D filter which will be the second stage. Thus a non-linear filtering problem has been converted to a linear filter design prob lem. The usefulness of separable filters for this type of image processing shows itself to be of limited value in this case, relative to circularly sym metric filters. Note that the subject image used here contains a great amount of man-made material. is well-known that images of houses, furni ture, and other man-made items generally contain a great deal of power in the horizontally and vertically oriented frequencies, compared to images of natural objects, such as mountains and trees. This was shown to be true in this case. The difficulty with the separable filters is that they are orientation-dependent in the processing of the 2D frequency plane. Specifi cally, the processing of x and y frequencies in the reflective component of the image is different than the other frequencies within this component. may be that separable filters would be more useful in enhancing images of natural scenes, where the x and y frequencies are not dominant. Clearly, the circularly symmetric filters exhibited superior perfor mance here, compared to the separable filters. Additionally, the symmetric filters, performed in the frequency domain, required less computation than
the execution of separable filters in the spatial domain with a frequency sampling structure filter network. 43 One result of this work that was rather surprising was the determi nation that the separation between the reflective and illuminative com ponents of the image occurred at a rather low frequency. The best subjec tive performance, in terms of simultaneous dynamic range reduction and contrast enhancement, was observed when the lower breakpoint of the fre quency response was placed at DC, and a steep transition between attenuation and amplification at higher frequencies was used. The best results were obtained when the gain=l crossover point was approximately = 2, where the frequency is = This frequency is slightly less than where is the Nyquist frequency. This value was lower than fIrst expected. Figures P7-PIO indicate that nearly all of the power in the reflec tive component is contained in frequencies below = 9, and a great deal of it is below = 4. This indicates that improved results could possibly be achieved with more attention given to the low-frequency characteristics of the 2D filters used. particular, a filter with a very steep transition between attenuation and amplification but retaining the linear phase charac teristic might prove useful. IT this filter were to be derived from the FSS filters used here, the impulse response length would have to be increased. Otherwise, the polynomial interpolation property of the nIter breaks down and discontinuities in the frequency response result. The resolution of digiti zation also affects this aspect of the problem. IT the resolution were increased by a factor of 2 x, the apparent frequency boundary between the reflective and illuminative components would appear to be half of their pre vious values.
44 These facts also are significant in regards to the possibility of alias ing. Although this possibility was entirely disregarded in connection with this work, images containing high frequency patterns will exhibit aliasing in the spatial domain . the image used as the basis for this work, both the reflective and illuminative components are contained in low enough frequen cies that aliasing was not observed. Avenues for further study. The ability to characterize the rela tionship between 2D spectral power content and subjective qualities of images is the fIrst step needed for the design of adaptive image processors. For instance, for a given class of images, say the "ideal" power content in the illuminative and reflective components of the image cd'uld be defIned. Then a system could adjust the linear stage of a homomorphic image processor to maintain ideal images as lighting conditions for a scene changed. It may even prove possible to extend the characterization of images to the point where the reflective and illuminative components can be automatically iden tified and the fIlter adapted This would be useful in remote data gathering systems such as mobile robots on other planets or the moon. may be that the image model used here is valid in concept, but does not go far enough. would be interesting to see if a model of an image as a product of more than two functions could be useful. Use of the homomorphic image described here would then be dependent on the assump tion that there are three or more different frequency bands that should be processed differently. The input and output stages would remain as herein, and the linear stage would have to be designed accordingly.
45 The design of 2D filters can be taken much further than it has been in this work. The circularly symmetric filters used here had real-valued 2D DFT's. may prove possible to use complex-valued filters that achieve higher performance far eliminating the observed ringing at high gains. Dedicated hardware to achieve high performance in the homomorphic image processor is another subject that might be investigated. has been shown that the use of a large number of very simple processors arranged in a highly parallelized architecture can lead to very high speed performance for specialized problems . Although the amount of memory required is almost surely great, the use of a highly parallel architecture using inexpensive microprocessors could prove cost-effective for this problem. The memory would probably be the most expensive item in such a system. The use of fIXed-point arithmetic with proper scaling in the calculations to avoid overflows could help reduce the costs of such a system, while improving per formance.
 BmLIOGRAPHY R. C. Gonzalez and P. Wintz, Addison Wesley Publishing Co., Menlo Park, California, 1987. W. K. Pratt, Wiley, New York, 1978. L. R. Rabiner and B. Gold, Prentice-Hall, Englewood Cliffs, New Jersey, 1975. 46 A. V. Oppenheim, R. W. Schafer, and T. G. Stockham, Jr., "Nonlinear Filtering of Multiplied and Convolved Signals", vol. 56, pp. 1264-1291, Aug. 1968. A. V. Oppenheim and R. W. Schafer, Prentice-Hall, Englewood Cliffs, New Jersey, 1975. C. M. Rader and B. Gold, "Digital Filter Design Techniques in the Frequency Domain", vol. 55, pp. 149-171, Feb. 1967. B. R. Bell and H. L. Johnson, "MAXIM-l: Description of Fine Grain Architecture for Solving Maxwell's Equations", Simu lation Councils, Inc., pp. 67-70, 1988.
APPENDIX A PARALLEL COMPUTING FOR 2D SIGNAL PROCESSING All of the processing for this thesis was performed using programs written in language and run on a Sequent multiprocessor computer. This computer consists of eight independent processors using global memory. The implementation of .. C" language on this machine includes extensions to allow relatively simple use of parallel processing techniques. The program to perform the circularly symmetric filtering opera tions made use of the parallel processing capability of the machine in order to speed up the calculations. Since the bulk of processor time is spent in performing the DFT and inverse DFT, these sections of the program were parallelized. The image array to undergo transformation consisted of 512 rows and 512 columns of pixels. Referring to Chapter II, it was shown that the 2D DFT of the image could be performed by fIrst transforming the rows, and then transforming the columns. The transformation of a row is not depen dent on the results of transforming any other row. The same holds true for columns. Therefore, the transforms within the row and column operations can be done in parallel. The actual parallelization was done as follows, using the row transforms as an example. The master processor, the one on which the pro-
48 gram was initiated, packed the image array and performed various memory allocation tasks. When the time came to do the DFT operations, four new processes were spawned on four other processors. These processes each formed a new pointer that referred to one quarter of the rows of the original image array. For instance, the fIrst new process's pointer referred to the fIrst row of the image array, while the second new process's pointer referred to row 129 of the image. Each new process then successively formed linear arrays of each row assigned to it and performed the FFT on the linear array, placing the result back into the image array. The use of independently formed pointers was necessary due to the use of cache memory in the Sequent computer architecture. this way, the execution of the 2D DFTs and inverse DFTs was speeded up by a factor of nearly 4.