- Email: [email protected]
Image enhancement using sliding histograms
COMPUTER VISION, GRAPHICS, AND IMAGE PROCESSING 44, 2 1 1 - 2 2 9 (1988)
NOTE Image Enhancement Using Sliding Histograms PAVEL A. CHOCHIA
Institute of Information Transmission Problems, The USSR Academy of Sciences, Moscow, USSR Received April 26, 1986; accepted January 7, 1988
Image enhancement methods based on local histogram analysis are discussed. It is shown that this problem is closely connected to that of contour preserving image smoothing and may be solved by a preliminary image decomposition. A method for image decomposition into background and texture components is proposed. The results obtained are used for image smoothing, enhancement, object and contour detection, and color image enhancement. © ]988 Academic Press, Inc.
1. INTRODUCTION
In this paper we consider image enhancement by local contrast amplification of both fine details and texture [12]. We show that this task is closely related to the image smoothing problem and may be solved by a preliminary image decomposition into two different information components. The aim of this paper is to develop an image smoothing method using the probability distribution function for a sliding fragment, and then to use it in some image processing tasks. This problem is based on robust estimation of sliding histogram parameters using order statistics [2-4]; the median filter, suggested by Tukey [2], is the most useful. Statistical properties and contour preserving facilities of the median are well known [5-7], and fast algorithms have been developed using sliding histograms to accelerate its computation [8]. For image smoothing and enhancement, a model of local image fragments is proposed. Using this model, we will show that both local averaging and local median methods achieve unsatisfactory image smoothing results. The new proposed image smoothing method is based on the M-type estimators used in order statistics [4]. This method makes it possible to decompose the image into two components. This operation leads us to a satisfactory solution of the image enhancement problem and some other problems. A digital monochromatic image is an M x N matrix X = [xm, ], and the image value x,~, represents the optical density (logarithm of the luminosity), quantized into K levels, i.e., 0 < x,,, < K - 1. In the following two different sets of elements in the vicinity of x will be used (Fig. 1). To differentiate them we use the terms "neighbourhood," V,,, (for inner) and "fragment," Wm, (for outer). Both of them are just moving rectangular windows, but they have different sizes. In real experiments they were squares from 3 × 3 up to 9 x 9 pixels (neighbourhood) and from 11 X 11 up to 65 X 65 pixels (fragment). The fragment histogram, Hmn = {hmn(k)), 211 0734-189X/88 $3.00 Copyright © 1988 by Academic Press, Inc. All rights of reproduction in any form reserved.
212
PAVEL A. CHOCHIA Wren
FIG. 1.
An image fragmentcontaining three domains (U1, U2, U3).
is a probability distribution for the magnitudes of the elements: E h , n n ( k ) = 1,
k=(0 ..... K-
1).
k
We assume that a discrete random quantity 4, 0 < ~ < K, has a normal distribution N(# a aS), if its distribution function F ( k ) coincides at all points k, where ~ can be nonzero, with the distribution function of a continuous normal random quantity ~/, having the same mean value and dispersion. 2. MODEL OF AN IMAGE FRAGMENT It is assumed that the original image x,,, is the sum of two statistically independent components [13, 14]: Xmn = Srn n -~ tmn ;
(1)
a background component Stun, which determines the mean intensities of large size image details, and a texture component, tmn , which contains information about the texture and small details of the image. In order to separate Stun, one has to know the difference between statistical properties of the background and texture components. Therefore before starting the analysis and the development of the transformation methods, we must choose a model of the image, enabling us to take into account the statistical correlations of the image elements within a fragment. We assume that the image consists of connected domains U 1, U 2. . . . . u N , . . . , separated by contour boundaries. According to the model of a neighbourhood [15], an image is a sum of two nonoverlapping sets: contour elements and internal elements of the domains. For elements of the neighbourhood x m , ~ Vm, (r = 1 . . . . . R), belonging to one and the same set, one can draw a plane (say, by means of the method of least squares) so that the relation holds Xmn
= ~mn
-]- ern(armnP r + }~mn) + (1 -- ern)(fl,~np" + ~/rn) ,
(2)
where e,n nr is a contour mask (e = 1 for the contour element and e = 0 for internal one), /X,nn is a constant determining the mean value of the neighbourhood elements, pr is the distance from the center (m, n) of the neighbourhood to the r t h element,
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
213
a r . and fl~, are random quantities which determine a slope of the plane in the direction of the rth point: a~,, (for contour) and B~, (for internal) are equal to the cosine of the angle between the vector from (m, n) to the rth point, and the plane gradient vector G,, n. ~ , and 7/~,, are independent random quantities having normal distributions N(0, o{) and N(0, ~ ) , respectively. The central element is Xrnn = ~mn "{- emn~mn +
(1
-
emn)~mn.
The one-component version of this model is close to model [16]. This model describes statistical correlations between elements within a single neighbourhood, limited to 3-5 discretization steps. However, if a contour line does not separate the points with coordinates (m, n) and ( m + e, n + 3), then there are statistical dependences between the values of /~m, and /~m+~,,+n for sufficiently large e, 3 (up to 20-40 and larger), since these points belong to the same image domain [17]. This fact can be used in constructing the model for a sufficiently large fragment. Let Wren be a fragment with its center at (m, n), which is divided by contour boundaries into J background domains U1,..., U J. Points of a domain U j belong to the same extended detail of the fragment and form a connected region x~ ~ U j c Wren (i = 1 , . . . , I J) (Fig. 1). Certain quantifies #i correspond to the elements x i. We shall assume that within a single domain U J, at distances not exceeding the fragment size, values of /~i are described with sufficient accuracy in terms of a piece-wise constant model ' + ~ i = S rnn
(3)
where s~, is the mean intensity of the domain U j inside the fragment Wrnn,and x/ are (perhaps, dependent) random quantities with a normal distribution N(0, o~). Such a model is quite reasonable for most real images. Hence an element x / o f the domain u J is given [15] by
X! = SJn ~- I£i q- ~i ~---SJmn+ T/.
(4)
The physical meaning of two added terms is quite different. The first component Ki determines the texture and small details of the image, the values of its elements are substantially correlated. The second component 7/i represents non-correlated Gaussian noise. As Ki and ~i are mutually independent random quantities with normal distributions N(0, a~) and N(0, o~), respectively, their sum ~'i = xi + ~1i is also a random quantity with the normal distribution N(0, 02), where o~ --- of + a~. The central element of the fragment in the domain U 1 is 1 1 1 _if_ I£mn 1 _]_ ~ 1 n = Sm n + ,r;nn. Xmn ~- Stun
(5)
It follows from (4) that if a fragment is in a single domain of the image, its histogram is unimodal. If the fragment covers parts of two or more ( J ) domains, the histogram is J-modal. Positions of the modes are determined by mean intensities of the domains ( S 1, S 2. . . . , SS), and powers of the modes are proportional to fractions of the area of the fragment occupied by corresponding domains (Fig. 2). Filling the
214
PAVEL A. CHOCHIA h(k)
k
S1
s3
S2
FIG. 2. The histogram of the image fragment (Fig. 1).
gaps between the modes depends on the mode dispersions (OrJ ) 2 and on intermodal distances, which determine overlapping of the distribution tails, and on a fraction of the contour elements in the whole distribution. Thus the problem of the evaluation of Sm. for image point (m, n) can be solved as an estimation of a mean intensity level s.,nl in the central domain U 1 of the fragment Wm.. This problem will be discussed in Section 5. We now review the known image enhancement methods. 3. IMAGE PROCESSING METHODS REVIEW
Processing methods providing an enhancement of local contrasts in an image can be characterized by two properties: linear/nonlinear methods and fixed/adaptive methods. The first characteristic means the linearity in the contrast amplification, and the second one is concerned with the varying of this amplification factor within the image field. All the methods are based upon the use of the fragment histogram H,,,. Assume that the fragment is a rectangle of size L × L, L = 2l + 1. The estimation employed is either the mean value of fragment elements K-1
A(Wm. ) = Y'. k h . , . ( k ) = L k~O
-2
m+l
n+l
Z
E
xij,
(6)
i=m--I j=n-I
or order statistics of the type R(q), 0 _< q < 1, defined by y-1
Rm.(q) = y,
y
if ~_, hm.(k) < q < ~_, hm.(k ). k=O
(7)
k=0
In some cases the median is used, namely m e d ( W ~ ) = Rm~(0.5 ).
(8)
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
215
1. Linear Fixed ( L F ) Methods The image is decomposed into two components: Xm. = Stun + tin.; various transformations are given by the formula : Ym. =atmn + bsm. + C.
(9)
The parameters a, b, c are constant throughout the image; in most cases c = (1 - b ) / 2 . Two methods for the evaluation of Sm, are used: the first is the average (6) over the sliding fragment [11, 12], and the second is the median (8) for the same fragment [12]. If a = b, then we obtain the standard linear correction formula Ymn = aXmn + C.
(10)
2. Linear Adaptive ( LA) Methods The formula used for the processing is similar to that of the LF methods (9), but the parameter a, which is the amplification factor of local contrasts, is determined from an estimate of variation of the fragment histogram Hmn [12]: Ymn = amntmn + bsrnn At- ¢"
(11)
Two variants are used to evaluate the parameter amn: (i) If Sin. is the average, then a,.. is determined by the dispersion of the element intensity distribution over the fragment am. = O r J ( % r + 8).
(12)
Here O2r is the dispersion of the intensity distribution for the original image in the fragment Wm.; ozs is the desired resultant dispersion; and 3 is a parameter restricting the contrast amplification. (ii) If Sm~ is the median, then amn is determined by the interquartile distance AR = R(0.75) - R(0.25), a,n . = ARres/(ARor + 3),
(13)
where R(0.25) and R(0.75) are magnitudes of the first and the third quartiles, respectively, and 3 is a parameter, as in (12). 3. Nonlinear Adaptive (?CA) Methods A known method of this type is the sliding histogram equalization method [9]. It is represented in terms of the order statistics Ym. = R m. (Zm.),
where z,.. = x . , . / K .
(14)
It is not difficult to extend this method, obtaining a sliding histogram modification. Let G ( k ) be the resulting probability distribution, and G - l ( z ) be the inverse
216
PAVEL A. CHOCHIA
function, i.e., G-I[G(k)] = k. Then the sliding histogram modification is given by
Ymn= G-t[Rmn(Zmn)].
(15)
The function G(k) can be chosen on the basis of certain requirements of the processing result. Variants employing a similar modification are known for the histogram of the whole image [10].
4. Nonlinear Fixed (NF) Method Methods of this type have not been found in the literature and are mentioned for completeness. They can be constructed, for example, with the formula
Ym, = f(tm,) + bsm, + c. Here f ( x ) can be an arbitrary nonlinear function. It is important to note that most methods considered use a preliminary separation of the image into two components, stun and tmn, for the following amplification of the component tin,. In the present work we restrict ourselves by studying the LF methods. As was mentioned above, the known LF methods differ in their determination of Sm~. Let us consider the accuracy in the determination of the mean intensity level Sm,l for the central domain U 1, containing the element Xm,, when the average or the median over the fragment W,,, is used. 4. AVERAGE AND MEDIAN OVER THE FRAGMENT
If the fragment Wm, is within a single domain of the image, then according to the model (4), the fragment histogram is normally distributed, and the expected values of the average M{A(W)} and the median M(med(W)} are equal to the mean intensity level S 1 for the image domain. If the fragment covers several image domains, the expected values M ( A ( W ) } and M{med(W)} will be equal to some intermediate quantities within the interval (Smi~, Smax). In this case the set of fragment elements can be regarded as a combination of the elements of J normal distributions, N( S 1, o2) .... , N( S J, o2). Let Pj be the probability that an arbitrary element of the fragment belongs to the domain UJ: EPj = 1. Then the expectation value of the average over the L × L square fragment is
M ( A ( W ) } = MIL-2~_,Xm. 1 = ~_,PiS/. \ W
"
(16)
j
The expected value of the median is Z0, which is the root of the equation
EP/F+(z) = 0 . 5 ;
(17)
J
Fj(z) is the distribution function for dements of the domain U/. Let us compare the results for J = 2. For simplicity we set o 1 = 02, AS = S 1 - S 2. The Pj dependences of M{A(W)} and M(med(W)} for AS = 20, 3o, and 40 are shown in Fig. 3a. It is seen, that for AS = 20 the average and the median have a
217
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS b A(W) med(W)
A(W) mecl(W)
4a
3a
3a 2a 2o
0
~ P2 0.5
I
0 ~
' 0.5
t 1
~- P2
FIG. 3. Dependences of the average (dashed line) and the median (solid fine) over a fragment for two-mode distribution: (a) a 1 = 02; I - AS = 2a, 2 - AS = 30, 3 - AS = 40; (b) AS = 30; 2 - o2 = 01,4--a 2=0.5al,5-o 2=0.
small difference. As AS is increased, M{med(W)} is almost invariably near the ends of the interval, where P2 < 0.3 or P2 > 0.7, but a step appears near P2 = 0.5, while the expected value of the average, M{ A(W)}, is always a linear function. The dependencies presented here have a clear physical meaning. Suppose that there are two domains in the image, separated by a contour difference, and intensities in each domain are described with the normal distributions N(S 1, 02) and N(S 2, 02). Then the graphs show variations of the average and the median for a fragment along a line perpendicular to the contour. Evidently, for % 4:02 the slope of the graph for the median is determined by the dispersion 02 for small P> and by o22 for large P2. The graph for the average over the fragment will not be changed (Fig. 3b). Results of processing for a real image (Fig. 4a) performed by means of (9), determining Stun from the average and the median, are shown in Figs. 4b, c. With a
FIG. 4. Results of the image processing by means of known methods: (a) a part of the original image (Fig. 10a); (b) processing by means of Eq. (9), Stun is the average over a sliding fragment; (c) the same, b u t stun is the median; (d) sliding histogram equalization by Eq. (14).
218
PAVEL A. CHOCHIA
FIt. 5. Test image processing (image dimensions 128 x 384 elements, fragment dimensions 33 × 33 elements): (a) the original image and the graph for one of its lines; (b) the average (6) and graph for the same line; (c) the median (8).
sufficiently high enhancement of the fine texture in flat background regions, a number of distortions have arisen during the processing: an over-contrast of the contour steps in the case of the average, and of the detail angles in the case of the median. In Fig. 4d we present for comparison the result obtained using the sliding histogram equalization algorithm (16); here is an over-contrast of the texture in flat image regions. Distortions are seen more dearly in the results of processing a test image (Fig. 5a). Magnitudes of the local average and the median are presented in Figs. 5b, c. It is seen that the transformations appreciably changed contours and shapes of the objects. The transformation using the average has smeared the contours uniformly, while the transformation using the median conserved extended linear contours but erased angles of the details. These distortions appear, because, in a fragment having the angle of the detail at its center, the total number of object dements is less than the number of background dements. Using such distorted estimates in (9) results in an over-contrast of contours, details, etc. Thus we can conclude that both the average and the median for a fragment lead to unsatisfactory estimates for the mean intensity level of the background in the central domain, Stun, 1 if the fragment covers two or more different image domains. 5. A METHOD FOR DETERMINATION OF BACKGROUND INTENSITY
Now we try to find a method suitable for determining the mean intensity level sin, , using a fragment histogram H,n,. A typical histogram for a fragment containing several domains UJ ( j = 1,..., J) is shown in Fig. 2. It consists of several modes, corresponding to domains in the fragment. An dement Xm, belongs to the central domain U 1 of the fragment, unless it is a contour element. The problem of evaluation of the background intensity, or equivalently the mean intensity level of the central domain S 1 is reduced to determination of the center of the correspond-
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
219
i n g distribution (mode) in the histogram. To solve this problem one has to: (i) choose the appropriate mode in Hm~ and (ii) determine the mean value of the distribution. By choice of any mode in the histogram, we determine an indication of the center of the corresponding distribution with some accuracy. It should be noted that by choosing the mode nearest to Xm. one does not always get correct results. Actually, in real images the summed dispersion o~ in (4) may be comparable to the difference of the mean intensities of adjacent domains IS g - SJl, and it is quite probable that the mode nearest to value x , ~ would be attributed to another domain. Analysis of real images shows that inside uniform domains o~ are less than % by a factor of 2 to 3; hence for choosing an appropriate mode one should use the value /*m~, Eq. (2), which is averaged over the elements of a neighbourhood Vm~, belonging to the same domain as xm~. For smoothing the image we may use a method, based on M-type estimators, used in order statistics [4, 20], which is a function of an ordered elements set ( ~ } = ~1 -~< ~2 ~
"'"
-~< i n ,
where w~ depends on the 4. If for any x i ~ Vmn we choose wi = 1 for Xmn < Xm~ + 8 and w~ = 0 for others x~, then we obtain the filter
Xmn =
WiXi
i=1
Wi.
li=l
(18)
This filter is used in [18] for image smoothing and is called a "sigma filter." Parameter ~ is recommended as 8 = 20, where o 2 is the noise dispersion. As was stated in [22], the sigma filter is one of the best edge-preserving noise cleaning methods. It follows from (5) that for a point (m, n) inside the domain U -i, Xmn has a distribution close to the normal one N(SJ, o~ + ~/Y.wi) , where S j is the mean intensity of the background component in UL To determine S 1 we shall construct an analogous sigma filter for a fragment. Fixing a value of A, we take points in the fragment W,,, (or equivalently, points of the histogram Hm,) which are within the window (Xmn -- A, Xmn "q- m), and calculate for them the average A(W; ~), or the median med(W; ~). Let us compare the accuracies of the estimates A(W; ~) and med(W; ~). Accordhag to the model (4), elements in one domain have a normal distribution N(S j, o~). Suppose that the values of the elements of neighbour domains are sufficiently remote and their influence is negligible. For simplicity we put a = 1 and use the notation I N T ( a , b) = (2~r) -1/2 ~bexp(-o2/2) do.
220
PAVEL A. CHOCHIA
/ .
med(W ;~)
A:
a
A=1.4o I
A=2a
0.5a
A=3a
,
.
0.5a
a
,
1.5a
.~ "~
2a
FIG. 6. The average A (W; ~) and the median med(W; Y) versus ~ (curves 1-4 and 5-8, respectively). The results have been obtained by means of the sigma filter (18) with A = a, 1.4o, 2a, 30.
The expected value for the average of elements within the window (~ - A, ~ + A) is
M{ A } = (2~r)-l/z ff+2v2exp(-vz/2) dv/INT(.~ =
(2~')-i/2[exp(- (.2×INT(£-
A)2/2)-
A, y + A)
exp(-(ff + A)2/2)]/
A, ff + A).
(19)
The expected value for the median, M{med}, for the same window is a root of the equation I N T ( ~ - A, M ( m e d } ) = INT(ff - A, y + A ) / 2 ,
(20)
In Fig. 6 we present the Y dependences of M(A} and M ( m e d ) for A = o, 1.4o, 20, 30. The ratios a = M{ A }/ff and a = M ( m e d } / Y determine the convergence velocity for the corresponding estimate. The comparison between the graphs shows that for all values of A, med(W; x) converge more rapidly than A(W; x). As one would expect, a increases with A from 1.4-1.5 at A = o to 3-5 at A = 20 and 15-30 at A = 30. It is noteworthy that if overlapping modes are present in the real distributions, the probability that an appreciable fraction of the second mode elements are inside the window (E - A, y + A) increases with A. Therefore it is unreasonable to take A > 20. An iterative procedure, where at the k t h step ~k+l is evaluated with the known xk, enables one to increase the initial convergence velocity a up to about a k, where k is the number of the iteration steps. For A = 20 and k = 3 we obtain a 3 = 100. In the determination of the mean intensity level in the central domain we have not taken into account its dimensions, i.e., the number of elements nl contained in U 1. It is evident, however, that with a decrease in n i a threshold value n e may be
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
221
attained, such that all the elements of U I must be, in fact, considered as a texture against the background of another domain UJ enclosing U 1. A choice of the threshold np is a nonformalized problem. It depends on the concrete subject, the image discretization step, and the final aim of processing. Therefore we propose to treat np as a parameter, proceeding in the following way. Let us suppose we have a fragment histogram Hmn, a parameter np, and an initial estimate for the mean intensity level of the background ~, calculated over a neighbourhood I'm,. The fragment histogram Hm~ can be used to determine the order statistics Rl(np/L 2) and R2(1 - np/L2). The basic value xl for the first iteration is given by the thresholding algorithm:
{
~,
X1 =
R1,
R2,
if R 1 < ~ < R 2, if X _< R1, if ~ > R 2.
(21)
This algorithm provides the desired result in the overwhelming majority of cases: if the fragment covers only two domains--the background (environment) U; and the central (small) domain U 1, as well as in situations where S 1 is either larger than the maximum or lower than the minimum intensity for the other domains of the fragment. Thus we derive the following algorithm for determination of the mean intensity stun for the central domain U 1 of fragment Wren. ALGORITHM 1 (A1). 1. Calculate histograms for the neighbourhood, Hm~,, centered at the point (m, n). 2. Evaluate ~ over the neighbourhood histogram
H°m~,and for the fragment, H~n, using the sigma filter
(18). 3. Evaluate X1 at a fixed rip, using algorithm (21). 4. Median (20) with the window (xk - A, ~k + A) is calculated by iteration from H,~,. The obtained quantity is considered as the desired s,~. Results of processing the test image by means of A1 at various np are shown in Fig. 7 (the test image is given in Fig. 5a). In the first case np was rather large, and
--1 F--L-
FIG. 7. Test image smoothing (original image is given in Fig. 5a) by means of A1. Neighbourhood dimensions 5 × 5 elements, fragment dimensions 33 × 33 elements: (a) n; = 80; (b) np = 10.
222
PAVEL A. CHOCHIA
FIG. 8. Distortion of contour step: (a) appearance of the double contour line when Stun is calculated by means of A1; (b) distortions disappear if modifiedvariant A2 is applied.
small details are filtered out of the image, as well as the noise. In the second case (at small np) the image is smoothed, but all details are retained. The accuracy of the reproduction of the original image (without noise), calculated with the r.m.s. criterion, was 0.3% after the first iteration. The graphs show a good smoothing of the noise, while the contour differences were conserved completely. 6. TEXTURE ENHANCEMENT Texture enhancement is performed by setting the difference tmn = Xm, -- S,~n in (9). Here Smn is the mean background intensity, obtained by means of A1. However, as experiments have shown, sometimes at domain boundaries a monotonous contour step is broken, and a false contour appears (Fig. 8a). One reason for these distortions is the fact that dements of the contour are ascribed by A1 to one of two neighbour domains of the image. Let us consider the origin of these distortions in more detail. Suppose we have two domains, U 1 and U 2, separated by a continuous contour step (Fig. 9a), their mean intensities are S 1 and S 2, and S 1 > S 2. The contour elements have intermediate intensities in the interval (S 2, $1). The result of the image processing by means of A1 will be as follows. Points to the left of the middleof the contour step will be ascribed to the region U 1 and after the processing their values will be close to S 1, while points to the right will have values close to S 2. The graph for this smoothed image (Sm,) is shown with dashed line (Fig. 9b). The same effect was discussed by Frieden [19]. Setting the values of Sin, and tin, into (9) with a = 3 we obtain the result shown in Fig. 9c. It is seen clearly that a " b u m p " appears in the contour, and, as a consequence, an additional line along the contour emerges in the transformed image. To avoid such distortions one should modify the procedure of determining Sm, for contour points. Evidently, if (m, n) is a point on the contour step, then the majority of fragment elements (in flat domains) will not be connected statistically with the d e m e n t Xm,. Therefore in the determination of Sm, for the contour point one should take into account only its nearest dements within the neighbourhood
223
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS ix(n)
2
s~ Ik l~.
Slt x(n) S2
FIG. 9. For the explanation of origin of the contour distortions: (1) the graph of the original contour step; (2) the value of Stunis calculated by means of A1; (3) the value of Ymn is calculated by means of Eq. (9).
Vmn. A sufficiently appropriate way to determine stun in this case is to take Xmn' calculated at the second step of A1. The proposed modification of A1 [21] is as follows. ALGORITHM 2 (A2). 1. Evaluate sin, by A1. 2. Identify contour points. To this end calculate an amplitude gm~ of the plane gradient G,~ used in model (2), by the formula gmn= [($m-l,n-lWSm-l,n+Sm-l,n+l--Sm+l,n-l--Sm+l,n--Sm+l,n+l)
"t"(Sm-1, n-1 "+ Sin, n-1 -1- Sm+l,n-1 --
Sm-l,n+l
--
Sin, n + l
--
2
Sm+l,n+I)2] 1/2
3. The value g,., is compared with the threshold & If gin, < 8, then Xm, is the dement inside a domain; then put g,~, = Sm,. If gin, > 8, then Xm, is the element of a contour overfall; then set g,,, = Xm," The obtained quantity g,,, is considered as the mean background intensity at the central point (m, n) 4. Evaluate the texture component, tm, = Xmn -- gmn" 5. Perform the image enhancement by transformation Ymn = atmn + bsmn + c.
(22)
Figure 10 shows the result of image enhancement by means of A2. It is seen that in the resulting image the contrast of texture and small details is enhanced considerably with respect to the background, while mean intensities of the extended details are not changed and significant distortions have not appeared. 7. COLOR IMAGE ENHANCEMENT
The method of local contrast enhancement presented above for black-and-white pictures can be extended by analogy to processing of color or multi-zonal images. In
224
PAVEL A. CHOCHIA
FIG. 10. The results of airphoto processing by means of A2: (a) the original image (512 × 512 elements); (b) local texture enhancement by means of A2.
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
225
aT s1
G
RE FIG. 11. Three-dimensional color histogram H(R, G, B) and multiplication of the original texture vector T up to a T.
this case any value of the image element must be considered as a point or a vector in a J-dimensional space of color coordinates (say, in the three-dimensional (R, G, B) space)
Xmn
Xmn ~- Stun "~ Kmn + ~mn = Stun + Tmn"
(23)
In this space a J-dimensional color histogram H ( k l , . . . , k j) for a fragment Wren is calculated at every point of which the probability of a given color value of an element is reflected. The points with non-zero probability condense into clusters (Fig. 11) corresponding to domains of a fragment (and to modes of the ordinary one-dimensional histogram). To define the background intensity Sin, of a color image (by analogy with the determination of background intensity of monochrome image by means of A1), it is necessary to separate a region v1 in the color histogram, corresponding to a central domain U 1 of the image fragment Win,,, and to find the location of its center (vector S 1, Fig. 11). This vector is considered as the value of the smoothed color image Stun at the point (m, n). Having determined the texture v e c t o r Trn n = X m n -- Stun , we obtain, in analogy to (22): Ymn = aTmn + bSmn --I- C.
(24)
Alternatively, conserving the original mean colors: Ym,, = aT,,,,, + Smn.
(25)
An exact computer realization of this method for analyzing the whole J-dimensional space is a time and memory consuming problem. The proposed approxima-
226
PAVEL A. CHOCHIA
FIG. 12. Color image enhancement: (a) the original image; (19) image enhancement by means of Eq. (25), a = 3; (c) a section of the original image with an enlargement; (d) the same section after enhancement.
tion of this method consists of an independent local contrast for each one of J color components by means of A2. The result of such processing of a color image is presented in Fig. 12. One can see that the image became dearer, colors became more saturated, and a number of texture features are revealed which were not observable before processing. 8. DETECTION OF THE DEFINED AREA OBJECTS
We consider the task of detection of the objects corresponding to domains with the area ( N 0 in one of the following interpretations: (i) NJ is greater than any value Q; (ii) inside the interval 0 1 < NJ < Q2; and (iii) NJ is less than Q. We should fimit our task by stipulating that the distance between domains must be greater than L, where L is the fragment size, and that the value of elements in the detected domain ( x ~ , ) is greater or less than all other elements x i ~ Wren \ U j in an analyzed fragment. The last condition is needed for the assumption that elements of required domains are found in the left end or in the right end of the fragment histogram H,,,. For NS > Q it is necessary to assume that the image have a number of small domains U 1. . . . , UU,... inside the large one U °. In this case one can use Algorithm A1 with the threshold np = Q. Then the background component, Sin,, will be a desired result, because it loses the small domains U j with N i < np. Figure 13a shows an X-ray image, obtained for a nondestructive test (here is a detail with some cavities). For detection of cavities with sizes greater than given Q, Algorithm A1 is used. In Fig. 13b, one can see the result. F o r detecting the objects with Q1 < NJ < Q2, one can define two different thresholds n 1 = Q1, n2 = Q2 and use their similar one np in Algorithm A1, with n 1 for a neighbourhood Vm, and n 2 for a fragment Win,. Using the threshold algorithm, like (21), for a neighbourhood I'm, (and its histogram H ~ , ) with R 1 = R ( n l / L ~ ) and R 2 = R(1 - nl/LZo), one can get new value Xm, instead of if,,,. Note, that £,,n lost the small domains with N i < Qx. Then 2,~, has been used in A2 with R 1 = R ( n E / L 2) and R E = R(1 - nE/L2w), obtained by the fragment
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
227
FIG. 13. Object detection with N j > Q: (a) the original X-ray image with cavities; (b) image processing by Algorithm A1. Detected cavities have areas greater than the given one.
histogram HWn. Here s,,n lost the middle domains with NJ < Q2. To detect the domains with Q1 < NJ < Q2 the following thresholding is required: 1,
if I~,~, - ~mnl > 6;
(26)
~Y-mn= O, ifl~Cmn--~mn I <3. Some examples are proposed in Fig. 14. The airphoto (Fig. 10a) is developed by (26) with different parameters (Q1, Q2)- In Fig. 14a Q1 and Q2 are chosen sufficiently small and detected objects also have small local areas. In Fig. 14b the parameters are increased and greater area objects are detected. Note, that the area of the domain is considered here as the "local" area, which is found inside the analyzed fragment.
FI6. 14. Object detection with small (a) and large (b) local areas in the airphoto (Fig. lOa).
228
PAVEL A. CHOCHIA
FXG.15. Contourdetection: (a) the originalimage; (b) the background component Smn;(C) contour detection in the originalimage; and (d) in the background component.
For detecting objects with N -/< Q, one can assume Q1 = 0. 9. CONTOUR DETECTION One of the traditional problems in image processing is contour detection. Many different methods for its solution are known, the majority of which are considered in [1]. Here we propose a new approach. As discussed in Section 2, the image consists of a set of domains separated by contours. In accordance with the image model the information about contours is imbedded in the background component Sm,, Eq. (1), which is smoothed inside the domains, and has the gaps in the contours. The component t,~, contains a texture, small details, and a noise. It can only bring distortions in contour detection. Thus it is evident that the image should be previously decomposed and only then the contour detection in the background component s,~, should be used. After picking out s,~n from the image (Section 5) the contour detection can be reached by any known contour operator. By varying the parameter np in algorithm A1 one can change the minimal sizes of the detected details. This approach corresponds to the one in Section 8. As an example, Figs. 15c, d show some results of detecting the contours by means of the well-known Robert's operator [1]
Ymn = [(Xmn- Xm+l,n+l) 2 @ (Xm,.+l- Xm+l,n)2] 1/2
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
229
in the original image (a) and in the background component, Smn , after decomposition (b). In image 15c one can see both the number of false pulses due to the noise and the destruction of the continuous lines. Image (d) is rather more legible. 1. 2. 3. 4. 5.
6.
7. 8. 9. 10. 11.
12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22.
REFERENCES W. K. Pratt, Digital Image Processing, Wiley, New York, 1978. J. W. Tukey, Exploratory Data Analysis, Addison-Wesley, Reading, MA, 1971. H. A. David, Order Statistics, Wiley, New York, 1970. P. J. Huber, Robust Statistics, Wiley, New York, 1981. B. I. Justusson, Median filtering: Statistical properties, in Two-Dimensional Digital Signal Processing. Part 2. Transforms and Median Filters (T. S. Huang, Ed.), Topics in Applied Physics Vol. 43, pp. 161-196, Springer-Veflag, New York, 1981. S. G. Tyan, Median filtering: Deterministic properties, Two-Dimensional Digital Signal Processing. Part 2. Transforms and Median Filters (T. S. Huang, Ed.), Topics in Appl. Phys. Vol. 43, pp. 197-217, Springer-Veflag, New York, 1981. N. C. Gallager and G. L. Wise, A theoretical analysis of the properties of median filters, IEEE Trans. Acoust. Speech Signal Process. ASSP-29, No. 6, 1981, 1136-1141. T. S. Huang, G. Y. Yang, and G. Y. Tang, A fast two-dimensional median filtering algorithm, IEEE Trans. Acoust. Speech Signal Process. ASSP-27, No. 1, 1979, 13-18. V. Kim and L. Yaroslavskii, Rank algorithms for picture processing, Comput. Vision Graphics Image Process. 35, 1986, 234-258. R. A. Hummel, Image enhancement by histogram transformation. Comput. Graphics Image Process. 6, No. 2, 1977, 184-194. T. P. Belikova, M. A. Kronrod, P. A. Chochia, and L. P. Yaroslavsldi, Digital processing of Mars surface photographs transmitted by automatic interplanetary probes "Mars-4" and "Mars-5," Space Res. 13, No. 6, 1975, 898-906. [Russian] I. Scollar, B. Weidner, and T. S. Huang, Image enhancement using the median and the interquartile distance, Comput. Vision Graphics Image Process. 25, No. 2, 1984, 236-251. J. K. Yah and D. J. Sakrison, Encoding of images based on a two-component source model, IEEE Trans. Commun. COM-25, No. 11, 1977, 1315-1322. P. A. Chochia, Digital pulse filtering on the television images, Commun. Technol. Set. Television Technol. No. 1, 1984, 26-36. [Russian] P. A. Chochia, Two-component statistical model of an image fragment, in Image Processing and Remote Sensing, Reports of All-Union Conf. Novosibirsk, 1984, pp. 60-61. [Russian] R. M. Haralic and L. Watson, A facet model for image data, Comput. Graphics Image Process. 15, No. 2, 1981, 113-129. S. Nishikava, R. J. Massa, and J. C. Mott-Smith, Area properties of television pictures, IEEE Trans. Inform. Transmiss. IT-11, No. 3, 1965, 348-352. J. S. Lee, Digital image smoothing and the sigma filter, Comput. Vision Graphics Image Process. 24, 1983, 255-269. B. R. Frieden, A new restoring algorithm for the preferential enhancement of edge gradients, J. Opt. Soc. Amer. 66, No. 3, 1976, 280-283. C. A. Pomalaza-Raez and C. D. McGillem, An adaptive nonlinear edge-preserving filter, IEEE Trans. Acoust. Speech Signal Process. ASSP-32, No. 3, 1984, 571-576. P. A. Chochia, Air-Space image enhancement methods using fragment histogram, Earth Res. from Space, No. 6, 1985, 66-78. [Russian] G. A. Mastin, Adaptive filters for digital image noise smoothing: An evaluation, Comput. Vision Graphics Image Process. 31, No. 1, 1985, 103-121.
NOTE Image Enhancement Using Sliding Histograms PAVEL A. CHOCHIA
Institute of Information Transmission Problems, The USSR Academy of Sciences, Moscow, USSR Received April 26, 1986; accepted January 7, 1988
Image enhancement methods based on local histogram analysis are discussed. It is shown that this problem is closely connected to that of contour preserving image smoothing and may be solved by a preliminary image decomposition. A method for image decomposition into background and texture components is proposed. The results obtained are used for image smoothing, enhancement, object and contour detection, and color image enhancement. © ]988 Academic Press, Inc.
1. INTRODUCTION
In this paper we consider image enhancement by local contrast amplification of both fine details and texture [12]. We show that this task is closely related to the image smoothing problem and may be solved by a preliminary image decomposition into two different information components. The aim of this paper is to develop an image smoothing method using the probability distribution function for a sliding fragment, and then to use it in some image processing tasks. This problem is based on robust estimation of sliding histogram parameters using order statistics [2-4]; the median filter, suggested by Tukey [2], is the most useful. Statistical properties and contour preserving facilities of the median are well known [5-7], and fast algorithms have been developed using sliding histograms to accelerate its computation [8]. For image smoothing and enhancement, a model of local image fragments is proposed. Using this model, we will show that both local averaging and local median methods achieve unsatisfactory image smoothing results. The new proposed image smoothing method is based on the M-type estimators used in order statistics [4]. This method makes it possible to decompose the image into two components. This operation leads us to a satisfactory solution of the image enhancement problem and some other problems. A digital monochromatic image is an M x N matrix X = [xm, ], and the image value x,~, represents the optical density (logarithm of the luminosity), quantized into K levels, i.e., 0 < x,,, < K - 1. In the following two different sets of elements in the vicinity of x will be used (Fig. 1). To differentiate them we use the terms "neighbourhood," V,,, (for inner) and "fragment," Wm, (for outer). Both of them are just moving rectangular windows, but they have different sizes. In real experiments they were squares from 3 × 3 up to 9 x 9 pixels (neighbourhood) and from 11 X 11 up to 65 X 65 pixels (fragment). The fragment histogram, Hmn = {hmn(k)), 211 0734-189X/88 $3.00 Copyright © 1988 by Academic Press, Inc. All rights of reproduction in any form reserved.
212
PAVEL A. CHOCHIA Wren
FIG. 1.
An image fragmentcontaining three domains (U1, U2, U3).
is a probability distribution for the magnitudes of the elements: E h , n n ( k ) = 1,
k=(0 ..... K-
1).
k
We assume that a discrete random quantity 4, 0 < ~ < K, has a normal distribution N(# a aS), if its distribution function F ( k ) coincides at all points k, where ~ can be nonzero, with the distribution function of a continuous normal random quantity ~/, having the same mean value and dispersion. 2. MODEL OF AN IMAGE FRAGMENT It is assumed that the original image x,,, is the sum of two statistically independent components [13, 14]: Xmn = Srn n -~ tmn ;
(1)
a background component Stun, which determines the mean intensities of large size image details, and a texture component, tmn , which contains information about the texture and small details of the image. In order to separate Stun, one has to know the difference between statistical properties of the background and texture components. Therefore before starting the analysis and the development of the transformation methods, we must choose a model of the image, enabling us to take into account the statistical correlations of the image elements within a fragment. We assume that the image consists of connected domains U 1, U 2. . . . . u N , . . . , separated by contour boundaries. According to the model of a neighbourhood [15], an image is a sum of two nonoverlapping sets: contour elements and internal elements of the domains. For elements of the neighbourhood x m , ~ Vm, (r = 1 . . . . . R), belonging to one and the same set, one can draw a plane (say, by means of the method of least squares) so that the relation holds Xmn
= ~mn
-]- ern(armnP r + }~mn) + (1 -- ern)(fl,~np" + ~/rn) ,
(2)
where e,n nr is a contour mask (e = 1 for the contour element and e = 0 for internal one), /X,nn is a constant determining the mean value of the neighbourhood elements, pr is the distance from the center (m, n) of the neighbourhood to the r t h element,
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
213
a r . and fl~, are random quantities which determine a slope of the plane in the direction of the rth point: a~,, (for contour) and B~, (for internal) are equal to the cosine of the angle between the vector from (m, n) to the rth point, and the plane gradient vector G,, n. ~ , and 7/~,, are independent random quantities having normal distributions N(0, o{) and N(0, ~ ) , respectively. The central element is Xrnn = ~mn "{- emn~mn +
(1
-
emn)~mn.
The one-component version of this model is close to model [16]. This model describes statistical correlations between elements within a single neighbourhood, limited to 3-5 discretization steps. However, if a contour line does not separate the points with coordinates (m, n) and ( m + e, n + 3), then there are statistical dependences between the values of /~m, and /~m+~,,+n for sufficiently large e, 3 (up to 20-40 and larger), since these points belong to the same image domain [17]. This fact can be used in constructing the model for a sufficiently large fragment. Let Wren be a fragment with its center at (m, n), which is divided by contour boundaries into J background domains U1,..., U J. Points of a domain U j belong to the same extended detail of the fragment and form a connected region x~ ~ U j c Wren (i = 1 , . . . , I J) (Fig. 1). Certain quantifies #i correspond to the elements x i. We shall assume that within a single domain U J, at distances not exceeding the fragment size, values of /~i are described with sufficient accuracy in terms of a piece-wise constant model ' + ~ i = S rnn
(3)
where s~, is the mean intensity of the domain U j inside the fragment Wrnn,and x/ are (perhaps, dependent) random quantities with a normal distribution N(0, o~). Such a model is quite reasonable for most real images. Hence an element x / o f the domain u J is given [15] by
X! = SJn ~- I£i q- ~i ~---SJmn+ T/.
(4)
The physical meaning of two added terms is quite different. The first component Ki determines the texture and small details of the image, the values of its elements are substantially correlated. The second component 7/i represents non-correlated Gaussian noise. As Ki and ~i are mutually independent random quantities with normal distributions N(0, a~) and N(0, o~), respectively, their sum ~'i = xi + ~1i is also a random quantity with the normal distribution N(0, 02), where o~ --- of + a~. The central element of the fragment in the domain U 1 is 1 1 1 _if_ I£mn 1 _]_ ~ 1 n = Sm n + ,r;nn. Xmn ~- Stun
(5)
It follows from (4) that if a fragment is in a single domain of the image, its histogram is unimodal. If the fragment covers parts of two or more ( J ) domains, the histogram is J-modal. Positions of the modes are determined by mean intensities of the domains ( S 1, S 2. . . . , SS), and powers of the modes are proportional to fractions of the area of the fragment occupied by corresponding domains (Fig. 2). Filling the
214
PAVEL A. CHOCHIA h(k)
k
S1
s3
S2
FIG. 2. The histogram of the image fragment (Fig. 1).
gaps between the modes depends on the mode dispersions (OrJ ) 2 and on intermodal distances, which determine overlapping of the distribution tails, and on a fraction of the contour elements in the whole distribution. Thus the problem of the evaluation of Sm. for image point (m, n) can be solved as an estimation of a mean intensity level s.,nl in the central domain U 1 of the fragment Wm.. This problem will be discussed in Section 5. We now review the known image enhancement methods. 3. IMAGE PROCESSING METHODS REVIEW
Processing methods providing an enhancement of local contrasts in an image can be characterized by two properties: linear/nonlinear methods and fixed/adaptive methods. The first characteristic means the linearity in the contrast amplification, and the second one is concerned with the varying of this amplification factor within the image field. All the methods are based upon the use of the fragment histogram H,,,. Assume that the fragment is a rectangle of size L × L, L = 2l + 1. The estimation employed is either the mean value of fragment elements K-1
A(Wm. ) = Y'. k h . , . ( k ) = L k~O
-2
m+l
n+l
Z
E
xij,
(6)
i=m--I j=n-I
or order statistics of the type R(q), 0 _< q < 1, defined by y-1
Rm.(q) = y,
y
if ~_, hm.(k) < q < ~_, hm.(k ). k=O
(7)
k=0
In some cases the median is used, namely m e d ( W ~ ) = Rm~(0.5 ).
(8)
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
215
1. Linear Fixed ( L F ) Methods The image is decomposed into two components: Xm. = Stun + tin.; various transformations are given by the formula : Ym. =atmn + bsm. + C.
(9)
The parameters a, b, c are constant throughout the image; in most cases c = (1 - b ) / 2 . Two methods for the evaluation of Sm, are used: the first is the average (6) over the sliding fragment [11, 12], and the second is the median (8) for the same fragment [12]. If a = b, then we obtain the standard linear correction formula Ymn = aXmn + C.
(10)
2. Linear Adaptive ( LA) Methods The formula used for the processing is similar to that of the LF methods (9), but the parameter a, which is the amplification factor of local contrasts, is determined from an estimate of variation of the fragment histogram Hmn [12]: Ymn = amntmn + bsrnn At- ¢"
(11)
Two variants are used to evaluate the parameter amn: (i) If Sin. is the average, then a,.. is determined by the dispersion of the element intensity distribution over the fragment am. = O r J ( % r + 8).
(12)
Here O2r is the dispersion of the intensity distribution for the original image in the fragment Wm.; ozs is the desired resultant dispersion; and 3 is a parameter restricting the contrast amplification. (ii) If Sm~ is the median, then amn is determined by the interquartile distance AR = R(0.75) - R(0.25), a,n . = ARres/(ARor + 3),
(13)
where R(0.25) and R(0.75) are magnitudes of the first and the third quartiles, respectively, and 3 is a parameter, as in (12). 3. Nonlinear Adaptive (?CA) Methods A known method of this type is the sliding histogram equalization method [9]. It is represented in terms of the order statistics Ym. = R m. (Zm.),
where z,.. = x . , . / K .
(14)
It is not difficult to extend this method, obtaining a sliding histogram modification. Let G ( k ) be the resulting probability distribution, and G - l ( z ) be the inverse
216
PAVEL A. CHOCHIA
function, i.e., G-I[G(k)] = k. Then the sliding histogram modification is given by
Ymn= G-t[Rmn(Zmn)].
(15)
The function G(k) can be chosen on the basis of certain requirements of the processing result. Variants employing a similar modification are known for the histogram of the whole image [10].
4. Nonlinear Fixed (NF) Method Methods of this type have not been found in the literature and are mentioned for completeness. They can be constructed, for example, with the formula
Ym, = f(tm,) + bsm, + c. Here f ( x ) can be an arbitrary nonlinear function. It is important to note that most methods considered use a preliminary separation of the image into two components, stun and tmn, for the following amplification of the component tin,. In the present work we restrict ourselves by studying the LF methods. As was mentioned above, the known LF methods differ in their determination of Sm~. Let us consider the accuracy in the determination of the mean intensity level Sm,l for the central domain U 1, containing the element Xm,, when the average or the median over the fragment W,,, is used. 4. AVERAGE AND MEDIAN OVER THE FRAGMENT
If the fragment Wm, is within a single domain of the image, then according to the model (4), the fragment histogram is normally distributed, and the expected values of the average M{A(W)} and the median M(med(W)} are equal to the mean intensity level S 1 for the image domain. If the fragment covers several image domains, the expected values M ( A ( W ) } and M{med(W)} will be equal to some intermediate quantities within the interval (Smi~, Smax). In this case the set of fragment elements can be regarded as a combination of the elements of J normal distributions, N( S 1, o2) .... , N( S J, o2). Let Pj be the probability that an arbitrary element of the fragment belongs to the domain UJ: EPj = 1. Then the expectation value of the average over the L × L square fragment is
M ( A ( W ) } = MIL-2~_,Xm. 1 = ~_,PiS/. \ W
"
(16)
j
The expected value of the median is Z0, which is the root of the equation
EP/F+(z) = 0 . 5 ;
(17)
J
Fj(z) is the distribution function for dements of the domain U/. Let us compare the results for J = 2. For simplicity we set o 1 = 02, AS = S 1 - S 2. The Pj dependences of M{A(W)} and M(med(W)} for AS = 20, 3o, and 40 are shown in Fig. 3a. It is seen, that for AS = 20 the average and the median have a
217
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS b A(W) med(W)
A(W) mecl(W)
4a
3a
3a 2a 2o
0
~ P2 0.5
I
0 ~
' 0.5
t 1
~- P2
FIG. 3. Dependences of the average (dashed line) and the median (solid fine) over a fragment for two-mode distribution: (a) a 1 = 02; I - AS = 2a, 2 - AS = 30, 3 - AS = 40; (b) AS = 30; 2 - o2 = 01,4--a 2=0.5al,5-o 2=0.
small difference. As AS is increased, M{med(W)} is almost invariably near the ends of the interval, where P2 < 0.3 or P2 > 0.7, but a step appears near P2 = 0.5, while the expected value of the average, M{ A(W)}, is always a linear function. The dependencies presented here have a clear physical meaning. Suppose that there are two domains in the image, separated by a contour difference, and intensities in each domain are described with the normal distributions N(S 1, 02) and N(S 2, 02). Then the graphs show variations of the average and the median for a fragment along a line perpendicular to the contour. Evidently, for % 4:02 the slope of the graph for the median is determined by the dispersion 02 for small P> and by o22 for large P2. The graph for the average over the fragment will not be changed (Fig. 3b). Results of processing for a real image (Fig. 4a) performed by means of (9), determining Stun from the average and the median, are shown in Figs. 4b, c. With a
FIG. 4. Results of the image processing by means of known methods: (a) a part of the original image (Fig. 10a); (b) processing by means of Eq. (9), Stun is the average over a sliding fragment; (c) the same, b u t stun is the median; (d) sliding histogram equalization by Eq. (14).
218
PAVEL A. CHOCHIA
FIt. 5. Test image processing (image dimensions 128 x 384 elements, fragment dimensions 33 × 33 elements): (a) the original image and the graph for one of its lines; (b) the average (6) and graph for the same line; (c) the median (8).
sufficiently high enhancement of the fine texture in flat background regions, a number of distortions have arisen during the processing: an over-contrast of the contour steps in the case of the average, and of the detail angles in the case of the median. In Fig. 4d we present for comparison the result obtained using the sliding histogram equalization algorithm (16); here is an over-contrast of the texture in flat image regions. Distortions are seen more dearly in the results of processing a test image (Fig. 5a). Magnitudes of the local average and the median are presented in Figs. 5b, c. It is seen that the transformations appreciably changed contours and shapes of the objects. The transformation using the average has smeared the contours uniformly, while the transformation using the median conserved extended linear contours but erased angles of the details. These distortions appear, because, in a fragment having the angle of the detail at its center, the total number of object dements is less than the number of background dements. Using such distorted estimates in (9) results in an over-contrast of contours, details, etc. Thus we can conclude that both the average and the median for a fragment lead to unsatisfactory estimates for the mean intensity level of the background in the central domain, Stun, 1 if the fragment covers two or more different image domains. 5. A METHOD FOR DETERMINATION OF BACKGROUND INTENSITY
Now we try to find a method suitable for determining the mean intensity level sin, , using a fragment histogram H,n,. A typical histogram for a fragment containing several domains UJ ( j = 1,..., J) is shown in Fig. 2. It consists of several modes, corresponding to domains in the fragment. An dement Xm, belongs to the central domain U 1 of the fragment, unless it is a contour element. The problem of evaluation of the background intensity, or equivalently the mean intensity level of the central domain S 1 is reduced to determination of the center of the correspond-
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
219
i n g distribution (mode) in the histogram. To solve this problem one has to: (i) choose the appropriate mode in Hm~ and (ii) determine the mean value of the distribution. By choice of any mode in the histogram, we determine an indication of the center of the corresponding distribution with some accuracy. It should be noted that by choosing the mode nearest to Xm. one does not always get correct results. Actually, in real images the summed dispersion o~ in (4) may be comparable to the difference of the mean intensities of adjacent domains IS g - SJl, and it is quite probable that the mode nearest to value x , ~ would be attributed to another domain. Analysis of real images shows that inside uniform domains o~ are less than % by a factor of 2 to 3; hence for choosing an appropriate mode one should use the value /*m~, Eq. (2), which is averaged over the elements of a neighbourhood Vm~, belonging to the same domain as xm~. For smoothing the image we may use a method, based on M-type estimators, used in order statistics [4, 20], which is a function of an ordered elements set ( ~ } = ~1 -~< ~2 ~
"'"
-~< i n ,
where w~ depends on the 4. If for any x i ~ Vmn we choose wi = 1 for Xmn < Xm~ + 8 and w~ = 0 for others x~, then we obtain the filter
Xmn =
WiXi
i=1
Wi.
li=l
(18)
This filter is used in [18] for image smoothing and is called a "sigma filter." Parameter ~ is recommended as 8 = 20, where o 2 is the noise dispersion. As was stated in [22], the sigma filter is one of the best edge-preserving noise cleaning methods. It follows from (5) that for a point (m, n) inside the domain U -i, Xmn has a distribution close to the normal one N(SJ, o~ + ~/Y.wi) , where S j is the mean intensity of the background component in UL To determine S 1 we shall construct an analogous sigma filter for a fragment. Fixing a value of A, we take points in the fragment W,,, (or equivalently, points of the histogram Hm,) which are within the window (Xmn -- A, Xmn "q- m), and calculate for them the average A(W; ~), or the median med(W; ~). Let us compare the accuracies of the estimates A(W; ~) and med(W; ~). Accordhag to the model (4), elements in one domain have a normal distribution N(S j, o~). Suppose that the values of the elements of neighbour domains are sufficiently remote and their influence is negligible. For simplicity we put a = 1 and use the notation I N T ( a , b) = (2~r) -1/2 ~bexp(-o2/2) do.
220
PAVEL A. CHOCHIA
/ .
med(W ;~)
A:
a
A=1.4o I
A=2a
0.5a
A=3a
,
.
0.5a
a
,
1.5a
.~ "~
2a
FIG. 6. The average A (W; ~) and the median med(W; Y) versus ~ (curves 1-4 and 5-8, respectively). The results have been obtained by means of the sigma filter (18) with A = a, 1.4o, 2a, 30.
The expected value for the average of elements within the window (~ - A, ~ + A) is
M{ A } = (2~r)-l/z ff+2v2exp(-vz/2) dv/INT(.~ =
(2~')-i/2[exp(- (.2×INT(£-
A)2/2)-
A, y + A)
exp(-(ff + A)2/2)]/
A, ff + A).
(19)
The expected value for the median, M{med}, for the same window is a root of the equation I N T ( ~ - A, M ( m e d } ) = INT(ff - A, y + A ) / 2 ,
(20)
In Fig. 6 we present the Y dependences of M(A} and M ( m e d ) for A = o, 1.4o, 20, 30. The ratios a = M{ A }/ff and a = M ( m e d } / Y determine the convergence velocity for the corresponding estimate. The comparison between the graphs shows that for all values of A, med(W; x) converge more rapidly than A(W; x). As one would expect, a increases with A from 1.4-1.5 at A = o to 3-5 at A = 20 and 15-30 at A = 30. It is noteworthy that if overlapping modes are present in the real distributions, the probability that an appreciable fraction of the second mode elements are inside the window (E - A, y + A) increases with A. Therefore it is unreasonable to take A > 20. An iterative procedure, where at the k t h step ~k+l is evaluated with the known xk, enables one to increase the initial convergence velocity a up to about a k, where k is the number of the iteration steps. For A = 20 and k = 3 we obtain a 3 = 100. In the determination of the mean intensity level in the central domain we have not taken into account its dimensions, i.e., the number of elements nl contained in U 1. It is evident, however, that with a decrease in n i a threshold value n e may be
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
221
attained, such that all the elements of U I must be, in fact, considered as a texture against the background of another domain UJ enclosing U 1. A choice of the threshold np is a nonformalized problem. It depends on the concrete subject, the image discretization step, and the final aim of processing. Therefore we propose to treat np as a parameter, proceeding in the following way. Let us suppose we have a fragment histogram Hmn, a parameter np, and an initial estimate for the mean intensity level of the background ~, calculated over a neighbourhood I'm,. The fragment histogram Hm~ can be used to determine the order statistics Rl(np/L 2) and R2(1 - np/L2). The basic value xl for the first iteration is given by the thresholding algorithm:
{
~,
X1 =
R1,
R2,
if R 1 < ~ < R 2, if X _< R1, if ~ > R 2.
(21)
This algorithm provides the desired result in the overwhelming majority of cases: if the fragment covers only two domains--the background (environment) U; and the central (small) domain U 1, as well as in situations where S 1 is either larger than the maximum or lower than the minimum intensity for the other domains of the fragment. Thus we derive the following algorithm for determination of the mean intensity stun for the central domain U 1 of fragment Wren. ALGORITHM 1 (A1). 1. Calculate histograms for the neighbourhood, Hm~,, centered at the point (m, n). 2. Evaluate ~ over the neighbourhood histogram
H°m~,and for the fragment, H~n, using the sigma filter
(18). 3. Evaluate X1 at a fixed rip, using algorithm (21). 4. Median (20) with the window (xk - A, ~k + A) is calculated by iteration from H,~,. The obtained quantity is considered as the desired s,~. Results of processing the test image by means of A1 at various np are shown in Fig. 7 (the test image is given in Fig. 5a). In the first case np was rather large, and
--1 F--L-
FIG. 7. Test image smoothing (original image is given in Fig. 5a) by means of A1. Neighbourhood dimensions 5 × 5 elements, fragment dimensions 33 × 33 elements: (a) n; = 80; (b) np = 10.
222
PAVEL A. CHOCHIA
FIG. 8. Distortion of contour step: (a) appearance of the double contour line when Stun is calculated by means of A1; (b) distortions disappear if modifiedvariant A2 is applied.
small details are filtered out of the image, as well as the noise. In the second case (at small np) the image is smoothed, but all details are retained. The accuracy of the reproduction of the original image (without noise), calculated with the r.m.s. criterion, was 0.3% after the first iteration. The graphs show a good smoothing of the noise, while the contour differences were conserved completely. 6. TEXTURE ENHANCEMENT Texture enhancement is performed by setting the difference tmn = Xm, -- S,~n in (9). Here Smn is the mean background intensity, obtained by means of A1. However, as experiments have shown, sometimes at domain boundaries a monotonous contour step is broken, and a false contour appears (Fig. 8a). One reason for these distortions is the fact that dements of the contour are ascribed by A1 to one of two neighbour domains of the image. Let us consider the origin of these distortions in more detail. Suppose we have two domains, U 1 and U 2, separated by a continuous contour step (Fig. 9a), their mean intensities are S 1 and S 2, and S 1 > S 2. The contour elements have intermediate intensities in the interval (S 2, $1). The result of the image processing by means of A1 will be as follows. Points to the left of the middleof the contour step will be ascribed to the region U 1 and after the processing their values will be close to S 1, while points to the right will have values close to S 2. The graph for this smoothed image (Sm,) is shown with dashed line (Fig. 9b). The same effect was discussed by Frieden [19]. Setting the values of Sin, and tin, into (9) with a = 3 we obtain the result shown in Fig. 9c. It is seen clearly that a " b u m p " appears in the contour, and, as a consequence, an additional line along the contour emerges in the transformed image. To avoid such distortions one should modify the procedure of determining Sm, for contour points. Evidently, if (m, n) is a point on the contour step, then the majority of fragment elements (in flat domains) will not be connected statistically with the d e m e n t Xm,. Therefore in the determination of Sm, for the contour point one should take into account only its nearest dements within the neighbourhood
223
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS ix(n)
2
s~ Ik l~.
Slt x(n) S2
FIG. 9. For the explanation of origin of the contour distortions: (1) the graph of the original contour step; (2) the value of Stunis calculated by means of A1; (3) the value of Ymn is calculated by means of Eq. (9).
Vmn. A sufficiently appropriate way to determine stun in this case is to take Xmn' calculated at the second step of A1. The proposed modification of A1 [21] is as follows. ALGORITHM 2 (A2). 1. Evaluate sin, by A1. 2. Identify contour points. To this end calculate an amplitude gm~ of the plane gradient G,~ used in model (2), by the formula gmn= [($m-l,n-lWSm-l,n+Sm-l,n+l--Sm+l,n-l--Sm+l,n--Sm+l,n+l)
"t"(Sm-1, n-1 "+ Sin, n-1 -1- Sm+l,n-1 --
Sm-l,n+l
--
Sin, n + l
--
2
Sm+l,n+I)2] 1/2
3. The value g,., is compared with the threshold & If gin, < 8, then Xm, is the dement inside a domain; then put g,~, = Sm,. If gin, > 8, then Xm, is the element of a contour overfall; then set g,,, = Xm," The obtained quantity g,,, is considered as the mean background intensity at the central point (m, n) 4. Evaluate the texture component, tm, = Xmn -- gmn" 5. Perform the image enhancement by transformation Ymn = atmn + bsmn + c.
(22)
Figure 10 shows the result of image enhancement by means of A2. It is seen that in the resulting image the contrast of texture and small details is enhanced considerably with respect to the background, while mean intensities of the extended details are not changed and significant distortions have not appeared. 7. COLOR IMAGE ENHANCEMENT
The method of local contrast enhancement presented above for black-and-white pictures can be extended by analogy to processing of color or multi-zonal images. In
224
PAVEL A. CHOCHIA
FIG. 10. The results of airphoto processing by means of A2: (a) the original image (512 × 512 elements); (b) local texture enhancement by means of A2.
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
225
aT s1
G
RE FIG. 11. Three-dimensional color histogram H(R, G, B) and multiplication of the original texture vector T up to a T.
this case any value of the image element must be considered as a point or a vector in a J-dimensional space of color coordinates (say, in the three-dimensional (R, G, B) space)
Xmn
Xmn ~- Stun "~ Kmn + ~mn = Stun + Tmn"
(23)
In this space a J-dimensional color histogram H ( k l , . . . , k j) for a fragment Wren is calculated at every point of which the probability of a given color value of an element is reflected. The points with non-zero probability condense into clusters (Fig. 11) corresponding to domains of a fragment (and to modes of the ordinary one-dimensional histogram). To define the background intensity Sin, of a color image (by analogy with the determination of background intensity of monochrome image by means of A1), it is necessary to separate a region v1 in the color histogram, corresponding to a central domain U 1 of the image fragment Win,,, and to find the location of its center (vector S 1, Fig. 11). This vector is considered as the value of the smoothed color image Stun at the point (m, n). Having determined the texture v e c t o r Trn n = X m n -- Stun , we obtain, in analogy to (22): Ymn = aTmn + bSmn --I- C.
(24)
Alternatively, conserving the original mean colors: Ym,, = aT,,,,, + Smn.
(25)
An exact computer realization of this method for analyzing the whole J-dimensional space is a time and memory consuming problem. The proposed approxima-
226
PAVEL A. CHOCHIA
FIG. 12. Color image enhancement: (a) the original image; (19) image enhancement by means of Eq. (25), a = 3; (c) a section of the original image with an enlargement; (d) the same section after enhancement.
tion of this method consists of an independent local contrast for each one of J color components by means of A2. The result of such processing of a color image is presented in Fig. 12. One can see that the image became dearer, colors became more saturated, and a number of texture features are revealed which were not observable before processing. 8. DETECTION OF THE DEFINED AREA OBJECTS
We consider the task of detection of the objects corresponding to domains with the area ( N 0 in one of the following interpretations: (i) NJ is greater than any value Q; (ii) inside the interval 0 1 < NJ < Q2; and (iii) NJ is less than Q. We should fimit our task by stipulating that the distance between domains must be greater than L, where L is the fragment size, and that the value of elements in the detected domain ( x ~ , ) is greater or less than all other elements x i ~ Wren \ U j in an analyzed fragment. The last condition is needed for the assumption that elements of required domains are found in the left end or in the right end of the fragment histogram H,,,. For NS > Q it is necessary to assume that the image have a number of small domains U 1. . . . , UU,... inside the large one U °. In this case one can use Algorithm A1 with the threshold np = Q. Then the background component, Sin,, will be a desired result, because it loses the small domains U j with N i < np. Figure 13a shows an X-ray image, obtained for a nondestructive test (here is a detail with some cavities). For detection of cavities with sizes greater than given Q, Algorithm A1 is used. In Fig. 13b, one can see the result. F o r detecting the objects with Q1 < NJ < Q2, one can define two different thresholds n 1 = Q1, n2 = Q2 and use their similar one np in Algorithm A1, with n 1 for a neighbourhood Vm, and n 2 for a fragment Win,. Using the threshold algorithm, like (21), for a neighbourhood I'm, (and its histogram H ~ , ) with R 1 = R ( n l / L ~ ) and R 2 = R(1 - nl/LZo), one can get new value Xm, instead of if,,,. Note, that £,,n lost the small domains with N i < Qx. Then 2,~, has been used in A2 with R 1 = R ( n E / L 2) and R E = R(1 - nE/L2w), obtained by the fragment
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
227
FIG. 13. Object detection with N j > Q: (a) the original X-ray image with cavities; (b) image processing by Algorithm A1. Detected cavities have areas greater than the given one.
histogram HWn. Here s,,n lost the middle domains with NJ < Q2. To detect the domains with Q1 < NJ < Q2 the following thresholding is required: 1,
if I~,~, - ~mnl > 6;
(26)
~Y-mn= O, ifl~Cmn--~mn I <3. Some examples are proposed in Fig. 14. The airphoto (Fig. 10a) is developed by (26) with different parameters (Q1, Q2)- In Fig. 14a Q1 and Q2 are chosen sufficiently small and detected objects also have small local areas. In Fig. 14b the parameters are increased and greater area objects are detected. Note, that the area of the domain is considered here as the "local" area, which is found inside the analyzed fragment.
FI6. 14. Object detection with small (a) and large (b) local areas in the airphoto (Fig. lOa).
228
PAVEL A. CHOCHIA
FXG.15. Contourdetection: (a) the originalimage; (b) the background component Smn;(C) contour detection in the originalimage; and (d) in the background component.
For detecting objects with N -/< Q, one can assume Q1 = 0. 9. CONTOUR DETECTION One of the traditional problems in image processing is contour detection. Many different methods for its solution are known, the majority of which are considered in [1]. Here we propose a new approach. As discussed in Section 2, the image consists of a set of domains separated by contours. In accordance with the image model the information about contours is imbedded in the background component Sm,, Eq. (1), which is smoothed inside the domains, and has the gaps in the contours. The component t,~, contains a texture, small details, and a noise. It can only bring distortions in contour detection. Thus it is evident that the image should be previously decomposed and only then the contour detection in the background component s,~, should be used. After picking out s,~n from the image (Section 5) the contour detection can be reached by any known contour operator. By varying the parameter np in algorithm A1 one can change the minimal sizes of the detected details. This approach corresponds to the one in Section 8. As an example, Figs. 15c, d show some results of detecting the contours by means of the well-known Robert's operator [1]
Ymn = [(Xmn- Xm+l,n+l) 2 @ (Xm,.+l- Xm+l,n)2] 1/2
IMAGE ENHANCEMENT USING SLIDING HISTOGRAMS
229
in the original image (a) and in the background component, Smn , after decomposition (b). In image 15c one can see both the number of false pulses due to the noise and the destruction of the continuous lines. Image (d) is rather more legible. 1. 2. 3. 4. 5.
6.
7. 8. 9. 10. 11.
12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22.
REFERENCES W. K. Pratt, Digital Image Processing, Wiley, New York, 1978. J. W. Tukey, Exploratory Data Analysis, Addison-Wesley, Reading, MA, 1971. H. A. David, Order Statistics, Wiley, New York, 1970. P. J. Huber, Robust Statistics, Wiley, New York, 1981. B. I. Justusson, Median filtering: Statistical properties, in Two-Dimensional Digital Signal Processing. Part 2. Transforms and Median Filters (T. S. Huang, Ed.), Topics in Applied Physics Vol. 43, pp. 161-196, Springer-Veflag, New York, 1981. S. G. Tyan, Median filtering: Deterministic properties, Two-Dimensional Digital Signal Processing. Part 2. Transforms and Median Filters (T. S. Huang, Ed.), Topics in Appl. Phys. Vol. 43, pp. 197-217, Springer-Veflag, New York, 1981. N. C. Gallager and G. L. Wise, A theoretical analysis of the properties of median filters, IEEE Trans. Acoust. Speech Signal Process. ASSP-29, No. 6, 1981, 1136-1141. T. S. Huang, G. Y. Yang, and G. Y. Tang, A fast two-dimensional median filtering algorithm, IEEE Trans. Acoust. Speech Signal Process. ASSP-27, No. 1, 1979, 13-18. V. Kim and L. Yaroslavskii, Rank algorithms for picture processing, Comput. Vision Graphics Image Process. 35, 1986, 234-258. R. A. Hummel, Image enhancement by histogram transformation. Comput. Graphics Image Process. 6, No. 2, 1977, 184-194. T. P. Belikova, M. A. Kronrod, P. A. Chochia, and L. P. Yaroslavsldi, Digital processing of Mars surface photographs transmitted by automatic interplanetary probes "Mars-4" and "Mars-5," Space Res. 13, No. 6, 1975, 898-906. [Russian] I. Scollar, B. Weidner, and T. S. Huang, Image enhancement using the median and the interquartile distance, Comput. Vision Graphics Image Process. 25, No. 2, 1984, 236-251. J. K. Yah and D. J. Sakrison, Encoding of images based on a two-component source model, IEEE Trans. Commun. COM-25, No. 11, 1977, 1315-1322. P. A. Chochia, Digital pulse filtering on the television images, Commun. Technol. Set. Television Technol. No. 1, 1984, 26-36. [Russian] P. A. Chochia, Two-component statistical model of an image fragment, in Image Processing and Remote Sensing, Reports of All-Union Conf. Novosibirsk, 1984, pp. 60-61. [Russian] R. M. Haralic and L. Watson, A facet model for image data, Comput. Graphics Image Process. 15, No. 2, 1981, 113-129. S. Nishikava, R. J. Massa, and J. C. Mott-Smith, Area properties of television pictures, IEEE Trans. Inform. Transmiss. IT-11, No. 3, 1965, 348-352. J. S. Lee, Digital image smoothing and the sigma filter, Comput. Vision Graphics Image Process. 24, 1983, 255-269. B. R. Frieden, A new restoring algorithm for the preferential enhancement of edge gradients, J. Opt. Soc. Amer. 66, No. 3, 1976, 280-283. C. A. Pomalaza-Raez and C. D. McGillem, An adaptive nonlinear edge-preserving filter, IEEE Trans. Acoust. Speech Signal Process. ASSP-32, No. 3, 1984, 571-576. P. A. Chochia, Air-Space image enhancement methods using fragment histogram, Earth Res. from Space, No. 6, 1985, 66-78. [Russian] G. A. Mastin, Adaptive filters for digital image noise smoothing: An evaluation, Comput. Vision Graphics Image Process. 31, No. 1, 1985, 103-121.