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
012/**
013 * <p>
014 * This is an implementation of the global "minimum error" thresholder proposed by Kittler and Illingworth in [1]. See
015 * Sec. 9.1.6 (Alg. 9.6) of [2] for a detailed description.
016 * </p>
017 * <p>
018 * [1] J. Illingworth and J. Kittler. "Minimum error thresholding". Pattern Recognition 19(1), 41–47 (1986). <br> [2] W.
019 * Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer (2022).
020 * </p>
021 *
022 * @author WB
023 * @version 2022/08/01
024 */
025public class MinErrorThresholder implements GlobalThresholder {
026        
027        static final double EPSILON = 1E-12;
028        
029        private double[] S2_0, S2_1;
030        private int N;
031        
032        public MinErrorThresholder() {
033                super();
034        }
035        
036        @Override
037        public float getThreshold(int[] h) {
038                int K = h.length;
039                makeSigmaTables(h);     // set up S2_0, S2_1, N
040
041                int n0 = 0, n1;
042                float qMin = NoThreshold;
043                double eMin = Double.POSITIVE_INFINITY;
044                for (int q = 0; q <= K-2; q++) {
045                        n0 = n0 + h[q];
046                        n1 = N - n0;
047                        if (n0 > 0 && n1 > 0) { 
048                                double P0 = (double) n0 / N;    // could use n0, n1 instead
049                                double P1 = (double) n1 / N;
050                                double e = // 1.0 + 
051                                        P0 * Math.log(S2_0[q]) + P1 * Math.log(S2_1[q])
052                                        - 2 * (P0 * Math.log(P0) + P1 * Math.log(P1));
053                                if (e < eMin) {         // minimize e;
054                                        eMin = e;
055                                        qMin = q;
056                                }
057                        }
058                }
059                return qMin; 
060        }
061        
062        private void makeSigmaTables(int[] h) {
063                int K = h.length;
064                final double unitVar = 1d/12;   // variance of a uniform distribution in unit interval
065                S2_0 = new double[K];
066                S2_1 = new double[K];
067                
068                int n0 = 0;
069                long A0 = 0;
070                long B0 = 0; 
071                for (int q = 0; q < K; q++) {
072                        long ql = q;    // need a long type to avoid overflows
073                        n0 = n0 + h[q];
074                        A0 = A0 + h[q] * ql;
075                        B0 = B0 + h[q] * ql*ql;
076                        S2_0[q] = (n0 > 0) ? 
077                                        unitVar + (B0 - (double)A0*A0/n0)/n0 : 0;       
078                }
079                
080                N = n0;
081                
082                int n1 = 0;
083                long A1 = 0;
084                long B1 = 0; 
085                S2_1[K-1] = 0;
086                for (int q = K-2; q >= 0; q--) {
087                        long qp1 = q+1;
088                        n1 = n1 + h[q+1];
089                        A1 = A1 + h[q+1] * qp1;
090                        B1 = B1 + h[q+1] * qp1*qp1;
091                        S2_1[q] = (n1 > 0) ? 
092                                        unitVar + (B1 - (double)A1*A1/n1)/n1 : 0;
093                }
094        }
095        
096}