Chapter 2. Histograms and Point Processes

2.1 Image Histograms

Distribution of Gray Levels. We model an image as an array of digital values {f(m,n): 0 m M-1, 0 n N-1}, where for any row m and column n, the pixel position (m,n) has pixel value f(m,n). Figure 2.1 presents the array in the dimensions x and y. There are L integer values that f(m,n) can assume. These L values represent gray levels from 0 (black) to L-1 (white, or highest intensity of light). L is most often 256, and this many shades of gray are universally displayable nowadays on computer graphics systems.

The pixel value f(m,n) can be considered to be the height in a third dimension above the point (m,n) in the xy-plane. Figure 2.2 shows a 3-dimensional graph of an image of the "+" symbol according to our model. The background is black, or 0 level, while the symbol pixel values are L-1 = 255 (white). In this simple case, the proportions of pixels at the gray levels is zero except at gray level 0 and gray level L-1. Let P be the number of pixels at gray level 0 and Q be the number of pixels at gray level L-1, so that P+Q is the total number MN of pixels. Then the distribution of proportions p_k of pixels at the various gray levels g_k, 0 k 255 is

g₀ = 0 (black): p₀ = P/(MN)

g_k = k: p_k = 0 (1 k 254)

g_L-1 = 255 (white): p₂₅₅ = Q/(MN)

Figure 2.1. An MxN Image Figure 2.2. A Simple Image of "+"

Figure 2.3 shows a graph of the proportional distribution over the gray levels 0 to L-1, where we have taken P = 160 and Q = 96. This means that P+Q = 256 (this small image is 16x16 pixels). In this case, the proportion of black pixels is 160/256 = 0.625 and the proportion of white pixels is 96/256 = 0.375.

Histograms. Now consider the image "lena256.pgm" from Chapter 1. The gray levels range from 0 to 255 and the image size is MxN = 256x256 = 65,535 pixels. Each of the gray levels from 0 to 255 can be represented by a byte (8 bits), so the image file contains a short header and the 65,535 pixel bytes. Let us take the counts of the number of pixels at each of the gray levels g₀ = 0,...,g₂₅₅ = 255 and denote these counts by c_k, for k = 0,...,255. To obtain the proportions of pixels at each gray level, we divide each count c_k by the total number of pixels to obtain h_k = c_k/65,535, k = 0,...,255. Rather than look at these 256 proportions, we display them in a graph form, as in Figure 2.3. Proportions are the probabilities that any randomly drawn pixel from an image will be a certain gray level. The image of Figure 2.2 has a graph of proportions that shows the Pr[black pixel] = 0.625 and Pr[white pixel] = 0.375.

In general, an MxN image has MN pixels. The counts of the number of pixels at each gray level k, c_k, k = 0,...,L-1 are divided by the total number of pixels to obtain the proportions h_k.

h_k = c_k/(MN), k = 0,...,L-1 ..................................................(2.1)

_(k=1,MN) h_k = (1/(MN)_(k=1,MN) c_k =

MN/MN = 1 ........................................................................(2.2)

Given an image {f(m,n): 0 m M-1, 0 n N-1} we call the graph of points {(k,h_k): k= 0,...,L-1} the histogram of the image.

A histogram is a summary that we can view to see some basic characteristics of the image. Sometimes it is convenient to average every two consecutive proportions and display them, so that there are 128 (L/2) rather than 256 (L) proportions. A graph of averages of r consecutive proportions at a time is called an approximate histogram.

If the nonzero proportions are over a small band of gray levels, then the gray levels are few and are close together so that there is little contrast and it will be difficult or impossible to see all features in the image. If the histogram is spread out across all gray levels with approximately equal heights (a rather uniform distribution) then the image will not only have good contrast but will represent all gray levels approximately the same. If the histogram has high spikes at a few gray levels and low values for most values, or if there are essentially no pixels at certain gray levels, then there are details hiding in these plentiful grays and we should spread out the pixel distribution more uniformly.

Figure 2.4 shows the approximate histograms for four images. Part (a) is too dark and has low contrast due to a narrow range of gray levels, while the second, Part (b), has good contrast because of a wide range of shades of gray from very dark to very light. Part (c) is low contrast and is too bright. Part (d) is too low in contrast and is devoid of blacks and whites. The image of Figure 1.1 of Chapter 1 has low contrast with pixels at midgray levels, and so would have the type of histogram as shown in Figure 2.4d. Figure 1.2 shows an image that has an equalized histogram similar to that of Figure 2.4b.

If all pixels of an image were at a single gray level, then the image would be uniformly gray at that shade of gray. It would contain no information. If the pixel values were distributed randomly over all shades of gray, that is, from 0 to 255, then there would also be no information. An arrangement of gray levels that captures a scene of objects contains information about that scene. To the extent that there are too few grays, or there is random error on too many pixels, called noise, that information is degraded. There are also other kinds of degradation such as motion of the camera or objects (see Chapter 11), or an unfocused lens.

Figure 2.4. Four Approximate Histograms

(a) Low Contrast Dark Image ......................................................................... (b) High Contrast Equalized Image

(c) Low Contrast Bright Image ........................................................... (d) Low Contrast Image with Middle Grays

A Histogram Algorithm. We present an efficient algorithm here for computing the histogram of an image. The algorithm is described in a high level pseudo language that is easily translated into C, C++ or Pascal. The histogram can be graphed and displayed on UNIX type systems using xgraph, but we can also use XV.

Algorithm 2.1: Computation of Histogram

for k = 0 to L-1 do /Initialize all counts/

.........c[k] 0; /c[k] = count of pixels at gray level k/

for m = 0 to M-1 do /For each row and for each/

.........for n = 0 to N-1 do /column in the image/

................c[f[m,n]] c[f[m,n]] + 1; /increment count at gray level f(m,n)/

for k = 0 to L-1 do /Proportionalize each gray level count/

.........h[k] c[k]/(M*N); /M*N = total pixel count, h[k] is proportion/

/Note: f(m,n) is an integer in [0,L-1]/

2.2 Grayscale Transformations

Each image has an average gray level µ and a variance ² computed from

µ = _(k=0,L-1) kh _k ................................................................................(2.3)

² = _(k=0,L-1) (k-µ)²h _k .........................................................................(2.4)

Without apriori information we would expect that µ is approximately L/2. However, many important details may be in the lower range (or in the upper range) of gray levels. In that case we want more pixels to be distributed over the gray level range of interest. It is sometimes useful to use the average µ over the entire image, or over a portion of interest, in a way that changes µ to a desired value. The gray levels can also be spread out to increase contrast.

It is possible to change the distribution of gray levels so as to dilate or contract desired intervals of gray levels. A linear transformation maps the interval of input gray levels to an interval of output gray levels in a proportional manner, that is, any gray level at a proportion p of the way through the input interval is mapped to a gray level at that same proportion p of the way through the output interval. A nonlinear transformation maps gray levels in the input interval disproportionately into the output interval. In the next sections we will see that transformations of pixel values changes the histogram. We use such transformations to change a histogram to a more desirable one for some purpose such as to expose details in a particular range of gray levels.

2.3 Linear Transformations of Image Grayscales

A linear transformation of an image operates on each pixel value, which is a gray level at a particular position in the image, to map it into another gray level at the same position. The input (argument) is a gray level f and the output is a new gray level g defined at position (m,n) by

g = af + b, or g(m,n) = L[f(m,n)] = af(m,n) + b .........................................(2.5)

The coefficient a is the slope (or gain) and b is the y-intercept (or bias).

Figure 2.5 presents a linear transformation that maps the gray levels of an input image into the gray levels of an output image. In this case the transformation dilates the input domain from a subinterval [f_min,f_max] onto the full interval [g_min,g_max] = [0,255].

A linear transformation of the input gray level interval [f_min,f_max] onto the output interval [g_min,g_max] has the form of Equation (2-5) above, where the slope is

a = (g_max- g_min)/(f_max- f_min) ...............................................................................(2.6)

and so

g = af + b = [(g_max- g_min)/(f_max- f_min)]f + b ......................................................(2.7)

When f = f_min we desire that g = g_min, so we can substitute the point (f_min,g_min) in Equation (2.7) to solve for b via

b = g_min - [(g_max- g_min)/(f_max- f_min)]f_min .............................................................(2.8)

Substituting for b in Equation (2.7) and collecting terms, we obtain the results

g = [(g_max- g_min)/(f_max- f_min)](f - f_min) + g_min ....................................................(2.9a)

g(m,n) = [(g_max- g_min)/(f_max- f_min)](f(m,n) - f_min) + g_min ...................................(2.9b)

Figure 2.5. Dilating Linear Transformation ..........................................Figure 2.6. Piecewise Linear Transformation

Equation (2.9) maps f_min into g_min, f_max into g_max, and everything a proportion p of the way between f_min and f_max into that proportion p of the way between g_minand g_max. The smaller [f_min,f_max] is mapped into the larger [0,255] in a linear fashion in Figure 2.5. This transformation may be used to map a smaller interval of gray values into a larger one, but also may be used to map a larger interval into a smaller one. Another use is to break the range of gray levels into subintervals and use a linear transformation on each subinterval. Figure 2.6 shows such a piecewise linear transformation. The only requirement is that the conditions 0 f_min f_max L-1 and 0 g_min g_max L-1 are met.

The algorithm given below performs linear transformations, where f_min, f_max, g_min and g_max are selected beforehand for the appropriate mapping. The input image is {f(m,n): 0 m M, 0 n N} and the output image is {g(m,n): 0 m M, 0 n N}.

Algorithm 2.2: Linear Transformation

a (g _max - g_min)/(f_max - f_min); /Compute slope for linear transformation/

for m = 0 to M-1 do

.......for n = 0 to N-1 do

......................g[m,n] a*(f[m,n] - f _min) + g_min; /Transform each pixel f[m,n] into g[m,n]/

2.4 Nonlinear Transformations

Logarithmic Transformations. A nonlinear transformation is usually done after a linear transformation has set the contrast and range of gray levels to that desired. It maps small equal intervals into nonequal intervals. Suppose that most of the pixels have values at the lower end of the gray scale and we want to spread them out to see the detail there, but that we don't care about the lighter (brighter) values in the upper range of grays. Then we want a small input interval at the low end to map to a larger interval at the low end for the output image. We also want 0 to map to 0 and 255 (or L-1) to map to 255 (or L-1). The function

g(m,n) = (c)log₂(f(m,n) + 1) ................................................................(2.10a)

spreads out the lower gray levels. It must map 0 to 0, and 0 = (c)log₂(1) = 0. It also must map 255 to 255, so that 255 = (c)log₂(255 + 1) = (c)log₂(256) = 8c. Thus c = 255/8 = 31.875, so we have

g(m,n) = (31.875)log₂(f(m,n) + 1) ......................................................(2.10b)

Figure 2.7 shows this type of mapping. For example, 128 maps to 31.875log₂(128 + 1) (31.875)(7.001) = 223.157, which is truncated to the integer 223. Thus the middle gray level is strongly dilated. We can also use

g(m,n) = (c)log_b(af(m,n) + 1) ..............................................................(2.10c)

where a > 0 is a constant and b > 1 is a logarithm base. X-ray images are known to satisfy the intensity function f(m,n) = f₀exp[-r(m,n)], where r(m,n) is the attenutation of the x-ray signal at (m,n) due to the density and thickness of the material. Therefore, we use logarithmic transformations (to the base e) to enhance the detail on x-ray images.

Exponential Transformations. Another case is where we are interested in spreading out the upper gray levels at the expense of the lower gray levels, which must be contracted. Again, we want the end points to map to end points. While a logarithmic function spreads out lower levels disproportionately, an exponential spreads out the upper levels disproportionately. If we use

g(m,n) = exp(af(m,n)) - 1 ..................................................................(2.11a)

then 0 maps into exp(0) - 1 = 0. To force 255 to map into 255, we have that 255 = exp(a255) - 1, so that 256 = exp(255a). Upon taking the natural logarithm of each side, we obtain ln(256) = 255a, or a = ln(256)/255 = 0.0217458. Then this mapping is

g(m,n) = exp(0.021746f(m,n)) - 1 ....................................................(2.11b)

Figure 2.8 shows an exponential type of transformation. As an example, 128 maps to exp(0.0217458(128)) - 1 = exp(2.78346) - 1 = 16.1749 - 1 = 15.1749, which is truncated to 15. Thus the middle gray level is strongly contracted. If we desire to stretch the middel gray levels, we can also use a piecewise linear transformation as shown in Figure 2.6, or fit a nonlinear map on the middle subinterval. We can use, where a, c > 0 are constants and b > 1 is the base, the map

g(m,n) = cb^af(m,n) - 1 .........................................................................(2.11c)

Figure 2.7. Logarithmic Mapping ..........................................................................Figure 2.8. Exponential Mapping

2.5 Thresholding

One or more threshold gray levels can be used to transform an image pointwise. Let f₁, f₂ and f₃ be thresholds. Suppose, for example, that we want to map all gray levels below f₁ into black (0), those between f₁ and f₂ into a larger range of lighter levels and everything betweem f₂ and f₃ into black. Let us map everything above f₃ into white. Figure 2.9 shows a function that performs this transformation. An unlimited variety of maps can be constructed in a similar manner. Note that the gray levels between f₁ and f₂ are dilated, while all other data are severely compressed. The use of thresholds is a powerful tool for exposing certain details in an image.

Figure 2.9. Using Threshold Mappings

2.6 Histogram Equalization

Cumulative Distributions. A histogram {h_k} for an image may have its nonzero proportions predominately in the lower, upper or middle part of the grayscale (see Figures 2.4a, 2.4c and 2.4d). Ideally, the image grays should cover the range [0,L-1] and not have too many or too few counts in any gray levels. A transformation that spreads out the gray levels used and also changes the proportions to be more uniform is called histogram equalization. Figure 2.4b shows an approximately equalized histogram.

Consider the distribution of gray levels in Figure 2.10 (shown as a continuous function for convenience rather than as a discrete one). The area under the curve is unity as it represents the total proportions over all gray levels designated here by f. This is the same as a probability density function (pdf) h_F(f) for the random variable F that assumes grayscale values f. The total probability is

L

h_F(f)df = 1............................................................................. (2.12)

0

For the purposes of this discussion, f is a continuous grayscale variable with normalized domain of [0,1].

The accumulated probability at any grayscale value f is

f

h_F(r)dr = H_F(f) = g ..................................................................(2.13)

0

The function H_F(f) that accumulates area up to each point f is called the cumulative distribution function (cdf) for the pdf h_F(f). Figure 2.11 shows the cdf H_F(f) for the pdf of Figure 2.10. Cdf's are strictly monotonic and thus one-to-one and have values between 0 and 1. We now use this function H_F(f) as a nonlinear transformation function on the gray levels. The objective of such a transformation is to obtain transformed gray levels g = H_F(f) that are uniformly distributed across the grayscale, that is, the pdf h_G(g) is a constant over all g.

We now show that the transformation g = H_F(f) maps the gray levels f to gray levels g such that h_G(g) is a uniformly distributed pdf for the random variable G that assumes the gray levels g. The subscripts on the pdf indicate what random variable it describes, F = f or G = g. From probability theory for the transformation of random variables F G, we have

g = H_F(f).................................................................................... (2.14)

h_G(g) = h_F(f)[df/dg].................................................................... (2.15)

The inverse transformation is

H_F^-1(g) = H_F^-1(H_F(f)) = f ..............................................................(2.16)

which exists as an inverse function because H_F(f) is strictly increasing and one-to-one and onto [0,1].

Figure 2.10. A Density Function Figure 2.11. A Cumulative Distribution

Upon differentiating Equation (2.13), we obtain

dg/df = dH_F(f)/df = h_F(f) ....................................................................(2.17)

By the one-to-one property and Equation (2.17)

df/dg = 1/(dg/df) = 1/h_F(f) .................................................................(2.18)

Upon substituting Equation (2.18) into Equation (2.14), we obtain

h_G(g) = h_F(f)[1/h_F(f)] = 1, 0 f 1 ...................................................(2.19)

so that the pdf h_G(g) (or histogram) of g is the constant 1 on the normalized gray levels [0,1]. While

this is true in the continuous case, it is only approximate in the discrete case. However, it does tend to spread out the gray levels and approximately equalize the proportions at the various gray levels. We have shown that g = H_F(f) is a grayscale transformation f g on [0,1] to [0,1] such that g has a uniform distribution. Because g satisfies 0 g 1, we need to multiply it by L-1 to obtain the gray levels, so that 0 (L-1)g L-1.

The Equalizing Transformation. Let {f_k: k = 1,...,L-1} be the set of gray levels with histogram proportions {h_k: k = 0,...,L-1}. Then the gray level transformation is given by

g_k = _(j=0,k) h_j ......................................................................................(2.20)

where the summation replaces the integral of Equation (2.13). An algorithm for histogram equalization is given below.

Algorithm 2.3: Histogram Equalization

sum 0.0; /Initialize sum to zero/

for k = 0 to L-1 do /For each gray level/

........sum sum + h[k]; /sum histogram proportions/

H[k] sum; /Collect cumulative values/

for m = 0 to M-1 do /For each row and each column pixel f/

........for n = 0 to N-1 do /position, compute new gray level g/

.................gray_level f[m,n]; /Get old gray level f as index/

................g[m,n] (L-1)*H[gray_level]; /Compute g, 0 g 1 and multiply the/cumulative proportion by max. gray level/

Notice that in the last line we could have put g[m,n] (L-1)*H[f[m,n]], but for clarity we have used an intermediate indexing variable gray_level. The transformed gray levels g[m,n] are approximately uniformly distributed across the grayscale.

2.7 Statistical Techniques

Statistics-based Linear Transformations. Let µ be the average pixel gray level over an image {f(m,n)} and let ² be the variance. If the overall image is too dark or too light, we may choose a desired mean µ_d. If the contrast is too low or too high, we may choose a desired variance _d², where _d > or _d < , respectively. We compute the slope (gain) for a linear transformation g = af + b via

a = _d/ .............................................................................................(2.21)

We can now compute the y-intercept, or bias, b per

b = µ_d - aµ .........................................................................................(2.22)

Figure 2.12. 5x5 Neighborhoods

The following linear transformation is called the statistics-based linear map and can easily be worked into Algorithm 2.2.

g(m,n) = [_d/]f(m,n) + (µ_d - [_d/]µ) = af + b ........................(2.23)

Statistical Differencing. For the next method, we need the concept of a neighborhood of a pixel. For our purposes here, a pxq neighborhood is a rectangular array of pixels of p rows and q columns, where p and q are odd integers. Thus we chould write it as a (2r+1)x(2s+1) neighborhood of the center pixel. Figure 2.12 shows a 5x5-neighborhood of a pixel at position (m,n).

Statistical differencing is a local adjustment method that tends to produce a similar contrast throughout the image and the degree of contrast is user selectable. It uses statistics over a large neighborhood, say 29x29 or 31x31, to adjust the center pixel. The statistics µ and are computed over the neighborhood pixels and applied to the center pixel value f(m_c,n_c) via

g(m_c,n_c) = µ + ß[f(m_c,n_c) - µ]............................................................ (2.24)

where ß = _d/ and _d is the desired parameter. We select _d > or _d < to respectively stretch (dilate) or compress (contract) the output grayscale.

The user control of this method can be improved by use of a desired mean parameter µ_d per

g(m_c,n_c) = µ_d + (1- )µ + ß[f(m_c,n_c) - µ].......................................... (2.25)

where 0 < < 1. This gives the user the ability to adjust the average levels µ on the neighborhoods up or down. Another adjustment is needed to prevent ß from being too large when is too close to zero over a neighborhood of similar values. We therefore use ß₀ in place of ß, where

ß₀ = r_d/(_d+r) ......................................................................................(2.26)

When = 0, ß₀ = r. By choosing r such that r > 1 or r < 1 (to move g(m_c,n_c) farther from, or closer to f(m_c,n_c)), we may dilate or contract the difference = f(m,n) - µ. The computation required is much greater than for the statistics-based global adjustment because the neighborhood statistics must be computed for every pixel in the image. When a pixel is near the boundary, we use only those neighborhood and mask entries that intersect the image. Algorithm 2.2 may be modified to use Equation (2.25) per

g(m_c,n_c) = µ_d + (1-)µ + [r_d/(_d+r..)][f(m_c,n_c) - µ] .............................(2.27)

Some rather dramatic enhancement effects can be achieved by the use of statistical differencing in the form of Equation (2.27).

2.8 Using XV for Histogram Transformations

Specifying Transformations with XV. The image processing program XV (xview) permits histogram equalization by left-clicking on a bar that says equalization, but it also has another powerful tool. It displays a graph of the identity transformation (the identity maps any pixel to itself, that is, g(m,n) = f(m,n)). However, there are four small square nodes on the line that is the transformation function. Two of the nodes are at the end points and the other two are in between. By left-clicking on any one of these nodes, holding down the button, and dragging the node up or down, the transformation function shape is changed to a curve that fits through the square nodes.

We can therefore specify nonlinear transformations of the gray levels by changing the location of the nodes, at which time the transformation curve instantly assumes a new shape through the nodes and the image on the screen changes immediately to reflect the new gray values of the pixels.

We nust first be logged into a UNIX based computer workstation and be in X-windows (see Section 1.5 in Chapter 1 for a description of how to do this). Then we load an image and display it via

/images$ xv shuttle.tif<ENTER>

When the image appears on the screen, right-click inside the image to obtain the xv controls window. From the center position along the top of the control window, bring the pointer down to the Window item, click and hold down the left mouse button and move down to the Color Edit item. Release the button. A new larger window now appears on the screen with xv color editor written in its top border. Move the image shuttle.tif up to the top of the monitor screen by left-clicking in its top border, holding the button down, and dragging the image upward and then releasing the button.

Now left-click inside the new xv color edit control window to bring it to the forefront. Move the pointer to the bottom central region, where the identity transformation graph resides with a diagonal line going from the lower left to the upper right. This is the histogram transformation function. Left-click on one of the square nodes, hold the button down, and move the node up or down. Upon releasing the button, a curve fits through the four square nodes. Left-click inside the image now to bring the image to the forefront. It changes to reflect the transformed gray levels of its pixels due to the displayed transformation. Click on the Reset button in the graph square to restore the original histogram.

Histogram Equalization with XV. To equalize the histogram of shuttle.tif, left-click the HistEq item on the lower left of the xv color edit window. The histogram transformation graph mentioned above and the image change immediately. To save the equalized image, left-click inside the xv controls window to bring it to the front of the screen and then left-click on the save button. Type in a new name at the bottom and select PGM, raw data format at the top of the pop-up menu and then left-click on OK. To quit XV, simply click on the Quit item in the xv controls window (at the bottom right). Figures 2.13 and 2.14 show the original and equalized versions of shuttle.tif.

Figure 2.13. Original "shuttle" ...................................................................................Figure 2.14. Equalized "shuttle"

2.9 Exercises

2.1 Find a linear transformation that maps [0,255] into [64,128]. What will this transformation do to the histogram? How will this affect the histograms of Figure 2.4?

2.2 Can you make a general statement of the effects any linear transformation will have on the histogram? How will the shape change? Will the proportions change? Explain.

2.3 Suppose that an image had a histogram that showed a high proportion of pixels in the middle gray range, and that here is where it is suspected that a lot of detail lies. Design a transformation that blacks out the lowest and highest gray levels and stretches the middle grays.

2.4 Use xview to display lena256.tif and then use the histogram transformation function by moving the tiny squares to map the gray levels in a similar fashion to the transformation shown in Figure 2.6 (the nonlinear transformation provided by xview will approximate the piecewise linear transformation of Figure 2.6). What is the result? What does this transformation do to the image to make it look different?

2.5 Use xview to convert both lena256.tif and shuttle.tif to the raw data files lena256.pgm and shuttle.pgm, respectively (see Chapter 1 for instructions on converting file formats). Now display lena256.pgm and shuttle.pgm on the screen using xview. What is the difference between, say, lena256.tif and lena256.pgm in the display?

2.6 Design a transformation of gray levels that stretches out the darkest and lightest shades of gray but compresses the gray levels between, say, 80 and 180.

2.7 Implement the design of Exercise 2.6 above with a nonlinear approximation using xview and the histogram transformation (graph) function in the color edit window.

2.8 Graph the linear transformation with slope of -1 that maps 0 to 255 and 255 to 0. What does this transformation do?

2.9 Use XV to display "lena256" and use the windows edit control window to adjust the transformation of the histogram to implement the transformation given in Exercise 2.8 above. Next, display the image. Is this the negative of the original?