Wilkinson, Michael H.F.

English

A multi-scale, moving-window method for local thresholding based on Robust Automatic Threshold Selection (RATS) is developed. Using a model for the noise response of the optimal edge detector in this context, the reliability of thresholds computed at different scales is determined. The threshold computed at the smallest scale at which the reliability is suffcient is used. The performance on 2-D images is evaluated on synthetic an natural images in the presence of varying background and noise. Results show the method deals better with these problems than earlier versions of RATS at most noise levels

ARTS repository - University of Groningen

   University of GroningenGaussian-weighted moving-window robust automatic threshold selectionWilkinson, Michael H.F.Published in:COMPUTER ANALYSIS OF IMAGES AND PATTERNS, PROCEEDINGSIMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite fromit. Please check the document version below.Document VersionPublisher's PDF, also known as Version of recordPublication date:2003Link to publication in University of Groningen/UMCG research databaseCitation for published version (APA):Wilkinson, M. H. F. (2003). Gaussian-weighted moving-window robust automatic threshold selection. In N.Petkov, & MA. Westenberg (Eds.), COMPUTER ANALYSIS OF IMAGES AND PATTERNS,PROCEEDINGS (pp. 369-376). (LECTURE NOTES IN COMPUTER SCIENCE; Vol. 2756). University ofGroningen, Johann Bernoulli Institute for Mathematics and Computer Science.CopyrightOther than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of theauthor(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-amendment.Take-down policyIf you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediatelyand investigate your claim.Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons thenumber of authors shown on this cover page is limited to 10 maximum.Download date: 06-06-2022Gaussian-Weighted Moving-Window RobustAutomatic Threshold SelectionMichael H. F. WilkinsonInstitute for Mathematics and Computing Science,University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands,michael@cs.rug.nlhttp://www.cs.rug.nl/~michael/Abstract. A multi-scale, moving-window method for local thresholdingbased on Robust Automatic Threshold Selection (RATS) is developed.Using a model for the noise response of the optimal edge detector in thiscontext, the reliability of thresholds computed at dierent scales is de-termined. The threshold computed at the smallest scale at which the re-liability is sucient is used. The performance on 2-D images is evaluatedon synthetic an natural images in the presence of varying backgroundand noise. Results show the method deals better with these problemsthan earlier versions of RATS at most noise levels.1 IntroductionIn all applications of thresholding, correct selection of the threshold is the keyissue, and many methods for automatic selection of optimal thresholds havebeen published[14, 8]. Ideally, thresholds should be computed locally, adaptingto the local image statistics, to deal properly with a locally varying background,or variations in the grey level of objects, both of which may occur in a singleimage [8]. An example is shown in Fig. 1, in which a local variations in objectintensity are compensated through local threshold selection.In this paper I will extend a local thresholding method, called Robust Au-tomatic Threshold Selection (RATS) [1]. A new, moving-window version of thealgorithm using Gaussian convolution will be developed. This versions will beextended to a multi-scale method, in which the smallest scale at which reliablethresholds can be computed is nally used. To do this, the eect of Gaussiannoise on the computed threshold is derived. Finally, the method is evaluated rston synthetic images with varying noise levels, and later on natural ones.The applications for the method include 2-D microscopic images of microor-ganisms, and 3-D angiograms.2 Robust Automatic Threshold SelectionRATS [1] is a method for bilevel thresholding of grey scale images, which hasbeen applied to images of bacteria [5, 6]. It is based on a simple image statistic,N. Petkov and M.A. Westenberg (Eds.): CAIP 2003, LNCS 2756, pp. 369-376, 2003c© Springer-Verlag Berlin Heidelberg 2003370(a) (b) (c) (d)Fig. 1. Local thresholding: (a) uorescence image of bacteria; (b) global thresholdingusing original RATS algorithm without noise correction; (c) global thresholding usingRATS algorithm with square Sobel gradient lter and noise correction; (d) locallythresholded result using RATS with the quad-tree approach from [5].which is the average of grey levels weighted by the edge strength at each point.Kittler et al. [1] show that the optimal threshold in a noise-free image T isT =∑e(x, y)p(x, y)∑e(x, y), (1)in which p(x, y) is the grey level at (x, y) and the edge strength e is given bye(x, y) = max(|gx(x, y)|, |gy(x, y))|, (2)withgx(x, y) = p(x−1, y)−p(x+1, y) and gy(x, y) = p(x, y−1)−p(x, y+1). (3)Initially the optimality of T was proved only for gradient operator in (2), andfor straight edges. It has since been shown that any edge detector with an evenresponse to a step edge at the origin will yield the same optimal result [7]. Inparticular, the gradient detectorg2(x, y) = g2x(x, y) + g2y(x, y) (4)shows no curvature bias, is rotation invariant, and has reduced noise bias. How-ever, the reduced noise bias comes at the expense of increased variance, whichcan be countered by using Sobel lter kernels to compute x and y derivatives[7]. The method readily extends to 3-D images, by replacing the edge detectorsto their 3-D counterparts.In the presence of noise T is biased towards the most common category inthe image (usually background) [1]. This noise bias is counteracted by using athreshold on the edge strength below which the pixels receive zero weight. Thestatistic now becomesTλ =∑wλ(x, y)p(x, y)∑w(x, y), (5)371withwλ(x, y) ={g2(x, y) if g2(x, y) > λ2η20 otherwise,(6)in which η is the standard deviation of the image noise, and λ is an adjustable pa-rameter, which depends on the actual edge strength used. For the edge strengthdened in (2) it was shown empirically that λ = 5 is a good choice for Gaussiannoise [7].2.1 Local Application of RATSRATS lends itself well to local application [1] for two reasons: (i) the statistic in(5) is robust against noise, and (ii) it is easy to check whether a region containsan edge by checking whether the denominator in (5) is above some threshold[1, 7]. Ideally, we want to compute the threshold in an isotropic surroundings ofeach pixel. This can be done using a moving-window version of RATS, whichcan be written as the ratio of two convolutionsTh(x, y) =(Πh ∗ (w · p))(x, y)(Πh ∗ w)(x, y), (7)in which ∗ denotes convolution and Πh(x, y) is given byΠh(x) ={1 if |x| ≤ h and |y| ≤ h,0 otherwise.(8)One problem with (7) is that Th is undened for all pixels where (Πh∗w)(x, y) =0. Thus wherever the distance between edges is greater than the width of thewindow, no threshold is computed. Besides, the square convolution kernels arenot isotropic. However, the convolution formalism allows generalization of thealgorithm to other convolution kernels, e.g. Gaussian. Gaussian kernels are sep-arable, and have innite impulse response (IIR), and so will contribute over theentire image. Besides, can be computed quickly using a recursive implementationwhich has an IIR [9]. We arrive atTσ(x, y) =(Gσ ∗ (w · p))(x, y)(Gσ ∗ w)(x, y), (9)with Gσ a Gaussian with zero mean and standard deviation of σ. A furtheradvantage (9) over (7) is that edges close to the current pixel are given higherweights than distant edge. However, despite their innite impulse response, Gσfalls o so rapidly that at large distances from edges thresholds may becomeunreliable, because a few remaining nearby noise edges may outweigh distanttrue edges. To counter this, multi-scale approach can be used, by computingTσ with σ = σ0, 2σ0, 4σ0, . . . , σmax, and using the lowest σ for which Tσ can becomputed reliably.3722.2 Selecting the Correct ScaleTo select the correct scale we need to understand how the noise inuences thestatistic Tσ. Let us assume that the noise in an image is Gaussian with zeromean and standard deviation η. Its distribution is simplyp(x)dx =1η√2πe−x2/2η2 (10)This means that the probability distribution pgx of gx (or gy or gz) ispgx(x)dx =12η√πe−x2/4η2 . (11)The probability distribution p1 of g2x isp1(x)dx = pgx(√x)d√x =12η√2πxe−x/4η2. (12)In 2-D, the probability distribution p2 of g2 in (4) isp2(x)dx = (p1 ∗ p1)(x)dx =14η2e−x/4η2, (13)The 3-D counterpart of g2 has a probability distribution p3 given byp3(x)dx = (p2 ∗ p1)(x)dx =√x4η3√2πe−x/4η2. (14)To select the lowest scale at which the denominator in (9) becomes signif-icantly dierent from the value expected from noise, we need both the meanvalue 〈wλ〉 and the variance σ2w of wλ. Using (6) and (13) the distribution pw ofwλ ispw(x) ={δ(x)∫ λ2η20p2(x′)dx′ if x < λ2η2,p2(x) otherwise.(15)It can be seen from (13) that the probability pλ that x > λ2η2 is e−λ2/4. There-fore, the expected value 〈wλ〉 is〈wλ〉 =∫ ∞λ2η2x4η2e−x/4η2dx = e−λ24(1 +λ24)4η2 (16)and the variance σ2w isσ2w = e−λ24((1− e−λ24 )(1 +λ24)2+ 1)16η4 ≈ e−λ24(1 +λ24)216η4 (17)in the 2-D case. Thus, the standard deviation σw is approximatelyσw ≈ e−λ28(1 +λ24)4η2. (18)373For λ ≥ 5 the approximation is accurate to within 6%. The mean noise responsein the Gaussian-convolved edge image is just 〈wλ〉, whereas the variance of thenoise response σ2Gw isσ2Gw = σ2w∫ ∞−∞∫ ∞−∞G2σi(x, y)dxdy =14πσ2iσ2w. (19)We can select the lowest scale σi for which the denominator of (9) is larger thana threshold TGw = 〈wλ〉+ 3σGw, or(Gσi ∗ w)(x, y) ≥ TGw = e−λ28(1 +λ24)(e−λ28 +32σi√π)4η2 (20)Note that as σi →∞, the variance σ2Gw → 0, and so TGw → 〈wλ〉. When usingthe Sobel kernels in our initial edge detector, the only thing that needs to bechanged in this calculation is to replace η with√5η/4 [7].3 AlgorithmLet array p contain the original image, array w store the weights, array wp theproduct wλp, array T the threshold, and a boolean array q the binary outputimage. All arrays are of the same size as the original image. Let the value Invdenote an invalid threshold (e.g., some value > maxx,y(p(x, y))). Finally, we haveTλ to store the global threshold according to (5).The multi-scale, Gaussian-weighted, moving-window RATS algorithm is sum-marized in Fig. 2. The only input the algorithm needs is the original image, thedesired value of λ, and the image noise η. After computing of wλ and storingin in w, we compute the product image wλp, and the global threshold Tλ. Weinitialize all elements of T to Inv, and convolve both w and wp with Gσ0 atthe lowest scale. After this initial phase, we loop through all scales but the last.During each loop we rst compute TGw at that scale, and then compute Tσifor all pixels (x, y) which have not yet been assigned a valid threshold (i.e.,T (x, y) = Inv), and for which w(x, y) ≥ TGw. We then compute Gσi+1 ∗w fromGσi ∗ w by convolving Gσi ∗ w with G√3σi , and likewise for Gσi+1 ∗ wp. At thelargest scale, a similar operation is performed, but here T (x, y) is set to theglobal threshold Tλ if no threshold can be computed at that scale. Finally, foreach pixel q(x, y) is set to true if p(x, y) > T (x, y) and to false otherwise.4 ResultsThe algorithm was implemented, using the Sobel gradient, λ = 7, and four scalesranging from σ = 2 to σ = 16, and tested on synthetic 2-D images of 256× 256pixels containing objects of dierent sizes and dierent or constant contrastwith respect to the local background. The constant object-intensity images (seeFig. 3(b)) served as reference segmentation for themselves and the correspond-ing variable object-intensity images as in Fig. 3(a).A background slope running3741. For all pixels (x, y) compute wλ(x, y) from input image p and store in w.2. Compute product image wλp and store in wp.3. Compute Tλ from wp and w.4. Set all values in T to Inv.5. Convolve w and wp with Gσ0 and store in w and wp respectively.6. For all scales i = 0, 1, . . . ,max−1 do7. Compute TGw for this scale8. For all pixels (x, y) with T (x, y) 6= Inv9. if w(x, y) ≥ TGw then10. T (x, y)← wp(x, y)/w(x, y).11. Convolve w and wp with G√3σi and store in w and wp respectively.12. Compute TGw for σmax13. For all pixels (x, y) with T (x, y) 6= Inv14. if w(x, y) ≥ TGw then15. T (x, y)← wp(x, y)/w(x, y).16. else17. T (x, y)← Tλ18. For all pixels (x, y)19. q(x, y)← p(x, y) > T (x, y)Fig. 2. The multi-scale, Gaussian-weighted, moving-window RATS algorithm. Notethat ← denotes assignment.from 0 at the left and height hs at the right, and Gaussian noise with standarddeviation η were added. The fraction of misclassied pixels was computed as afunction of η and hs. The results are shown in Fig. 3(c)-(h). Using all four scalesthe performance is generally good up to an η of about 8 (corresponding to anS/N-ratio of about 8 for the faintest objects) for objects of varying intensity(Fig. 3(c), (d) and (f)). Fig. 3(d) shows how omission of the lowest scale resultsin segmentation errors where faint objects lie close to bright. Fig. 3(f) shows acomparison of the new method with global thresholding by Otsu's method [2].The latter performs badly at all noise levels, because it classies fainter fore-ground objects as background. For the former method the error fraction risessharply beyond η = 8 to about 7% at η = 32. An example is shown in Fig. 3(e).The slope hs has no impact on segmentation quality (not shown). Segmentationof Fig. 3(b) is easier, the errors never exceeding 0.5 % (not shown). Using justa single scale σ results in rather poor segmentation, unless the smallest scalesare used Fig. 3(g). Apparently, the lowest values of σ have a high impact onthe quality of segmentation. Finally, Fig 3(h) shows a comparison between theearlier quad-tree method and the new algorithm for hs = 16. Only at η = 16does the old method perform better, suggesting that the way noise is dealt withusing (20) could be improved, perhaps by giving less weight to the lowest scale.The results in Fig. 3(g) also suggest that using σ = 4 at η = 16 is better thanthe multi-scale approach.375(a) (b) (c) (d) (e)0 5 10 15 20 25 30 3500.010.020.030.040.050.060.070.08Noise η in grey levelsFraction of misclassified pixelsMS RATSGlobal Otsu0 5 10 15 20 25 30 3500.010.020.030.040.050.060.070.08Noise η in grey levelsFraction of missclassified pixelsσ = 2σ = 4σ = 8σ = 160 5 10 15 20 25 30 3500.010.020.030.040.050.060.070.08Noise η in grey levelsFraction of missclassified pixelsnew method4 levels5 levels6 levels(f) (g) (h)Fig. 3. Segmentation of synthetic images of ellipses: (a) synthetic image of ellipses ofdiering intensities; (b) same as (a) but with constant intensities; (c) segmentationresult of (a) with noise η = 1 and slope hs = 0 added; (d) as (c) but without usingscale σ = 2; (e) as (c) but η = 32 and hs = 32; (f) fraction of correctly classied pixelsof (a) as a function of η, for multi-scale RATS with hs = 0, using 4 scales plus globalthreshold Tλ, compared to global thresholding according to Otsu [2]; (g) same as (f)but using just a single scale σ and Tλ; (h) same as (f) but comparing the new methodto the quad-tree method [5] for dierent numbers of levels in the quad-tree.The method was also tested on images of bacteria, with σ0 ranging from 2to 8. As can be seen in Fig. 4, when three scales are used, the method detects afaint object skipped by the quad-tree approach shown in Fig. 1. If σ0 = 2, themethod detects parts of the diraction halos around the brighter objects.5 DiscussionA new, multi-scale version of RATS has been developed which can adapt well tovariations in both background and object intensity. A framework to select theappropriate scale has been developed. The experiments show that the methodworks with a modest number of scales, and with therefore a modest computa-tional cost. A 512× 512 image takes just 0.16 s if four scales are used, whereas a2483× 3508 image takes 6.03 s on a Pentium 4 at 1.9 GHz. At a single scale thetimings are 0.06 s and 1.83 s, respectively. Extensions to 3-D should be straight-forward, provided a 3-D equivalent of (20) is derived. The selection of the lowestscale may need more work. The experiment with synthetic images yielded thebest results with σ0 = 2, whereas the experiment using images of bacteria yieldedthe best results with σ0 = 4. It might well be that the Gaussian assumptionsused in (20) do not hold for small σ. Solving this problem requires analytical376(a) (b) (c) (d)Fig. 4. Segmentation of image of bacteria from Fig. 1(a) with σmax = 16, λ = 5 andη = 2.3: (a) σ0 = 8; (b) σ0 = 4; (c) σ0 = 2; (d) contrast stretched original showingfaint object and diraction halo around brighter objects.solutions or numerical approximations of the distribution of the denominator of(9). This will be studied in future work.One drawback of this and many other implementations of RATS is that weneed an estimate of the image noise. Ideally, we would like to be able to determinethe noise from the image itself, rather than rely on external calibration data,which might be absent. If we assume that the lowest gradient pixels (or voxels)are attributable to noise only, we could estimate η by tting the histogram ofthe edge ltered image at these values to the distribution in (13) in 2-D or (14)in 3-D. Work is in progress for just such an extension.References1. J. Kittler, J. Illingworth, and J. Föglein. Threshold selection based on a simpleimage statistic. Comp. Vision Graph. Image Proc. , 30:125147, 1985.2. N. Otsu. A threshold selection method from gray-level histogram. IEEE Trans.Syst. Man Cybern., 9:6266, 1979.3. P. K. Sahoo, S. Soltani, A. K. C. Wong, and Y. C. Chen. A survey of thresholdingtechniques. Comp. Vision Graph. Image Proc. , 41:233260, 1988.4. Ø. D. Trier and A. K. Jain. Goal-directed evalution of binarization methods. IEEETrans. Image Proc., 17(12):11911201, 1995.5. M. H. F. Wilkinson. Rapid automatic segmentation of uorescent and phase-contrast images of bacteria. In J. Slavik, editor, Fluorescence Microscopy and Flu-orescent Probes, pages 261266. Plenum Press, New York, 1996.6. M. H. F. Wilkinson. Automated and manual segmentation techniques in imageanalysis of microbes. In M. H. F. Wilkinson and F. Schut, editors, Digital ImageAnalysis of Microbes , pages 135171. John Wiley and Sons, Ltd, Chichester, UK,1998.7. M. H. F. Wilkinson. Optimizing edge detectors for robust automatic thresholdselection: coping with edge curvature and noise. Graph. Mod. Image Proc. , 60:385401, 1998.8. Y. Yang and H. Yan. An adaptive logical method for binarization of degradeddocument images. Pattern Recognition, 33:787807, 2000.9. I. T. Young and L. J. van Vliet. Recursive implementation of the Gaussian lter.Signal Processing, 44:139151, 1995.

Gaussian-weighted moving-window robust automatic threshold selection

https://pure.rug.nl/ws/files/3022787/2003LNCSWilkinson.pdf

Gaussian-weighted moving-window robust automatic threshold selection

Abstract

Similar works

Full text

Available Versions

ARTS repository - University of Groningen