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 * Minimum Error thresholder after Kittler and Illingworth (1986).
016 */
017public class MinErrorThresholder extends GlobalThresholder {
018        
019        static final double EPSILON = 1E-12;
020        
021        private double[] S2_0, S2_1;
022        private int N;
023        
024        public MinErrorThresholder() {
025                super();
026        }
027        
028        public int getThreshold(int[] h) {
029                int K = h.length;
030                makeSigmaTables(h);     // set up S2_0, S2_1, N
031
032                int n0 = 0, n1;
033                int qMin = -1;
034                double eMin = Double.POSITIVE_INFINITY;
035                for (int q = 0; q <= K-2; q++) {
036                        n0 = n0 + h[q];
037                        n1 = N - n0;
038                        if (n0 > 0 && n1 > 0) { 
039                                double P0 = (double) n0 / N;    // could use n0, n1 instead
040                                double P1 = (double) n1 / N;
041                                double e = // 1.0 + 
042                                        P0 * Math.log(S2_0[q]) + P1 * Math.log(S2_1[q])
043                                        - 2 * (P0 * Math.log(P0) + P1 * Math.log(P1));
044                                if (e < eMin) {         // minimize e;
045                                        eMin = e;
046                                        qMin = q;
047                                }
048                        }
049                }
050                return qMin; 
051        }
052        
053        private void makeSigmaTables(int[] h) {
054                int K = h.length;
055                final double unitVar = 1d/12;   // variance of a uniform distribution in unit interval
056                S2_0 = new double[K];
057                S2_1 = new double[K];
058                
059                int n0 = 0;
060                long A0 = 0;
061                long B0 = 0; 
062                for (int q = 0; q < K; q++) {
063                        long ql = q;    // need a long type to avoid overflows
064                        n0 = n0 + h[q];
065                        A0 = A0 + h[q] * ql;
066                        B0 = B0 + h[q] * ql*ql;
067                        S2_0[q] = (n0 > 0) ? 
068                                        unitVar + ((double)B0 - (double)A0*A0/n0)/n0 : 0;       
069                }
070                
071                N = n0;
072                
073                int n1 = 0;
074                long A1 = 0;
075                long B1 = 0; 
076                S2_1[K-1] = 0;
077                for (int q = K-2; q >= 0; q--) {
078                        long qp1 = q+1;
079                        n1 = n1 + h[q+1];
080                        A1 = A1 + h[q+1] * qp1;
081                        B1 = B1 + h[q+1] * qp1*qp1;
082                        S2_1[q] = (n1 > 0) ? 
083                                        unitVar + ((double)B1 - (double)A1*A1/n1)/n1 : 0;
084                }
085        }
086        
087}