June 1995
Pattem$;;gnition ELSEVIER
Pattern Recognition Letters 16 (1995) 653-666
A fast thresholding selection procedure for multimodal and unimodal histograms Du-Ming Tsai Department of Industrial Engineering, Yuan-ZeInstitute of Technology, 135 Yuan-TungRoad, Nei-Li, Taiwan, R.O.C. Received 10 October 1994; revised 7 December 1994
Abstract
In this paper, a simple and efficient histogram-based approach is presented for multi-level thresholding. It uses Gaussian kernel smoothing to detect peaks and valleys in a multimodal histogram, and uses a local maximum curvature method to detect points of discontinuity in a u&nodal histogram. The computational time will decrease as the desired number of thresholding levels increases. The performance of the proposed algorithm is compared with those of the widely applied between-class variance and entropy methods. Keywords: Multi-level thresholdiog; Gaussian smoothing; Local maximum curvature; Multimodal histogram; Unimcdal histogram
1. Introduction Thresholding has been a popular tool used in image segmentation. It is useful in separating objects from background, or discriminating objects from objects that have distinct gray-levels. Sahoo et al. (1988) have presented a thorough survey of a variety of thresholding techniques. Among those techniques, global, histogram-based algorithms (Glasbey, 1993) are widely used to determine the threshold, and they can be classified as parametric and nonparametric approaches. In the parametric approaches (Snyder, 1990; Weszka and Rosenfeld, 19791, the gray-level distribution of each class is assumed to have a probability density function; it is usually assumed to be a Gaussian distribution. One attempts to find an estimate of the parameters of the distribution that will best fit the given histogram data in the leastsquares sense. The result is typically a nonlinear optimization problem that is computationally expen-
sive and time-consuming to find the solution. In the nonparametric approaches, one is to find the thresholds that separate the gray-level regions of an image in an optimum manner according to some discriminate criteria such as the between-class variance (Otsu, 1979), entropy (Kapur et al., 19851, and cross entropy (Li and Lee, 1993). The nonparametric approaches are computationally efficient and simple to implement, compared to the parametric methods. The between-class variance is the most widely applied method due to its simplicity and consistency of performance on the range of applications considered. However, its performance becomes unreliable under the condition of highly unequal population sixes (Kittler and Illingworth, 1985; Li and Lee, 1993). Lee and Chung (1990) have also found that the between-class variance and entropy methods do not work properly if object sixes become too small. In bilevel thresholding, the existing nonparametric methods are robust and computationally fast for
0167-8655/95/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSDIO167-&3655(95)00005-4
654
D.-M. Tsai/Pattem
Recognition Letters I6 (1995) 653-666
time-critical applications. However, the computational complexity of those methods is exponentially increased and the selected thresholds generally become less credible as the number of classes to be separated increases. In this paper, a fast nonparametric histogrambased algorithm is developed for multi-level thresholding. It assumes each desired class in an image can be represented by an approximately hill-shaped distribution in the gray-level histogram. The fluctuations of an original histogram are smoothed by recursively convoluting the histogram with a Gaussian kernel so that the desired peaks and valleys at varying levels of detail can be obtained. The detected valleys in the smoothed histogram indicate the locations of the thresholds. The proposed algorithm does not need any a priori knowledge about the image other than the gray-level histogram itself. The desired number of classes is the only user-specified parameter to the algorithm. The proposed algorithm is computationally simple and efficient. Actually, the computational time will decrease as the desired number of classes increases. It can detect the proper locations of valleys in a multimodal histogram, even though the population sixes are unequally distributed. In the case that images involve low-contrast intensities, the gray-level distributions may overlap so that the histogram has a single peak. For a unimodal histogram, we cannot find the valley from the smoothed histogram. Instead, we detect the discontinuity of the smoothed histogram curve by the local maximum curvature method. The point of discontinuity on the curve represents the intersection of two overlapping distributions, and is selected as the threshold. This approach is also applied to the image with multimodal histogram, where the resulting smoothed histogram has less peaks than the desired number of classes. The paper is organized as follows: In Section 2, the Gaussian smoothing procedure to detect the peaks and valleys in a multimodal histogram is presented. In Section 3, the maximum curvature method to detect the point of discontinuity in unimodal histograms is described. Section 4 presents the experimental results. The performances of the proposed algorithms and the standard between-class variance and entropy methods are also compared in this section. Section 5 concludes this work.
2. Gaussian smoothing for multimodal histograms The proposed algorithm assumes that a peak shown in the gray-level histogram corresponds to a homogeneous region of the image, and a valley exists between two neighboring peaks. The challenge is to locate the bottoms of the valleys that best separate the classes. Let h,(i) represent the number of pixels in the original image with gray-level i for i=O, 1 , . . . , L, where L is the maximum gray-level. In a gray-level histogram, if h,(i)
> h,( i - 1)
and
h,(i)
>h,(i+
l),
then a peak at gray-level i is detected. Similarly, if h,(j)
1)
and
h,(j)
%(j+
and
h,(j)=0
1)
or h,(j)
and there exist two peaks at gray-levels pr and p2 such that p1 j, then a valley at graylevel j is selected. The definition of peaks and valley above assumes that the peaks and valleys will not be present at gray-levels 0 and L. The fluctuation of the original gray-level histogram may generate many false peaks and valleys. In order to find the proper peaks and valleys at varying levels of thresholding, we use a Gaussian kernel to smooth the histogram. The degree of smoothing is manipulated by the width of the Gaussian kernel and the number of convolutions. Function h,(t), the number of pixels with gray-level t, is convoluted with a one-dimensional Gaussian kernel g(t, a) of width u: g(t,
a) =
&
exp( -rz/2az).
H(t, a), the convolution of h,(t) and the Gaussian kernel, is defined as
=la
h,(u).g(t-u, --co
a) du.
For digital implementation, the digital Gaussian kernel in (Rattarangsi and Chin, 1992) with a window size of W = 3 is used to generate smoothing functions at various values of (T, and it is given by g( - 1) = 0.2261, g(0) = 0.5478, g( 1) = 0.2261
D.-M. Tsai/Pattern
Recognition Letters 16 (1995) 653466
and
for t=[W/2],
Cd4 = 1. u
The discrete convolution of h,(t) and the digital Gaussian kernel g(u) is defined as [W/21 w,
w> =
c u=
hJ(t+
u) *g(u)
-[W/21
where [a] represents the integer truncation operation. For a digital Gaussian smoothing function with higher value of u, the above W = 3 kernel is used in a repeated convolution process. A 2(n + 1) + 1 digital smoothing function is created by repeating the self-convolution process n times. For instance, a W = 5 kernel is obtained by convoluting the W = 3 Gaussian kernel with itself once, and it is given by g( -2)
= g(2) = 0.0511,
g( - 1) = g(1) = 0.2477,
g(0) = 0.4024.
The larger the window size of the Gaussian kernel, the flatter the smoothed curve. The smoothing resulting from a smaller number of convolutions with a large window size kernel can be achieved by a large number of convolutions with a small window size kernel. A Gaussian kernel with large window size may oversmooth the histogram and skip desired peaks. In contrast, a small window size kernel increases the number of convolutions and, hence, reduces the computational efficiency to reach the desired degree of smoothing. Based on the discussion above, the detailed Gaussian smoothing procedure for detecting peaks and valleys in the gray-level histogram is described as follows. The Gaussian smoothing algorithm Denote N = the desired number of threshold levels, h,(t) = the number of pixels with gray-level t at the kth iteration of convolution; h,(t) represents the original histogram. Step 1. Select a large Gaussian kernel size W > 3. Let the current number of iterations be k = 1. Step 2. Convolute h,_ ,(t) with the Gaussian kernel g(u) [W/21 h,(t)=
c u= -[W/21
L(t++g(u),
[W/21+1
655
,..., L-[W/2].
For t < [W/2], h,_ I(t) is convoluted with the Gaussian kernel of size W = 2 t + 1. For t > L W/21, hdt) is convoluted with the Gaussian kernel of size W = 2(L - t) + 1. For t = 0 or L, g(0) = 1. Step 3. Detect the number of peaks in h,(t) as described previously. Denote the detected number of peaks by $C (a) If N > ZV,then let k = k + 1 and repeat Step 2. (b) If i
3. Maximum curvature for &modal
histograms
If the smoothed gray-level histogram of an image is unimodal (or the detected number of peaks with a Gaussian kernel size W = 3 is less than the desired number of threshold levels), the Gaussian smoothing algorithm described in the previous section cannot find the valley to separate the classes. A unimodal histogram is typically obtained when the image consists of low-contrast intensities between distinct classes. Therefore, the distributions may overlap so that the smoothed histogram has a single peak. The distribution of each class will be a smooth curve after the convolutions with the Gaussian kernel. However, the smoothed gray-level histogram may be discontinuous at the point of intersection of
D.-M. Tsai/Pattern Recognition Letters 16 (1995) 653-666
656
is not one of the gray-levels smaller than R and larger than L -R. For a smoothed unimodal histogram, the threshold T is selected at the point of maximum curvature, excluding those points in the neighborhood of the peak. Hence, K,=
Gray level
max{K,, tlD,(p), t
Rgt&L-R},
where L&(p) represents the neighboring points of the peak at gray-level p, i.e.,
Fig. 1. Overlapping of hvo gray-level distributions.
two overlapping distributions, as illustrated in Fig. 1. The dotted lines shown in the figure represent the distributions of two distinct classes, and the solid line indicates the corresponding gray-level histogram. The point of discontinuity in the smoothed histogram is selected as the threshold to segment the image. We employ a local maximum curvature method to detect the point of discontinuity for a unimodal histogram (or a multimodal histogram, where the detected number of peaks is less than the desired number of threshold levels). The curvature K at a point t on the curve is defined as the instantaneous rate of change of +, which is the angle of the tangent at point r with the x-axis, with respect to the arc length s, given by Kt=g The discrete version of curvature for a smoothed histogram h,(t) at the kth convolution is defined as
K,=;
’ Ii
for t=R,
For a smoothed multimodal histogram, where the detected number of peaks (A) is less than the desired number of threshold levels (N), the remaining (N - $1 thresholds are selected at the top (iV - A> curvature points separated by a spacing of at least R points. This discontinuity detection procedure can be incorporated in Step 3(b) of the Gaussian smoothing algorithm aforementioned. The maximum curvature method will be carried out only if the Gaussian smoothing algorithm fails to detect the desired number of peaks. Therefore, the curvature measurement for each gray-level is performed on the smoothed histogram that has less detected peaks than the desired number of threshold levels. Given a distribution in the smoothed histogram, the curve to the left of the distribution’s peak is monotone increasing, and the curve to the right of the peak is monotone decreasing. The computation of tangent angle I,+ is, therefore, less affected by noise. Since a large region of support R may oversmooth the points of discontinuity, a small R (say, R Q 3) is generally sufficient to determine the average tangent angle and curvature at each gray-level.
4. Experimental results I+t+j-tit-j19
R-II ,..., L-R.
In the above equations, R specifies the region of support. It is used as a smoothing factor to compute the mean tangent angle 1(1,and the mean curvature K, at gray-level t. It assumes that the threshold value
The performance of the proposed algorithms is compared to two well-known and most widely applied methods, namely the between-class variance (Otsu, 1979) and entropy (Kapur et al., 1985) methods. For the sake of completeness, these two methods are briefly described as follows.
D.-M, Tsai / Pattern Recognition Letters 16 (I 995) 653-666
657
G
/ 109
j 0.63
j
B
1 107
/ 0.30
1 0.894
E
143 I
II
1.13
0.891
0.578
II
(e) G : Gaussian smoothing B : Between-class variance E : Entropy
Cd) Fig. 2. Bilevel duesholding of a PCB: (a) original image, (b) Gaussian smoothing method, (c) between-class method, (e) comparison of statistical data.
variance method, (d) entropy
D.-M. Tsai/Pattem
658
Recognition Letters 16 (1995) 653-666
Denote the probability of gray-level i in an image by Pit and Pi=W)/
i
i=O, 1,2 ,..., L.
&Wt
j=O
Between-class variance method m,=- ’ ci.p,
Wo=Cpi,
for O
The optimal criteria used in bilevel thresholding (T, = L) and trilevel thresholding, respectively, are E,(T,)
=H, +H,,
E,(T,,
T,) =H,+H,
+H2.
Since there is no guarantee that V,, V,, E2, or E3 does not have a local maximum, the exhaustive search is used to find the best threshold values Tl and T,.
WO
wl=Cpi,
m,=Lzi*p,
for T,
Wl
w2=
cpi,
m2= i
ci.p,
for T,
The optimal criterion used in bilevel thresholding is: V,(T,) = w. * w1 - (m, - mo)2. The threshold value Tl is the gray-level that maximixes the between-class variance V, CT2= L in bilevel thresholding). To eliminate the computation of the second-order statistical moments, the optimal criterion used in trilevel thresholding is given by (VISILOG Programmer’s Guide, 1988): V,(T,, T,) = wo . ~1 * ~2 -
(m, - mo)2( m2 - m,)2( m2 - mo)“.
Tl and T2 are the threshold values that maximize V,. Entropy method wo= cp,,
Ho= -ce
. In pi ( WO 1
for 0 Q i < Tl w1=
CP,,
H,=_q.ln
(pi 1 Wl
for Tl < i Q T, w2=cpi, for T,
H,=
-c$*
In 2 (
1
Pe$ormance The images under test are of size 512 X 360 pixels with 256 gray-levels. The algorithms are implemented on a 486-based personal computer using C. All images are taken under natural room lighting without the support of any special light sources. The Gaussian kernel size W used in the experiment is 3 throughout the entire smoothing procedure so that the worst computational time can be observed. The region of support R used to evaluate the curvature is 1. Experimental data of each observed object is based on the average of three image samples. The uniformity measure adopted from Levine and Nazif (1985) is used to evaluate the performance of the methods. It is given by u=1-ff;/oT’,
OgU
where u..~= within-class variance for given threshold values, and UT”= total gray-level variance of the image. A well segmented image will have the uniformity measure close to 1. In bilevel thresholding, the Gaussian smoothing algorithm and the between-class variance generate similar thresholding results based on the experimental results of six distinct images, although the computational time of the Gaussian smoothing algorithm is not as competitive as that of the between-class variance method. Fig. 2 shows an example of a printed circuit board (PCB). The Gaussian smoothing algorithm and the between-class variance method result in a similar score of uniformity measure, as seen in Fig. 2(e). In trilevel thresholding, the distinction of the resulting threshold values between these three methods is significant. The computational time of the Gaussian smoothing algorithm is decreased, compared to
D.-M. Tsai / Pattern Recognition Letters 16 (I 995) 653-666
1
51
101
659
151
201
251
151
201
251
(e)
200 150 100 50 0 UA-
(b)
viethod1’1kesholc 1 Time (sec.)
measure
114, 193
0.36
0.937
97, 196
11.50
0.781
153,201I 123.63
0.867
values G
(4
LJniformity
iaussian smoothing Ietween-class variance E : Entropy
Fig. 3. Trilevel thesholding of a PC3 with small population of white intensity characters: (a) original image, (b) Gaussian smoothing method, (C) between-class variance method, (d) entropy method, (e) original histogram of (a), (f) smoothed histogram with three peaks, (g) comparison of statistical data.
660
D.-M. Tsai / Pattern Recognition Letters 16 (I 995) 653-666
r
300 250 200 150 100 50
0 I-AA1
51
101
151
201
251
51
101
151
201
251
(4 250 200 150
U
1
Method IIreshold values
Time
uniformity
kc.)
measure
G
61, 112
0.33
0.900
B
80, 130
11.17
0.702
E
67,146
95.00
0.358
L (4
(g) G : Gaussian smoothing B : Between-class variance E : Entropy
Fig. 4. Trilevel thresholding of a PC23 with small population of black intensity (pin holes): (a) original image, (b) Gaussian smoothing method, (c) between-class variance method, (d) entropy method, (e) original histogram of (a), (f) smoothed histogram with three peaks, (g, comparison of statistical data.
D.-M. Tsai / Pattern Recognition Letters 16 (1995) 653-666
661
200 r
100 [50. 0 1
(a)
51
101
151
201
251
(d) 150 100 50 0 1
51
1
51
101
151
201
151
201
251
(e)
Co)
150
100 50 0
(c)
101
251
(0 Thresholding level Four-level Five-level
Threshold
Time
Uniformity
values
(see.)
measure
0.77
0.905
0.14
0.921
102,143,201 42,102,142,201
(g) Fig. 5. Four-level and five-level thresholding of a PCB: (a) original image, (b) four-level segmented image, (c) five-level segmented image, (d) original histogram (a), (e) smoothed histogram with four peaks, (f) smoothed histogram with five peaks, (g) experimental data.
D.-M. Tsai/Pattern Recognition Letters I6 (1995) 653-666
662
(e)
Cd) Fig. 6. An aluminum cap image with unimodal histogram: (a) original image, (b), (c) and (d) bilevel tbresholding results of the Gaussian smmtig, between*lass variance and entropy methods, respectively, (e), (f) and
D.-M. Tsai/Pattern
663
Recognition Letters 16 (1995) 653466
that for bilevel thresholding, since less convolutions will be required for a larger number of threshold levels. In contrast, the computational times of the between-class variance and entropy methods are dramatically increased. The entropy method is extremely time-consuming since it involves the evaluation of the logarithm function. In an experiment of six distinct images, the average computational times using the Gaussian smoothing, between-class variance and entropy methods are 0.4, 11.1 and 107.1 seconds, respectively. Figs. 3(b)-3(d) show the trilevel thresholding results of a PCB using the three methods. Figs. 3(e) and 3(f) illustrate the original histogram of the image in Fig. 3(a) and the resulting smoothed histogram with three peaks, respectively. Note that the characters (the white parts in Fig. 3(a)) printed on the PCB possess only a small portion of pixels with respect to the entire image area. The corresponding peaks in the histogram are extremely unequal in height as seen in Fig. 3(e). However, the Gaussian smoothing algorithm can still locate the valley and segment the white characters from the background properly. Figs. 4(bMd) show another trilevel thresholding result of a PCB using the three methods. The electric circuit and the background are segmented in a scattering way using the between-class variance method. The electric circuit is missing using the entropy method. Figs. 4(e) and 4(f) illustrate the original histogram of the image in Fig. 4(a) and the resulting smoothed histogram, respectively. Note that the population sizes of the pin holes (the black parts in Fig. 4(a)), electric circuit and background are distributed in a highly unequal way. However, the Gaussian smoothing algorithm preserves the pin holes distribution in the smoothed histogram and properly locates the valley to separate the pin holes from the background. Figs. 5(b) and 5(c) show the four-level and fivelevel thresholding results for a complex PCB, respectively. Figs. 5(d)-5(f) illustrate the original histogram, the smoothed histograms with four peaks and five peaks, respectively. In four-level thresholding, the pin holes and the background are segmented into the same class. The white characters (1 and Ul shown in Fig. 5(a)) which possess only an extremely small population size are properly detected from the four-peak histogram. The pin holes are separated
from the background using the five-peak histogram, even though they involve only a very small population size. As seen in Figs. 3(g), 4(g) and 5(g), the Gaussian smoothing algorithm gives high scores of uniformity measures around 0.9. Under the Gaussian smoothing, a small peak can persist if it is isolated while a significant peak will be merged with its neighboring peak if these two peaks are very close to each other. This situation is observed in Figs. 5(e) and 5(f). Under four-level thresholding, the small peak corresponding to the white characters (1 and Ul shown in Fig. 5(a)) persists and the distributions of the pin holes and the background are merged into the same class. If separating the pin holes from the background is desired, we may conclude that four-level thresholding is not appropriate for the image shown in Fig. 5(a). Therefore, the ideal number of peaks for this image should be selected to be five so that the white characters, the pin holes and the background can be properly segmented. Alternatively, we may consider the white characters as noise in the image due to their small population, and neglect the existence of the isolated, small peak in the histogram. Therefore, we can identify three valleys (for four-level thresholding) at gray-levels 40, 102, and 142 from the histogram shown in Fig. 5(f). This implementation needs a priori knowledge about the distribution of noise. We have also selected two low-contrast texture images to obtain unimodal histograms so that the performance of the maximum curvature method can be evaluated. Fig. 6(a) shows the picture of an aluminum cap. In the image, the mark pressed on the cap has identical intensity with the surface of the cap except its shaded edges. Fig. 7 illustrates the corresponding unimodal histogram after Gaussian smoothing. The largest and the second largest curvatures are located at the gray-level 89 and 99, respectively. The
61
76
91
106
121
Fig. 7. The unimodal histogram of the image shown in Fig. 6(a).
664
D.-M. Tsai/Pattem
Recognition Letters 16 (1995) 653-666
Cd) Fig. 8. A leather image with unimodal histogram:(a) original image, (b), (c) and (d) bilevel thresholdingresults of the Gaussian smoothing, between-class variance and entropy methods, respectively, (e), (f) and (g> trilevel tbresholding results of the three methods, respectively.
D.-M. Tsai/Pattem
665
Recognition Letters I6 (1995) 653-666
Table 1 Comparison of statistical data for the image in Fig. 6(a) Trilevel thresholding
Bilevel thresholding Method Gaussian smoothing Between-class variance Entropy
Threshold value 89 93 86
Uniformity measure
Threshold value
Uniformity measure
0.474 0.549 0.372
89,99 81,100 86,103
0.688 0.654 0.613
Table 2 Comparison of statistical data for the imaee in Fin. 8(a) Trilevel thresholding
Bilevel thresholding Method Gaussian smoothing Between-class variance Entropy
Threshold value 70 74 100
Uniformity measure
Threshold value
Uniformity measure
0.426 0.621 0.136
71,85 68,82 87,100
0.697 0.626 0.369
bilevel thresholding results of the Gaussian smoothing with thresholding value 89, between-class variance and entropy methods are presented in Figs. 6(b), 6(c) and 6(d), respectively. The trilevel thresholding results of the three methods are shown in Figs. 6(e), 6(f) and 6(g). Table 1 summarizes the selected threshold values and corresponding uniformity measures of the three methods. The second low-contrast texture image is a polo figure pressed on a leather wallet as shown in Fig. 8(a). Figs. @b&8(d) and Figs. S(e)-S(g) present the bilevel and trilevel thresholding results of the leather image using the three methods, respectively. Table 2 presents the statistical data of the three methods. The threshold values selected by the maximum curvature method are similar to those determined by the behveen-class variance method. The segmentation results of these two low-contrast images are very sensitive to the threshold values. A difference of four gray-levels may result in 13% - 30% deviation of the uniformity measure between the maximum curvature and between-class variance methods, as seen in Tables 1 and 2 for the bilevel thresholding case.
valleys in a multimodal histogram, and uses local maximum curvature to detect points of discontinuity in a unimodal histogram. The proposed algorithms require only a user-specified parameter, namely the desired number of classes. The threshold values selected by the proposed algorithms are close to those determined by the between-class variance method. The computational time of the Gaussian smoothing algorithm will decrease as the desired number of classes increases. It becomes more attractive when the desired levels of thresholding is greater than two. This makes the proposed methods appropriate for time-critical applications, where the images involve multimodal or unimodal histograms.
References Glasbey, C.A. (1993). An analysis of histogram-based thresholding algorithms. CVGIP: Graphical Models and Image Processing 55, 532-537.
Kapur, J., P. Sahoo and A. Wong (1985). A new method for gray-level picture thresholding using the entropy of the histogram. Computer Vision, Graphics, and Image Processing 29, 273-285.
5. Conclusions A simple and efficient histogram-based approach has been presented for multi-level thresholding. It uses Gaussian kernel smoothing to detect peaks and
Kittler, J. and J. Illingworth (1985). On threshold selection using clustering criteria. IEEE Trans. Cyst. Man Cyberner. 15, 652655. Lee, S.U. and S.Y. Chung (1990). A comparative performance
study of several global thresholding techniques for segmentation. Computer Vision, Graphics, and Image Processing 52, 171-190.
666
D.-M. Tsai/Pattern Recognition Letters 16 (1995) 653-666
Levine, M.D. and A.M. Nazif (1985). Dynamic measurement of computer generated image segmentations. IEEE Trans. Puttern Anal. Mach. Intell. 7, 155-164. Li, C.H. and C.K Lee (1993). Minimum cross entropy thresholding. Pattern Recognition 26, 617-625. Otsu, N. (1979). A threshold selection method for gray-level histogram. IEEE Trans. Syst. Man Cybernet. 9, 62-66. Rattarangsi, A. and R.T. Chin (1992). Scale-based detection of comers of planar curves. IEEE Trans. Pattern Anal. Mach. Intell. 14, 430-449.
Sahoo, P., S. Soltani and A. Wong (1988). A survey of thresholding techniques. Computer Vision, Graphics, and Image Processing 41, 233-260. Synder, W. (1990). Optimal thresholding - a new approach. Pattern Recognition Len. 11, 803-810. VISILOG Programmer’s Guide (1988). Noesis S.A.R.L., France. Weszka, J. and A. Rosenfeld (1979). Histogram modifications for threshold selection. IEEE Trans. Syst. Mun Cybernet. 9, 3852.