- Email: [email protected]
Image enhancement by intensity-dependent spread functions
CVGIP:
GRAPHICAL
MODELS
Vol. 54, No. 1, January,
AND
IMAGE
PROCESSING
pp. 45-55, 1992
Image Enhancement by Intensity-Dependent
Spread Functions
GREG REESE Applied Math and Physics Section, Sverdrup Technology,
Incorporated,
Stennis Space Center, Mississippi
39529
Received August 27, 1990; accepted September 10, 1991
If the illumination on a reflectance edge is constant, that illumination cancels out of the ratio of the two intensities, and the peak height of the IDS response to that edge will always be the same. This property implies that IDS can aid in detecting a physical characteristic of objects in a scene-their reflectance. Thus if the same object is present several times in one image under very different constant illuminations, for example, once in sunlight and once in shadow, IDS transforms both objects so that their outputs look quite similar to each other. As some later pictures show, the outputs are not exactly the same be-
This paper describes a method of enhancing images based on a new model of human vision. The model is an adaptive, nonlinear system that enhances edges, nonlinearly compresses dynamic ranges, automatically adjusts itself to local variations in intensity, can (but need not) be carried out on a parallel processor, and has only one adjustable parameter that requires no intricate tuning. It is particularly useful in images that have regions of low contrast, such as brightly lit areas and those covered by shadows. The text explains how to convert the original mathematical theory to one useful for digital image processing and illustrates the output of the system to a variety of images. The focus is on qualitative enhancement of images, but the system would also be useful in machine vision. 0 1992 Academic Press, Inc.
1.
cause images
of objects
get more smeared
as the intensi-
ties get lower. Nevertheless, the result can still be a substantial improvement over the unenhanced image. The IDS output at a point depends on the input intensities near that point and the system adapts itself to changes in these intensities. There have been numerous studies of adaptive systems (for a review see [8]). Although the performance of adaptive systems is often substantially better than that of their nonadaptive counterparts, they do suffer from a number of drawbacks, in general, being complex, computationally expensive, and subject to the vagaries of window size and shape. Since IDS has no windows it does not suffer from the latter problems and although it requires many calculations, the computations are very simple. The next section of the paper is a review of the mathematical formulation of IDS as a system that operates on inputs defined everywhere on the real plane. The work after that shows how to make the transition to circumstances appropriate to digital image processing, namely, an input of finite size and consisting of a finite number of quantized intensities. The last part of the paper contains applications of IDS to some “real-world” images and shows how it changes and improves the originals.
INTRODUCTION
The human visual system is very powerful and still far exceeds the performance of artificial visual systems in almost all areas. For this reason, researchers have used properties or models of how people see to improve image processing systems and in turn have studied those systems to learn more about human vision. The purpose of this paper is to apply a new model of part of the human visual system to the enhancement of digital images. The model, invented by Cornsweet [3], is known as IDS (an abbreviation for Intensity-Dependent Spread functions or Intensity-Dependent spatial Summation) and proposes a means by which cells in the retina interact. The mechanism of IDS is exceedingly simple, but it has been successful in explaining many phenomena of brightness perception. In image processing, the principal benefits of using IDS are that it enhances contrast, nonlinearly compresses the dynamic range, and enhances edges. IDS creates Mach bands at steps in intensities and the peak amplitude of this IDS output depends only on the ratio of the two intensities. To see the importance of this, note that a simplistic (but widely applicable) relationship between the intensity Z(X, y), illumination L(x, y), and reflectance R(n, y) is [7]
2.
IDS
Cornsweet and Yellott presented the mathematical description of IDS in [41. The input Z(X, y) is defined on the x-y plane, where x and y are the two independent spatial variables and Z(x, y) is positive (greater than zero) every-
m, Y> = ux, Y> * mx, Y>. 45
1049-%52/92 $3.00 Copyright 0 1992 by Academic hess, Inc. All rights of reproduction in any form reserved.
46
GREG REESE
where on the plane. Z(x, y) is the input, the input in?ensity, or the input image. IDS associates a spread function with every input point (p, 4) on the plane. This spread function has the follow-
ing properties: 1. It is nonnegative (greater than or equal to zero). 2. It has circular symmetry about the point (p, q). 3. Its volume is one. 4. No height of the spread function is greater than that at its center point (p, q) and this height is proportional to Z(p, q), the intensity at that center point. (Although Cornsweet and Yellott discuss a useful extension of IDS in which the height is proportional to some power of Z(p, q) [4, Section 51, this paper assumes the basic case of a height proportional to Z(p, q) itself.)
function is nonzero is the region of support.) Spread functions like Gaussian or exponential curves extend over the entire plane. The output from IDS at some point (x, y) is simply the integral of the values of all the spread functions at the point (x, y); that is, 0(x- Y) = I:= Ir, z (P, 4)WP,
+
(Y
qNx
- q)211dP 4.
- PI2 (2)
This is the fundamental equation of IDS. Using this equation, Cornsweet and Yellott derived the following properties of IDS: Rotation invariance. If the input is rotated by any amount, the output is unchanged except for a rotation of The third and fourth properties are crucial to understand- that same amount. ing spread functions and how they form the IDS output. Translation invariance. If the input is translated by any High intensities produce tall spread functions (by prop- amount, the output is unchanged except for a translation erty 4), but these spread functions must be narrow in of that same amount. order to keep the volume equal to one (by property 3). Uniform field. The output to any uniform field is one. Similarly, low intensities produce short, wide spread The property about uniform fields (inputs that are confunctions. stant everywhere) shows that IDS is a nonlinear system. A general expression for the spread function is Moreover, because of the finite (or essentially finite) radii of support of spread functions, large regions of constant (1) intensity are the same as uniform fields for output points Z(P, CIWNP, 4)[(x - Pj2 + (Y - q)211, well inside those regions. Thus, in general, wherever where S specifies a particular function to be used and x there are large patches of constant intensity in the input, and y are the coordinates at which the spread function is there will be regions in the output whose value is one, to be evaluated. Table 1 lists four common spread func- regardless of the values of the input intensities. tions. Spread functions like cylinders or cones are nonFigure 1 shows a graphical example of IDS in which zero only for a finite distance from their center-their the particular type of spread function is the Gaussian radius of support. (The region over which the spread listed in Table 1. The input is in Fig. la, the formation of IDS is in Fig. lb, and the output is in Fig. lc. (Although all inputs to IDS are defined on the plane, if the input pattern varies in only one direction, the output also varTABLE 1 ies only along that direction. All graphs in this paper that CommonIDS &m?adFunctions show only one independent variable actually represent a function of two variables that changes in only one direcFormula Name tion, taken to be along the x-axis.) The spread functions Gaussian on the right of the origin in Fig. lb are taller and narrower than those on the left because the input intensity on the right is greater than that on the left. Note that on both Exponential sides, far from the jump in intensities, the output is one, as the property about uniform fields predicts. 1 Figure la is also an example of a very important input, for r 5 m, 4) xhaTG3 step. A step Z(x) has the equation the Cylinder otherwise
,,,.“,, Cone
[l - ppr]
h-r<
d&
Z(x) =
Z, for x 2 0 { ZL for x < 0.
(3)
otherwise
0 r = d/(x - py +
(Y
- qY
For simplicity, the intensity on the left (IL) will be lower than the high intensity on the right (Z,).
IMAGE ENHANCEMENT
47
BY IDS
3. DISCRETE IDS
Intensity 2.5.. 21.5..
0.5--
-4
-2
Intensity
5
4
2
x
a
-5
and Spread
Functions
2.5
-2.5
The definition of IDS in Section 2 has both its input and output on the real plane. I refer to this form of the theory as “continuous IDS” and the version of IDS used for image processing as “discrete IDS.” The first part of this section illustrates the transition from continuous IDS to discrete IDS and the second part examines some of the computational requirements for doing IDS on a computer. 3.1 Transition to Discrete IDS
5
b IDlS, loutPut T
Let the input and output images consist of rectangular arrays of pixels. For simplicity, assume that the arrays are actually N x N squares, with N typically being 512. To relate the input array to continuous IDS, imagine placing it on the plane as shown in Fig. 2. The center of each pixel is on an integer point (i, j), where both i and j run from 1 to N. Each pixel’s side is one unit of length. Equation (2) gives the fundamental relationship between the input and the output in continuous IDS. For discrete IDS, the equation to use is w,
Y> = Ir,
+
/ym Z’(P, q)sv’(P,
4)b
- PI2
- s)21)4-J &,
(4)
I’(P,4) = cm, 4)
(5)
(Y
where C FIG. I. Cross section of IDS: (a) step input, (b) formation of output with Gaussian spread functions, and (c) IDS output to step. INPUT IMAGE
The step is important because it is a model, albeit a primitive one, of the boundary between two objects, and one of the fundamental goals of low-level vision is to fina such boundaries. Figure lc shows the IDS response to a step input. It has odd symmetry about the point (0, l), approaches one far from the origin, and has extrema whose values depend only on the ratio of the two intensities. Furthermore, the response is one exactly at the jump, less than one to the left, and greater than one to the right. There have been several demonstrations of image compression and reconstruction by IDS, with the methods based to a great extent on these properties of the step response 11, 2, 51. Another advantage of IDS is that it compresses the range of the input, but in a nonlinear fashion, so that certain important features of inputs, especially steps, are faithfully represented in the output. (Figure 5 is an illustration of this nonlinear compression.)
J N-
l
-
..32l-
I
I I 12
I
I
I
I
I
I
3
.
..N
I
t I
)X
FIG. 2. Input image superimposed on the x-y plane. Each pixel is centered on a point with integer coordinates.
48
GREG REESE
and c is a positive constant. A later part of this section gives the reason for using a scaled version of the original input Z(P, 4). The first step in the transition from continuous IDS to discrete IDS is to divide the entire real plane into unit squares, each centered on a point with integer coordinates. The next step is to rewrite the integral in Eq. (4) as
Z’(P,4UV’(P, 4Mx -PI2 + (Y - q)211&&, (6)
where m and n are integers. Equations (4) and (6) are identical but the integration in Eq. (6) takes place square by square. The first approximation in discrete IDS is due to the input having a finite extent. Equation (4) assumes that the input covers the entire plane, but in discrete IDS there is only an N x N section of the plane available. Thus for discrete IDS Eq. (6) becomes
Z’(P, q)W’(p, qW - p12+ (Y - q121) dp 4. (7) The output points near the border of the array will feel the effect of this truncated input most harshly because spread functions decrease with increasing distance and the value of the spread functions that would normally be in the rest of the plane is higher at the border of the output array than in its middle. Furthermore, since spread functions are always nonnegative, the effect at the edges will be one of decreasing the output because less of the spread functions is being added in. Thus the result of truncating the input is that there may be a slight darkening around the edge of the output image. The second approximation also involves the truncation of an infinite quantity. This time, it is the (possibly) infinite extent of the spread function. Although it is possible to have these spread functions reach all of the output pixels, it is impractical and unnecessary to do so. It is not practical because it sharply increases the computation time compared to that of spread functions with much smaller radii of support. It is not necessary because if the radius of support after truncation is large enough, the value of the spread function that is ignored is negligible and would contribute little if added to the output points. The finite extent of spread functions also affects the discrete IDS output near the border of an image. If the
largest radius of support is R,,,, output pixels within a distance of R,,, from an edge of the image may not have the same value as they would for continuous IDS because if the input image were made larger by adding a border of pixels, the spread functions from these added pixels might reach the output pixels near the original border and change their values. Conversely, output pixels that are farther than R,, from the edge of the image would remain unaffected by the addition of more input pixels. Just what do the input pixels represent? Figure 2 shows a grid of finite size overlaid on the x-y plane. Each square of the grid stands for 1 pixel in the input image. The intensity at any pixel was produced by a sensor that erased any variation of intensity over its own surface and replaced it with one averaged value. Since any intensity fluctuations that occurred within a pixel have been averaged out and are irretrievable, for purposes of working with IDS I define the original input (over the entire plane) as being composed of a grid of constant-intensity squares. Each input pixel is one of the squares of the grid and represents a square of constant input intensity. Thus if the input were to continuous IDS, it would still be defined everywhere on the entire plane but would be made up of squares of constant intensity. There are two advantages to this definition. The first is that I can use the mathematics for continuous IDS to calculate the output and compare it to that made by discrete IDS. This procedure is valid because the input to continuous IDS is defined to be a mosaic of squares that covers all points on the real plane. The second advantage is that many important intensity patterns, such as a uniform field, a step, or a square on a uniform background, can be broken down exactly into small squares. For these kinds of patterns, the discrete IDS output will match that of continuous IDS, permitting the use of two forms of IDS in a complementary fashion, with analytical results from continuous IDS verifying the output of discrete IDS and discrete IDS providing calculated values for the many input patterns that are too difficult to handle by theory. An N x N grid will also represent the output image and will be exactly aligned with the grid of the input image. As before (see Fig. 2), each square of the grid is centered on a point (i, j), where i and j are integers running from 1 to N. However, for the output image the value of the pixel at (i, j) will be the discrete IDS output at the point (i, j). One advantage of calculating the output at a point (i,j) is that continuous IDS also calculates the output at individual points. (To see this, replace (x, y) on the left side of Eq. (4) with (i, j).) Thus the definition of an output pixel allows a comparison of the output of continuous IDS to the output of discrete IDS. Another advantage of this definition is that it makes it easy to compute the
IMAGE ENHANCEMENT
output at points within a pixel. Such subpixel calculations are not normally important for qualitative enhancement, but they may matter in other applications of IDS, such as robotic vision or image reconstruction. There is another important difference between continuous IDS and discrete IDS. Consider the output at a point (i - 1, j) due to one of the four nearest input squares, say the square at (i, j). This output is given by the integral in the sum of Eq. (7), with x set to i - 1, y set toj, and Z’(p, 4) replaced by Zb, its constant value over the pixel at (i, j). This integral is j+ l/2 i j-
112 i
using cl@, y) as the input instead of just Z(X, y). For example, if the intensity is at its maximum value cZ,,~, the second equation in Table 2 implies that the minimum radius of support of the cylindrical spread function is
&in = q-~.
GS{Z[jNi
-
1 -
p12
+
(j
-
qj21}
dp
which is equivalent to the volume over the square at (i I ,j) and under a spread function whose intensity is Z; and is centered at (i, j). Now suppose that the radius of support of the spread function with intensity Z$ is less than or equal to one-half. This will cause the integrand to be zero everywhere except along the edge of the square at i - l/2, implying that the integral will also be zero. Since the entire spread function from the input at (i,j) is inside the square centered there and the volume under any spread function is one, the contribution to the discrete IDS output from that input is (as given by Eq. (8)) one for the output at (i, j) and zero elsewhere. If Zg were the minimum input intensity, all other intensities would also have a radius of support less than one-half, each output point would only feel the effect of the input square beneath it, and each output point would have a value of one, regardless of the input pattern. This deleterious effect is known as safuru~ion. It does not happen in continuous IDS because the output is calculated at every point on the plane, not just those points with integer coordinates. The way to avoid saturation is to set the scaling parameter c in Eq. (5) so that the minimum radius of support (which occurs in the spread function of the maximum intensity) is larger than one-half. This is the reason for TABLE 2 Properties of the Cylindrical Spread Function Name
Formula
I Spread function
0
I for r s 7% otherwise
Radius
1 r=TG
Position of peak in step response
C-C-
(8)
dq,
(9)
Combining Eq. (9) with Rmin> i, the requirement that eliminates saturation, yields
i+ 112 i-,,2
49
BY IDS
4 (10) 4mlx~
For &bit images, I,,
is 255, so c is bounded by
4 ' < 255~~' Different types of spread functions will have different values of c, but the method of determining the particular value is the same as that illustrated above. Eq. (9) shows that lowering c increases the radius of support, which in turn lengthens the computation time because the spread function is nonzero over more output pixels. The response of discrete IDS to a step shows additional effects of decreasing c. For a cylindrical spread function, Table 2 gives the position p’ of the maximum of the output as
where Z, is the sum of the two intensities. This equation shows that the location of the maximum of the edge response is proportional to the inverse of the square root of c. Decreasing c causes the maximum to appear farther from the jump in the two intensities. This broadens the edge, making it look wider and causing the outputs from two nearby edges to merge together more quickly. Figures 3 and 4 summarize the work up to this point. Figure 3 is a series of steps. The bars from left to right have intensities of 8, 16, 32, 64, 128, and 255. Figure 4 shows the discrete IDS output to this input and Figure 5 is a plot of one row in the middle of Fig. 4. (I have scaled the output in both of these figures to lie in the range
O-255.)
50
GREG REESE
FIG. 3. Input of five steps. Values of bars from left to right are 8, 16, 32, 64, 128, 255.
First, note that the large portions of each of the regions between all of the jumps in intensities have approximately the same output value. Second, note that the dynamic range has been nonlinearly compressed because the peak intensity of the output to each step is the same. Third, although the peak intensities are equal, an output peak appears closer to the jump as the intensities in a step increase. Thus for a fixed ratio of intensities, steps with higher intensities look crisper in the output than those with lower values. (An alternate interpretation of Fig. 3 is that the intensities of the steps are the same but the scaling factor doubles when moving from one step to the step on its right. This illustrates that the higher the scaling factor, the crisper the edge outputs.),Finally, note that parts of the border of the picture, especially near the lower input intensities, are dark, but as Fig. 4 shows, this phenomenon is not too intrusive. The location of the maximum closest to the origin, pki,, occurs for a step in which Is, the sum of the two intensities, is a maximum. (See the second equation in Table 2.) For 8-bit images, this happens when the low intensity is 254 and the high intensity is 255. The sum of the two, Ismor,is 509. Since maxima in discrete IDS appear only if they are at integer distances from the step, phi, must be made to be at least a distance of one from the step if it is to be present in the output. Eq. (12) and the requirement pki, 2 1 imply
2b0
360
rho
FIG. 5. Cross section of discrete IDS output in Fig. 4.
The second equation in Table 2 states that for the highest value of c in Eq. (13) (c = l/50971-)and the highest intensity (I = 255), the minimum radius of the spread function is 1.41. This is acceptable because it is greater than onehalf and thus does not cause saturation. Manipulation of Eq. (12) shows that the next detectable maximum (at a distance of 2 from the step) occurs when Is is 509/4, or about 127. This means that all steps whose intensities sum to more than 127 but less than 509 have their maximum between 1 and 2 units from the edge. These maxima are all undetectable because they’fall between the integer values where the output is calculated. (This illustrates more forcefully why one might want to do subpixel calculations for the output.) In particular, note that all steps whose low intensity is 127 or greater are in this category. However, just because the maxima of bright steps do not appear in the output does not mean that the steps are not emphasized at all. If a maximum falls fairly close to an integer point, the height at that point is almost the same as the maximum height. Table 3 illustrates this for Fig. 4. Although decreasing c makes the edge outputs wider, it also has a beneficial effect. It separates the locations of the peaks, because these locations are inversely proportional to the square root of c (see Eq. (12)). This is espe-
TABLE 3 Comparison between Peak Values in Continuous IDS and Discrete IDS for the Output to the Steps iu Fig. 3 Continuous Law intensity
a
FIG. 4. Discrete IDS output to Fig. 3.
16 32 64 128
High intensity 16 32 64 I28 255
Peak location 4.60 3.25 2.30 I .63 1.15
IDS Peak value I.108 1.108 I.108 I.108 I.108
Discrete Peak location 5 3 2 2 I
IDS Peak value
Percentage error
I.106 I.107 1.105 1.090 1.104
0.181 0.090 0.271 1.625 0.361
IMAGE
ENHANCEMENT
BY IDS
51
of such points within the region of support of any spread function as well as a finite number of intensities, a lookup table (LUT) can hold the values of the above integral for all intensities and all points within the region of support. To process an image, the program loads the precomputed LUT for the desired type of spread function. Then, given the input intensity at (m, n), it finds the output at (i, j) by referencing the LUT via that intensity and the distance between (m, n) and (i,j) and adds the output value given by the LUT to whatever is already stored as the output at (i, j). Repeating this procedure for each output pixel within reach of a given input pixel and for each input pixel in the original image creates the final discrete IDS output. Note that the actual discrete IDS calculation involves only additions, although implicitly there are calculations associated with determining the distance between two points (m, n) and (i, j) as well as many references to computer memory or disk storage. An additional benefit of using a LUT is that it lets the 1. The input to discrete IDS is not the entire plane. double integral in Eq. (7) be calculated ahead of time and 2. The output is a sampling of points of one region of so be numerically evaluated as accurately as desired the entire plane. without affecting the speed of the discrete IDS computa3. The regions of support of the spread functions may tion itself. A typical LUT size is that for a cylindrical have been truncated. spread function and a scaling parameter c 7 1/5097~, 4. The volume under the spread functions may not be which at 2374 4-byte entries is 9496 bytes long. exactly one. Even though the actual processing only involves addi5. The numerical integration needed to calculate the tion, it is still time consuming. It takes about 20 s to 2 min output may not be exact. (See Section 3.2 below.) to complete a discrete IDS run on an Apollo DNIOOOO 6. The real discrete IDS output is scaled to integers workstation, with the longer times being for the darker by a linear stretch, and there are attendant errors from images. In these trials, the input image was 512 x 512 quantization and the discarding of extreme values. pixels, 8 bits per pixel, and the output was the same size 7. Input intensities of zero are changed to one. but had 32 bits per pixel. This process can be done much faster in hardware. In fact, there is a commercially available hardware processor that accepts the same size input, 3.2 Computational Considerations yields a 512 x 512 x 16-bit output, and operates at 30 Suppose that the input is 512 x 512 pixels, the spread images per second [6]. function is the cylinder, and c is 1/5097r. Table 2 implies Discrete IDS would be wonderfully suited for parallel that the minimum radius is 1.41 and the maximum radius processing because all input pixels operate independently is 22.56. Because the radius and intensity are related,by of each other. Each processor would be very simple ber = m (see Table 2), for &bit images, most of the cause it only needs to determine the distance between radii are about 2. To get a rough feel for the amount of two points and add a value to an output point. This felicicomputation discrete IDS requires, suppose that each intous relationship between IDS and parallel processing is put pixel produces output at any output pixel within 2 a reflection of the theory’s original conception as a model distance units. Setting a circle of radius 2 on the’middle of of the retina in the eye. any square in the center of the grid in Fig. 2 shows that each input pixel infh.tences 21 output pixels. The calculation of the influence of an input pixel on an output pixel 4. RESULTS takes, place 21 x 5 12 x 5 12 times per image, or apgroximately 52 million times. Each of these calculations is the This section shows the results of using,discrete IDS on evaluation of the double integral in ‘Eq. (7). .This is a real and artificial images. All input images are 5 12 x 5 12 tremendous amount of computation. pixels with 8 bits per pixel, the spread function is the Fortunately, however, it is possible to completely cylinder, and unless otherwise noted, the scaling parameavoid calculating the above integral during the actual dis-d ter’ c is 115097f. Crete IDS computation. Since discrete IDS computes the Figure 6 shows six crosses, all of which are the same output only at integer points and there is a finite number size. The ratio of intensities of any cross to its back- _ cially true for steps of high intensity. In general, details in bright regions of the input will stand out more if c is decreased. Section 4 demonstrates these effects with a series of images processed with different values of c. The last step in the transition from continuous IDS to discrete IDS is that the output of discrete IDS, a collection of real numbers (often between 0.8 and 1.2), must be converted to integers for use with standard image processing software. The method of scaling that I use is the standard linear stretch [7] between some high and low points in the discrete IDS output (e.g., the 2 and 98% points in the cumulative histogram) to the g-bit minimum and maximum display values. In sum, the procedure in this section for doing discrete IDS makes the processed images come out very close to what the analytic output would be if it were calculated with continuous IDS. Nevertheless, there are a number of differences that arise from the following causes:
52
GREG REESE
approximately the same gray level- 127. They are not all equal to exactly 127 for the reasons explained at the end of Section 3.1, but they do come quite close, in Fig. 7 varying from 126 to 128. The more important phenomenon is the distortion in the output due to the different intensities of the six crosses and their backgrounds. Even though the ratio of the cross to its background is the same in all six subimages, the output is crisper when the absolute values of those intensities are higher. Note that edges in both the vertical and the horizontal directions are enhanced-enhanced in the sense that they are now perturbations from a uniform background. In fact, because IDS is rotationand translation-invariant (see Section 2), it enhances edges running in any direction and at any location. All of the edges display Mach bands, with the bands being narrower (though of the same peak amplitude) for the crosses with higher intensities. Although Fig. 7 is the discrete IDS output to an artificial input, it demonstrates a number of phenomena that occur in outputs to images of real scenes: FIG. 6. Six crosses. Values of crosses in the top row from left to right are 4, 8, and 16. Values of crosses in the bottom row from left to right are 32, 64, and 128. Backgrounds are one-half of cross values.
ground is 2 : 1. The crosses in the top row have intensities 4, 8, and 16 (from left to right) and those in the bottom row have intensities 32, 64, and 128. Figure 7 is the discrete IDS output to the image in Fig. 6. Note that the uniform areas far from changes in intensities come out at
FIG. 7. IDS output to crosses in Fig. 6.
Uniform regions are set to the same constant output. Mach bands occur at edges. For a fixed ratio of intensities, edges with lower intensities produce wider Mach bands, even though the peak heights of the bands stay the same. Similar figures in the input scene emerge fairly alike in the output, but the figures with lower intensities are distorted more than those with higher intensities. l l
l
Figure 8, a more realistic image, is a chest X ray and, as often happens in these kinds of pictures, the spinal column does not come out well, appearing as a vertical white region. Figure 9 is the output of discrete IDS to
FIG. 8. Chest X ray.
IMAGE ENHANCEMENT
53
FIG. 11. IDS output to X ray in Fig. 8 with c = 1115OOm.
FIG. 9. IDS output to X ray in Fig. 8 with c = 1/509~.
Fig. 8. A number of features of the original image are more evident in the processed output:
BY IDS
Figure 9 shows the discrete IDS output with the scaling parameter c set to l/50971-.The next three figures show the effects of decreasing c. Figures 10, 11, and 12 are the outputs with c decreased to 1/10007r, 1/1500~, and l/200077. Some interesting characteristics of this series of images are:
The edges of the ribs stand out more. Various lines in the picture are now seen to actually cross the spinal column, especially the line running from the top right of the image to the left side of the middle of the spinal column. The large, dark, oval blob in the upper part of the right lung becomes apparent. There is now structure in the spinal column, namely an alternating series of dark and light inverted “v” shapes. These shapes disappear in the set of images shown next.
As c decreases, the edges become wider. As c decreases, more detail appears in the regions that were white in the original image. In particular, note that the inverted v’s in the region of the spinal column disappear and the edges of the column, especially the left edge, become visible. The oval blob shows up more clearly in all of the
FIG. 10. IDS output to X ray in Fig. 8 with c = l/lOOOn.
FIG. 12. IDS output to X ray in Fig. 8 with c = 1/2OOOn.
l
l
l
l
l l
l
54
GREG REESE
processed images than in the original. Note that discrete IDS simultaneously enhances the input in dark areas (the upper left and right) and in white areas (the spinal column). The images in Fig. 9 through Fig. 12 demonstrate that as the scaling parameter ‘c decreases, details in high-intensity regions become more visible, but at the same time edges become wider. Thus the user makes a particular choiCe of c based on a trade-off between having details in high-intensity areas or having crisp edges. In all cases, c is the only parameter in discrete IDS, so there is no need to worry about trying to relate and fine-tune a plethora of parameters. Figure 13 is a satellite photograph of an area in the arctic region. The jagged, vertical lines descending from the top of the picture are ice leads, fractures or passageways through sea ice. Unfortunately, the leads are often covered by clouds or haze and are difficult to see in unenhanced images. Figure 14 is the discrete IDS output of Fig. 13. In the processed image: The large leads are more visible and they extend further down the picture. The small leads branching off of the large leads are now clearer. . There is a long horizontal line in the middle of the discrete IDS image that is barely visible in the original picture. This line is a horizontal continuation of one of the vertical leads on the left side of the picture. Details in the lower left stand out more. l
FIG. 14. IDS output to ice leads in Fig. 13.
l
l
Discrete IDS does a good job of looking through hazy clouds and seeing the underlying scene, in this case the ice leads. It works well because the haze acts like a fairly uniform shadow and as Fig. 6 and Fig. 7 demonstrate, discrete IDS is fairly insensitive to uniform changes in illumination. However, if the cloud cover is too thick to see through, no amount of enhancement will bring out the scene beneath the clouds. In this situation, discrete IDS will enhance the clouds themselves (see the middle and top right of Fig. 13) and such improvements may be useful to meteorologists. 5.
CONCLUSION
IDS is a theory that is a good model of how people see patterns of light and dark, but it is also a good tool for improving digital images. The version of IDS for the computer, discrete IDS, enhances the edges, nonlinearly compresses the dynamic range, and automatically adjusts itself to the local intensities of an image, all without the use of windows or statistics. IDS is particularly useful for bringing out detail in regions of images that have low contrast. If discrete IDS is calculated with a look-up table, the only computations involved are additions and references to the table. Discrete IDS would be very well suited for operation on a parallel computer. ACKNOWLEDGMENT FIG. 13. Ice leads in the Lincoln Sea.
My thanks to Rich Goldkamp for his help in acquiring images and patiently shuttling them among various computer systems.
IMAGE ENHANCEMENT
REFERENCES 1. R. Alter-Gartenberg, F. 0. Huck, and R. Narayanswamy, Image recovery from edge primitives, J Opt. Sot. Am. Ser. A 7(5), May 1990, 898-911. 2. R. Alter-Gartenberg and R. Narayanswamy, Image coding by edge primitives, in Proceedings, Visual Communications and Image Processing ‘90, SPIE, Vol. 1360, pp. 1608-1619, Laussanne, Switzerland, 1990. 3. T. N. Cornsweet, Prentice award lecture: A simple retinal mechanism that has complex and profound effects on perception, Am. .J. Opt. Physiol.
Opt. 62(7), 1985, 427438.
4. T. N. Cornsweet and J. I. Yellott, Jr., Intensity-dependent spatial summation, J. Opt. Sot. Am. Ser. A 2(10), Oct. 1985, 1769-1786.
BY IDS
55
5. R. Narayanswamy, R. Alter-Gartenberg, and F. 0. Huck. Robust image coding with a model of adaptive retinal processing, in Optics, Illumination and Image Sensing for Machine Vision V, SPIE, Vol. 1385, pp. 83-103, Boston, MA, 1991. 6. AC-200 Series Adaptive Imager Products, Odetics Artificial Intelligence Center, 1515 S. Manchester Ave., Anaheim, CA 928022907.
7. A. Rosenfeld and A. C. Kak, Digital Picture Processing, Academic Press, Orlando, 1982. 8. V. T. Tom, Adaptive filter techniques for digital image enhancement, in Proceedings, SPIE, Digital Image Processing, 198.5, Vol. 528, pp. 29-42.
GRAPHICAL
MODELS
Vol. 54, No. 1, January,
AND
IMAGE
PROCESSING
pp. 45-55, 1992
Image Enhancement by Intensity-Dependent
Spread Functions
GREG REESE Applied Math and Physics Section, Sverdrup Technology,
Incorporated,
Stennis Space Center, Mississippi
39529
Received August 27, 1990; accepted September 10, 1991
If the illumination on a reflectance edge is constant, that illumination cancels out of the ratio of the two intensities, and the peak height of the IDS response to that edge will always be the same. This property implies that IDS can aid in detecting a physical characteristic of objects in a scene-their reflectance. Thus if the same object is present several times in one image under very different constant illuminations, for example, once in sunlight and once in shadow, IDS transforms both objects so that their outputs look quite similar to each other. As some later pictures show, the outputs are not exactly the same be-
This paper describes a method of enhancing images based on a new model of human vision. The model is an adaptive, nonlinear system that enhances edges, nonlinearly compresses dynamic ranges, automatically adjusts itself to local variations in intensity, can (but need not) be carried out on a parallel processor, and has only one adjustable parameter that requires no intricate tuning. It is particularly useful in images that have regions of low contrast, such as brightly lit areas and those covered by shadows. The text explains how to convert the original mathematical theory to one useful for digital image processing and illustrates the output of the system to a variety of images. The focus is on qualitative enhancement of images, but the system would also be useful in machine vision. 0 1992 Academic Press, Inc.
1.
cause images
of objects
get more smeared
as the intensi-
ties get lower. Nevertheless, the result can still be a substantial improvement over the unenhanced image. The IDS output at a point depends on the input intensities near that point and the system adapts itself to changes in these intensities. There have been numerous studies of adaptive systems (for a review see [8]). Although the performance of adaptive systems is often substantially better than that of their nonadaptive counterparts, they do suffer from a number of drawbacks, in general, being complex, computationally expensive, and subject to the vagaries of window size and shape. Since IDS has no windows it does not suffer from the latter problems and although it requires many calculations, the computations are very simple. The next section of the paper is a review of the mathematical formulation of IDS as a system that operates on inputs defined everywhere on the real plane. The work after that shows how to make the transition to circumstances appropriate to digital image processing, namely, an input of finite size and consisting of a finite number of quantized intensities. The last part of the paper contains applications of IDS to some “real-world” images and shows how it changes and improves the originals.
INTRODUCTION
The human visual system is very powerful and still far exceeds the performance of artificial visual systems in almost all areas. For this reason, researchers have used properties or models of how people see to improve image processing systems and in turn have studied those systems to learn more about human vision. The purpose of this paper is to apply a new model of part of the human visual system to the enhancement of digital images. The model, invented by Cornsweet [3], is known as IDS (an abbreviation for Intensity-Dependent Spread functions or Intensity-Dependent spatial Summation) and proposes a means by which cells in the retina interact. The mechanism of IDS is exceedingly simple, but it has been successful in explaining many phenomena of brightness perception. In image processing, the principal benefits of using IDS are that it enhances contrast, nonlinearly compresses the dynamic range, and enhances edges. IDS creates Mach bands at steps in intensities and the peak amplitude of this IDS output depends only on the ratio of the two intensities. To see the importance of this, note that a simplistic (but widely applicable) relationship between the intensity Z(X, y), illumination L(x, y), and reflectance R(n, y) is [7]
2.
IDS
Cornsweet and Yellott presented the mathematical description of IDS in [41. The input Z(X, y) is defined on the x-y plane, where x and y are the two independent spatial variables and Z(x, y) is positive (greater than zero) every-
m, Y> = ux, Y> * mx, Y>. 45
1049-%52/92 $3.00 Copyright 0 1992 by Academic hess, Inc. All rights of reproduction in any form reserved.
46
GREG REESE
where on the plane. Z(x, y) is the input, the input in?ensity, or the input image. IDS associates a spread function with every input point (p, 4) on the plane. This spread function has the follow-
ing properties: 1. It is nonnegative (greater than or equal to zero). 2. It has circular symmetry about the point (p, q). 3. Its volume is one. 4. No height of the spread function is greater than that at its center point (p, q) and this height is proportional to Z(p, q), the intensity at that center point. (Although Cornsweet and Yellott discuss a useful extension of IDS in which the height is proportional to some power of Z(p, q) [4, Section 51, this paper assumes the basic case of a height proportional to Z(p, q) itself.)
function is nonzero is the region of support.) Spread functions like Gaussian or exponential curves extend over the entire plane. The output from IDS at some point (x, y) is simply the integral of the values of all the spread functions at the point (x, y); that is, 0(x- Y) = I:= Ir, z (P, 4)WP,
+
(Y
qNx
- q)211dP 4.
- PI2 (2)
This is the fundamental equation of IDS. Using this equation, Cornsweet and Yellott derived the following properties of IDS: Rotation invariance. If the input is rotated by any amount, the output is unchanged except for a rotation of The third and fourth properties are crucial to understand- that same amount. ing spread functions and how they form the IDS output. Translation invariance. If the input is translated by any High intensities produce tall spread functions (by prop- amount, the output is unchanged except for a translation erty 4), but these spread functions must be narrow in of that same amount. order to keep the volume equal to one (by property 3). Uniform field. The output to any uniform field is one. Similarly, low intensities produce short, wide spread The property about uniform fields (inputs that are confunctions. stant everywhere) shows that IDS is a nonlinear system. A general expression for the spread function is Moreover, because of the finite (or essentially finite) radii of support of spread functions, large regions of constant (1) intensity are the same as uniform fields for output points Z(P, CIWNP, 4)[(x - Pj2 + (Y - q)211, well inside those regions. Thus, in general, wherever where S specifies a particular function to be used and x there are large patches of constant intensity in the input, and y are the coordinates at which the spread function is there will be regions in the output whose value is one, to be evaluated. Table 1 lists four common spread func- regardless of the values of the input intensities. tions. Spread functions like cylinders or cones are nonFigure 1 shows a graphical example of IDS in which zero only for a finite distance from their center-their the particular type of spread function is the Gaussian radius of support. (The region over which the spread listed in Table 1. The input is in Fig. la, the formation of IDS is in Fig. lb, and the output is in Fig. lc. (Although all inputs to IDS are defined on the plane, if the input pattern varies in only one direction, the output also varTABLE 1 ies only along that direction. All graphs in this paper that CommonIDS &m?adFunctions show only one independent variable actually represent a function of two variables that changes in only one direcFormula Name tion, taken to be along the x-axis.) The spread functions Gaussian on the right of the origin in Fig. lb are taller and narrower than those on the left because the input intensity on the right is greater than that on the left. Note that on both Exponential sides, far from the jump in intensities, the output is one, as the property about uniform fields predicts. 1 Figure la is also an example of a very important input, for r 5 m, 4) xhaTG3 step. A step Z(x) has the equation the Cylinder otherwise
,,,.“,, Cone
[l - ppr]
h-r<
d&
Z(x) =
Z, for x 2 0 { ZL for x < 0.
(3)
otherwise
0 r = d/(x - py +
(Y
- qY
For simplicity, the intensity on the left (IL) will be lower than the high intensity on the right (Z,).
IMAGE ENHANCEMENT
47
BY IDS
3. DISCRETE IDS
Intensity 2.5.. 21.5..
0.5--
-4
-2
Intensity
5
4
2
x
a
-5
and Spread
Functions
2.5
-2.5
The definition of IDS in Section 2 has both its input and output on the real plane. I refer to this form of the theory as “continuous IDS” and the version of IDS used for image processing as “discrete IDS.” The first part of this section illustrates the transition from continuous IDS to discrete IDS and the second part examines some of the computational requirements for doing IDS on a computer. 3.1 Transition to Discrete IDS
5
b IDlS, loutPut T
Let the input and output images consist of rectangular arrays of pixels. For simplicity, assume that the arrays are actually N x N squares, with N typically being 512. To relate the input array to continuous IDS, imagine placing it on the plane as shown in Fig. 2. The center of each pixel is on an integer point (i, j), where both i and j run from 1 to N. Each pixel’s side is one unit of length. Equation (2) gives the fundamental relationship between the input and the output in continuous IDS. For discrete IDS, the equation to use is w,
Y> = Ir,
+
/ym Z’(P, q)sv’(P,
4)b
- PI2
- s)21)4-J &,
(4)
I’(P,4) = cm, 4)
(5)
(Y
where C FIG. I. Cross section of IDS: (a) step input, (b) formation of output with Gaussian spread functions, and (c) IDS output to step. INPUT IMAGE
The step is important because it is a model, albeit a primitive one, of the boundary between two objects, and one of the fundamental goals of low-level vision is to fina such boundaries. Figure lc shows the IDS response to a step input. It has odd symmetry about the point (0, l), approaches one far from the origin, and has extrema whose values depend only on the ratio of the two intensities. Furthermore, the response is one exactly at the jump, less than one to the left, and greater than one to the right. There have been several demonstrations of image compression and reconstruction by IDS, with the methods based to a great extent on these properties of the step response 11, 2, 51. Another advantage of IDS is that it compresses the range of the input, but in a nonlinear fashion, so that certain important features of inputs, especially steps, are faithfully represented in the output. (Figure 5 is an illustration of this nonlinear compression.)
J N-
l
-
..32l-
I
I I 12
I
I
I
I
I
I
3
.
..N
I
t I
)X
FIG. 2. Input image superimposed on the x-y plane. Each pixel is centered on a point with integer coordinates.
48
GREG REESE
and c is a positive constant. A later part of this section gives the reason for using a scaled version of the original input Z(P, 4). The first step in the transition from continuous IDS to discrete IDS is to divide the entire real plane into unit squares, each centered on a point with integer coordinates. The next step is to rewrite the integral in Eq. (4) as
Z’(P,4UV’(P, 4Mx -PI2 + (Y - q)211&&, (6)
where m and n are integers. Equations (4) and (6) are identical but the integration in Eq. (6) takes place square by square. The first approximation in discrete IDS is due to the input having a finite extent. Equation (4) assumes that the input covers the entire plane, but in discrete IDS there is only an N x N section of the plane available. Thus for discrete IDS Eq. (6) becomes
Z’(P, q)W’(p, qW - p12+ (Y - q121) dp 4. (7) The output points near the border of the array will feel the effect of this truncated input most harshly because spread functions decrease with increasing distance and the value of the spread functions that would normally be in the rest of the plane is higher at the border of the output array than in its middle. Furthermore, since spread functions are always nonnegative, the effect at the edges will be one of decreasing the output because less of the spread functions is being added in. Thus the result of truncating the input is that there may be a slight darkening around the edge of the output image. The second approximation also involves the truncation of an infinite quantity. This time, it is the (possibly) infinite extent of the spread function. Although it is possible to have these spread functions reach all of the output pixels, it is impractical and unnecessary to do so. It is not practical because it sharply increases the computation time compared to that of spread functions with much smaller radii of support. It is not necessary because if the radius of support after truncation is large enough, the value of the spread function that is ignored is negligible and would contribute little if added to the output points. The finite extent of spread functions also affects the discrete IDS output near the border of an image. If the
largest radius of support is R,,,, output pixels within a distance of R,,, from an edge of the image may not have the same value as they would for continuous IDS because if the input image were made larger by adding a border of pixels, the spread functions from these added pixels might reach the output pixels near the original border and change their values. Conversely, output pixels that are farther than R,, from the edge of the image would remain unaffected by the addition of more input pixels. Just what do the input pixels represent? Figure 2 shows a grid of finite size overlaid on the x-y plane. Each square of the grid stands for 1 pixel in the input image. The intensity at any pixel was produced by a sensor that erased any variation of intensity over its own surface and replaced it with one averaged value. Since any intensity fluctuations that occurred within a pixel have been averaged out and are irretrievable, for purposes of working with IDS I define the original input (over the entire plane) as being composed of a grid of constant-intensity squares. Each input pixel is one of the squares of the grid and represents a square of constant input intensity. Thus if the input were to continuous IDS, it would still be defined everywhere on the entire plane but would be made up of squares of constant intensity. There are two advantages to this definition. The first is that I can use the mathematics for continuous IDS to calculate the output and compare it to that made by discrete IDS. This procedure is valid because the input to continuous IDS is defined to be a mosaic of squares that covers all points on the real plane. The second advantage is that many important intensity patterns, such as a uniform field, a step, or a square on a uniform background, can be broken down exactly into small squares. For these kinds of patterns, the discrete IDS output will match that of continuous IDS, permitting the use of two forms of IDS in a complementary fashion, with analytical results from continuous IDS verifying the output of discrete IDS and discrete IDS providing calculated values for the many input patterns that are too difficult to handle by theory. An N x N grid will also represent the output image and will be exactly aligned with the grid of the input image. As before (see Fig. 2), each square of the grid is centered on a point (i, j), where i and j are integers running from 1 to N. However, for the output image the value of the pixel at (i, j) will be the discrete IDS output at the point (i, j). One advantage of calculating the output at a point (i,j) is that continuous IDS also calculates the output at individual points. (To see this, replace (x, y) on the left side of Eq. (4) with (i, j).) Thus the definition of an output pixel allows a comparison of the output of continuous IDS to the output of discrete IDS. Another advantage of this definition is that it makes it easy to compute the
IMAGE ENHANCEMENT
output at points within a pixel. Such subpixel calculations are not normally important for qualitative enhancement, but they may matter in other applications of IDS, such as robotic vision or image reconstruction. There is another important difference between continuous IDS and discrete IDS. Consider the output at a point (i - 1, j) due to one of the four nearest input squares, say the square at (i, j). This output is given by the integral in the sum of Eq. (7), with x set to i - 1, y set toj, and Z’(p, 4) replaced by Zb, its constant value over the pixel at (i, j). This integral is j+ l/2 i j-
112 i
using cl@, y) as the input instead of just Z(X, y). For example, if the intensity is at its maximum value cZ,,~, the second equation in Table 2 implies that the minimum radius of support of the cylindrical spread function is
&in = q-~.
GS{Z[jNi
-
1 -
p12
+
(j
-
qj21}
dp
which is equivalent to the volume over the square at (i I ,j) and under a spread function whose intensity is Z; and is centered at (i, j). Now suppose that the radius of support of the spread function with intensity Z$ is less than or equal to one-half. This will cause the integrand to be zero everywhere except along the edge of the square at i - l/2, implying that the integral will also be zero. Since the entire spread function from the input at (i,j) is inside the square centered there and the volume under any spread function is one, the contribution to the discrete IDS output from that input is (as given by Eq. (8)) one for the output at (i, j) and zero elsewhere. If Zg were the minimum input intensity, all other intensities would also have a radius of support less than one-half, each output point would only feel the effect of the input square beneath it, and each output point would have a value of one, regardless of the input pattern. This deleterious effect is known as safuru~ion. It does not happen in continuous IDS because the output is calculated at every point on the plane, not just those points with integer coordinates. The way to avoid saturation is to set the scaling parameter c in Eq. (5) so that the minimum radius of support (which occurs in the spread function of the maximum intensity) is larger than one-half. This is the reason for TABLE 2 Properties of the Cylindrical Spread Function Name
Formula
I Spread function
0
I for r s 7% otherwise
Radius
1 r=TG
Position of peak in step response
C-C-
(8)
dq,
(9)
Combining Eq. (9) with Rmin> i, the requirement that eliminates saturation, yields
i+ 112 i-,,2
49
BY IDS
4 (10) 4mlx~
For &bit images, I,,
is 255, so c is bounded by
4 ' < 255~~' Different types of spread functions will have different values of c, but the method of determining the particular value is the same as that illustrated above. Eq. (9) shows that lowering c increases the radius of support, which in turn lengthens the computation time because the spread function is nonzero over more output pixels. The response of discrete IDS to a step shows additional effects of decreasing c. For a cylindrical spread function, Table 2 gives the position p’ of the maximum of the output as
where Z, is the sum of the two intensities. This equation shows that the location of the maximum of the edge response is proportional to the inverse of the square root of c. Decreasing c causes the maximum to appear farther from the jump in the two intensities. This broadens the edge, making it look wider and causing the outputs from two nearby edges to merge together more quickly. Figures 3 and 4 summarize the work up to this point. Figure 3 is a series of steps. The bars from left to right have intensities of 8, 16, 32, 64, 128, and 255. Figure 4 shows the discrete IDS output to this input and Figure 5 is a plot of one row in the middle of Fig. 4. (I have scaled the output in both of these figures to lie in the range
O-255.)
50
GREG REESE
FIG. 3. Input of five steps. Values of bars from left to right are 8, 16, 32, 64, 128, 255.
First, note that the large portions of each of the regions between all of the jumps in intensities have approximately the same output value. Second, note that the dynamic range has been nonlinearly compressed because the peak intensity of the output to each step is the same. Third, although the peak intensities are equal, an output peak appears closer to the jump as the intensities in a step increase. Thus for a fixed ratio of intensities, steps with higher intensities look crisper in the output than those with lower values. (An alternate interpretation of Fig. 3 is that the intensities of the steps are the same but the scaling factor doubles when moving from one step to the step on its right. This illustrates that the higher the scaling factor, the crisper the edge outputs.),Finally, note that parts of the border of the picture, especially near the lower input intensities, are dark, but as Fig. 4 shows, this phenomenon is not too intrusive. The location of the maximum closest to the origin, pki,, occurs for a step in which Is, the sum of the two intensities, is a maximum. (See the second equation in Table 2.) For 8-bit images, this happens when the low intensity is 254 and the high intensity is 255. The sum of the two, Ismor,is 509. Since maxima in discrete IDS appear only if they are at integer distances from the step, phi, must be made to be at least a distance of one from the step if it is to be present in the output. Eq. (12) and the requirement pki, 2 1 imply
2b0
360
rho
FIG. 5. Cross section of discrete IDS output in Fig. 4.
The second equation in Table 2 states that for the highest value of c in Eq. (13) (c = l/50971-)and the highest intensity (I = 255), the minimum radius of the spread function is 1.41. This is acceptable because it is greater than onehalf and thus does not cause saturation. Manipulation of Eq. (12) shows that the next detectable maximum (at a distance of 2 from the step) occurs when Is is 509/4, or about 127. This means that all steps whose intensities sum to more than 127 but less than 509 have their maximum between 1 and 2 units from the edge. These maxima are all undetectable because they’fall between the integer values where the output is calculated. (This illustrates more forcefully why one might want to do subpixel calculations for the output.) In particular, note that all steps whose low intensity is 127 or greater are in this category. However, just because the maxima of bright steps do not appear in the output does not mean that the steps are not emphasized at all. If a maximum falls fairly close to an integer point, the height at that point is almost the same as the maximum height. Table 3 illustrates this for Fig. 4. Although decreasing c makes the edge outputs wider, it also has a beneficial effect. It separates the locations of the peaks, because these locations are inversely proportional to the square root of c (see Eq. (12)). This is espe-
TABLE 3 Comparison between Peak Values in Continuous IDS and Discrete IDS for the Output to the Steps iu Fig. 3 Continuous Law intensity
a
FIG. 4. Discrete IDS output to Fig. 3.
16 32 64 128
High intensity 16 32 64 I28 255
Peak location 4.60 3.25 2.30 I .63 1.15
IDS Peak value I.108 1.108 I.108 I.108 I.108
Discrete Peak location 5 3 2 2 I
IDS Peak value
Percentage error
I.106 I.107 1.105 1.090 1.104
0.181 0.090 0.271 1.625 0.361
IMAGE
ENHANCEMENT
BY IDS
51
of such points within the region of support of any spread function as well as a finite number of intensities, a lookup table (LUT) can hold the values of the above integral for all intensities and all points within the region of support. To process an image, the program loads the precomputed LUT for the desired type of spread function. Then, given the input intensity at (m, n), it finds the output at (i, j) by referencing the LUT via that intensity and the distance between (m, n) and (i,j) and adds the output value given by the LUT to whatever is already stored as the output at (i, j). Repeating this procedure for each output pixel within reach of a given input pixel and for each input pixel in the original image creates the final discrete IDS output. Note that the actual discrete IDS calculation involves only additions, although implicitly there are calculations associated with determining the distance between two points (m, n) and (i, j) as well as many references to computer memory or disk storage. An additional benefit of using a LUT is that it lets the 1. The input to discrete IDS is not the entire plane. double integral in Eq. (7) be calculated ahead of time and 2. The output is a sampling of points of one region of so be numerically evaluated as accurately as desired the entire plane. without affecting the speed of the discrete IDS computa3. The regions of support of the spread functions may tion itself. A typical LUT size is that for a cylindrical have been truncated. spread function and a scaling parameter c 7 1/5097~, 4. The volume under the spread functions may not be which at 2374 4-byte entries is 9496 bytes long. exactly one. Even though the actual processing only involves addi5. The numerical integration needed to calculate the tion, it is still time consuming. It takes about 20 s to 2 min output may not be exact. (See Section 3.2 below.) to complete a discrete IDS run on an Apollo DNIOOOO 6. The real discrete IDS output is scaled to integers workstation, with the longer times being for the darker by a linear stretch, and there are attendant errors from images. In these trials, the input image was 512 x 512 quantization and the discarding of extreme values. pixels, 8 bits per pixel, and the output was the same size 7. Input intensities of zero are changed to one. but had 32 bits per pixel. This process can be done much faster in hardware. In fact, there is a commercially available hardware processor that accepts the same size input, 3.2 Computational Considerations yields a 512 x 512 x 16-bit output, and operates at 30 Suppose that the input is 512 x 512 pixels, the spread images per second [6]. function is the cylinder, and c is 1/5097r. Table 2 implies Discrete IDS would be wonderfully suited for parallel that the minimum radius is 1.41 and the maximum radius processing because all input pixels operate independently is 22.56. Because the radius and intensity are related,by of each other. Each processor would be very simple ber = m (see Table 2), for &bit images, most of the cause it only needs to determine the distance between radii are about 2. To get a rough feel for the amount of two points and add a value to an output point. This felicicomputation discrete IDS requires, suppose that each intous relationship between IDS and parallel processing is put pixel produces output at any output pixel within 2 a reflection of the theory’s original conception as a model distance units. Setting a circle of radius 2 on the’middle of of the retina in the eye. any square in the center of the grid in Fig. 2 shows that each input pixel infh.tences 21 output pixels. The calculation of the influence of an input pixel on an output pixel 4. RESULTS takes, place 21 x 5 12 x 5 12 times per image, or apgroximately 52 million times. Each of these calculations is the This section shows the results of using,discrete IDS on evaluation of the double integral in ‘Eq. (7). .This is a real and artificial images. All input images are 5 12 x 5 12 tremendous amount of computation. pixels with 8 bits per pixel, the spread function is the Fortunately, however, it is possible to completely cylinder, and unless otherwise noted, the scaling parameavoid calculating the above integral during the actual dis-d ter’ c is 115097f. Crete IDS computation. Since discrete IDS computes the Figure 6 shows six crosses, all of which are the same output only at integer points and there is a finite number size. The ratio of intensities of any cross to its back- _ cially true for steps of high intensity. In general, details in bright regions of the input will stand out more if c is decreased. Section 4 demonstrates these effects with a series of images processed with different values of c. The last step in the transition from continuous IDS to discrete IDS is that the output of discrete IDS, a collection of real numbers (often between 0.8 and 1.2), must be converted to integers for use with standard image processing software. The method of scaling that I use is the standard linear stretch [7] between some high and low points in the discrete IDS output (e.g., the 2 and 98% points in the cumulative histogram) to the g-bit minimum and maximum display values. In sum, the procedure in this section for doing discrete IDS makes the processed images come out very close to what the analytic output would be if it were calculated with continuous IDS. Nevertheless, there are a number of differences that arise from the following causes:
52
GREG REESE
approximately the same gray level- 127. They are not all equal to exactly 127 for the reasons explained at the end of Section 3.1, but they do come quite close, in Fig. 7 varying from 126 to 128. The more important phenomenon is the distortion in the output due to the different intensities of the six crosses and their backgrounds. Even though the ratio of the cross to its background is the same in all six subimages, the output is crisper when the absolute values of those intensities are higher. Note that edges in both the vertical and the horizontal directions are enhanced-enhanced in the sense that they are now perturbations from a uniform background. In fact, because IDS is rotationand translation-invariant (see Section 2), it enhances edges running in any direction and at any location. All of the edges display Mach bands, with the bands being narrower (though of the same peak amplitude) for the crosses with higher intensities. Although Fig. 7 is the discrete IDS output to an artificial input, it demonstrates a number of phenomena that occur in outputs to images of real scenes: FIG. 6. Six crosses. Values of crosses in the top row from left to right are 4, 8, and 16. Values of crosses in the bottom row from left to right are 32, 64, and 128. Backgrounds are one-half of cross values.
ground is 2 : 1. The crosses in the top row have intensities 4, 8, and 16 (from left to right) and those in the bottom row have intensities 32, 64, and 128. Figure 7 is the discrete IDS output to the image in Fig. 6. Note that the uniform areas far from changes in intensities come out at
FIG. 7. IDS output to crosses in Fig. 6.
Uniform regions are set to the same constant output. Mach bands occur at edges. For a fixed ratio of intensities, edges with lower intensities produce wider Mach bands, even though the peak heights of the bands stay the same. Similar figures in the input scene emerge fairly alike in the output, but the figures with lower intensities are distorted more than those with higher intensities. l l
l
Figure 8, a more realistic image, is a chest X ray and, as often happens in these kinds of pictures, the spinal column does not come out well, appearing as a vertical white region. Figure 9 is the output of discrete IDS to
FIG. 8. Chest X ray.
IMAGE ENHANCEMENT
53
FIG. 11. IDS output to X ray in Fig. 8 with c = 1115OOm.
FIG. 9. IDS output to X ray in Fig. 8 with c = 1/509~.
Fig. 8. A number of features of the original image are more evident in the processed output:
BY IDS
Figure 9 shows the discrete IDS output with the scaling parameter c set to l/50971-.The next three figures show the effects of decreasing c. Figures 10, 11, and 12 are the outputs with c decreased to 1/10007r, 1/1500~, and l/200077. Some interesting characteristics of this series of images are:
The edges of the ribs stand out more. Various lines in the picture are now seen to actually cross the spinal column, especially the line running from the top right of the image to the left side of the middle of the spinal column. The large, dark, oval blob in the upper part of the right lung becomes apparent. There is now structure in the spinal column, namely an alternating series of dark and light inverted “v” shapes. These shapes disappear in the set of images shown next.
As c decreases, the edges become wider. As c decreases, more detail appears in the regions that were white in the original image. In particular, note that the inverted v’s in the region of the spinal column disappear and the edges of the column, especially the left edge, become visible. The oval blob shows up more clearly in all of the
FIG. 10. IDS output to X ray in Fig. 8 with c = l/lOOOn.
FIG. 12. IDS output to X ray in Fig. 8 with c = 1/2OOOn.
l
l
l
l
l l
l
54
GREG REESE
processed images than in the original. Note that discrete IDS simultaneously enhances the input in dark areas (the upper left and right) and in white areas (the spinal column). The images in Fig. 9 through Fig. 12 demonstrate that as the scaling parameter ‘c decreases, details in high-intensity regions become more visible, but at the same time edges become wider. Thus the user makes a particular choiCe of c based on a trade-off between having details in high-intensity areas or having crisp edges. In all cases, c is the only parameter in discrete IDS, so there is no need to worry about trying to relate and fine-tune a plethora of parameters. Figure 13 is a satellite photograph of an area in the arctic region. The jagged, vertical lines descending from the top of the picture are ice leads, fractures or passageways through sea ice. Unfortunately, the leads are often covered by clouds or haze and are difficult to see in unenhanced images. Figure 14 is the discrete IDS output of Fig. 13. In the processed image: The large leads are more visible and they extend further down the picture. The small leads branching off of the large leads are now clearer. . There is a long horizontal line in the middle of the discrete IDS image that is barely visible in the original picture. This line is a horizontal continuation of one of the vertical leads on the left side of the picture. Details in the lower left stand out more. l
FIG. 14. IDS output to ice leads in Fig. 13.
l
l
Discrete IDS does a good job of looking through hazy clouds and seeing the underlying scene, in this case the ice leads. It works well because the haze acts like a fairly uniform shadow and as Fig. 6 and Fig. 7 demonstrate, discrete IDS is fairly insensitive to uniform changes in illumination. However, if the cloud cover is too thick to see through, no amount of enhancement will bring out the scene beneath the clouds. In this situation, discrete IDS will enhance the clouds themselves (see the middle and top right of Fig. 13) and such improvements may be useful to meteorologists. 5.
CONCLUSION
IDS is a theory that is a good model of how people see patterns of light and dark, but it is also a good tool for improving digital images. The version of IDS for the computer, discrete IDS, enhances the edges, nonlinearly compresses the dynamic range, and automatically adjusts itself to the local intensities of an image, all without the use of windows or statistics. IDS is particularly useful for bringing out detail in regions of images that have low contrast. If discrete IDS is calculated with a look-up table, the only computations involved are additions and references to the table. Discrete IDS would be very well suited for operation on a parallel computer. ACKNOWLEDGMENT FIG. 13. Ice leads in the Lincoln Sea.
My thanks to Rich Goldkamp for his help in acquiring images and patiently shuttling them among various computer systems.
IMAGE ENHANCEMENT
REFERENCES 1. R. Alter-Gartenberg, F. 0. Huck, and R. Narayanswamy, Image recovery from edge primitives, J Opt. Sot. Am. Ser. A 7(5), May 1990, 898-911. 2. R. Alter-Gartenberg and R. Narayanswamy, Image coding by edge primitives, in Proceedings, Visual Communications and Image Processing ‘90, SPIE, Vol. 1360, pp. 1608-1619, Laussanne, Switzerland, 1990. 3. T. N. Cornsweet, Prentice award lecture: A simple retinal mechanism that has complex and profound effects on perception, Am. .J. Opt. Physiol.
Opt. 62(7), 1985, 427438.
4. T. N. Cornsweet and J. I. Yellott, Jr., Intensity-dependent spatial summation, J. Opt. Sot. Am. Ser. A 2(10), Oct. 1985, 1769-1786.
BY IDS
55
5. R. Narayanswamy, R. Alter-Gartenberg, and F. 0. Huck. Robust image coding with a model of adaptive retinal processing, in Optics, Illumination and Image Sensing for Machine Vision V, SPIE, Vol. 1385, pp. 83-103, Boston, MA, 1991. 6. AC-200 Series Adaptive Imager Products, Odetics Artificial Intelligence Center, 1515 S. Manchester Ave., Anaheim, CA 928022907.
7. A. Rosenfeld and A. C. Kak, Digital Picture Processing, Academic Press, Orlando, 1982. 8. V. T. Tom, Adaptive filter techniques for digital image enhancement, in Proceedings, SPIE, Digital Image Processing, 198.5, Vol. 528, pp. 29-42.