001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 *  image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause 
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 
006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009
010package imagingbook.pub.threshold.global;
011
012import imagingbook.pub.threshold.global.GlobalThresholder;
013
014/**
015 * Thresholder as described in N. Otsu, "A threshold selection method from gray-level histograms",
016 * IEEE Transactions on Systems, Man, and Cybernetics 9(1), 62-66 (1979).
017 */
018public class OtsuThresholder extends GlobalThresholder {
019        
020        private int[] h = null;
021        private double[] M0 = null;             // table of background means
022        private double[] M1 = null;             // table of foreground means
023        private int N = 0;                              // number of image pixels
024        
025        public OtsuThresholder() {
026                super();
027        }
028        
029        public int getThreshold(int[] hist) {
030                h = hist;
031                int K = h.length;
032                N = makeMeanTables(h);
033
034                double sigma2Bmax = 0;
035                int qMax = -1;
036                int n0 = 0;
037                
038                // examine all possible threshold values q:
039                for (int q = 0; q <= K-2; q++) {
040                        n0 = n0 +  h[q]; 
041                        int n1 = N - n0;
042                        if (n0 > 0 && n1 > 0) {
043                                double meanDiff = M0[q] - M1[q];
044                                double sigma2B =  meanDiff * meanDiff * n0 * n1; // (1/N^2) has been omitted
045                                if (sigma2B > sigma2Bmax) {
046                                        sigma2Bmax = sigma2B;
047                                        qMax = q;
048                                }
049                        }
050                }
051                
052                return qMax;
053        }
054        
055        private int makeMeanTables(int[] h) {
056                int K = h.length;
057                M0 = new double[K];
058                M1 = new double[K];
059                int n0 = 0;
060                long s0 = 0;
061                for (int q = 0; q < K; q++) {
062                        n0 = n0 + h[q];
063                        s0 = s0 + q * h[q];
064                        M0[q] = (n0 > 0) ? ((double) s0)/n0 : -1;
065                }
066                
067                int N = n0;
068                
069                int n1 = 0;
070                long s1 = 0;
071                M1[K-1] = 0;
072                for (int q = h.length-2; q >= 0; q--) {
073                        n1 = n1 + h[q+1];
074                        s1 = s1 + (q+1) * h[q+1];
075                        M1[q] = (n1 > 0) ? ((double) s1)/n1 : -1;
076                }
077                
078                return N;
079        }
080}