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-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009
010package imagingbook.common.threshold.global;
011
012import static imagingbook.common.math.Arithmetic.sqr;
013
014/**
015 * <p>
016 * This is an implementation of the global thresholder proposed by Otsu [1]. See Sec. 9.1.4 (Alg. 9.4 and 9.3) in [2]
017 * for a detailed description.
018 * </p>
019 * <p>
020 * [1] N. Otsu, "A threshold selection method from gray-level histograms", IEEE Transactions on Systems, Man, and
021 * Cybernetics 9(1), 62-66 (1979). <br> [2] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic
022 * Introduction</em>, 3rd ed, Springer (2022).
023 * </p>
024 *
025 * @author WB
026 * @version 2022/08/21
027 */
028public class OtsuThresholder implements GlobalThresholder {
029
030        private double[] M0;            // table of background means
031        private double[] M1;            // table of foreground means
032        private int N;                          // number of image pixels
033        
034        public OtsuThresholder() {
035                super();
036        }
037        
038        @Override
039        public float getThreshold(int[] h) {
040                int K = h.length;
041                makeMeanTables(h);
042
043                double sigma2Bmax = 0;
044                float qMax = NoThreshold;
045                int n0 = 0;
046                
047                // examine all possible threshold values q:
048                for (int q = 0; q <= K-2; q++) {
049                        n0 = n0 +  h[q];                        // # of background pixels for q
050                        int n1 = N - n0;                        // # of foreground pixels for q
051                        if (n0 > 0 && n1 > 0) {         // both sets must be non-empty
052                                double sigma2B =  sqr(M0[q] - M1[q]) * n0 * n1;         // factor (1/N^2) omitted
053                                if (sigma2B > sigma2Bmax) {
054                                        sigma2Bmax = sigma2B;
055                                        qMax = q;
056                                }
057                        }
058                }
059                
060                return qMax;
061        }
062        
063        private void makeMeanTables(int[] h) {
064                int K = h.length;
065                this.M0 = new double[K];
066                this.M1 = new double[K];
067                int n0 = 0;
068                long s0 = 0;
069                for (int q = 0; q < K; q++) {
070                        n0 = n0 + h[q];
071                        s0 = s0 + q * h[q];
072                        M0[q] = (n0 > 0) ? ((double) s0) / n0 : -1;
073                }
074                
075                this.N = n0;
076                
077                int n1 = 0;
078                long s1 = 0;
079                M1[K-1] = 0;
080                for (int q = h.length-2; q >= 0; q--) {
081                        n1 = n1 + h[q+1];
082                        s1 = s1 + (q+1) * h[q+1];
083                        M1[q] = (n1 > 0) ? ((double) s1)/n1 : -1;
084                }
085        }
086}