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 * Maximum entropy thresholder modeled after Kapur et al. (1985).
016 */
017public class MaxEntropyThresholder extends GlobalThresholder {
018        
019        static final double EPSILON = 1E-12;
020        
021        private double[] H0array = new double[256];     // only used for reporting
022        private double[] H1array = new double[256];     // only used for reporting
023        private double[] H01array = new double[256]; // only used for reporting
024        
025        private double[] S0 = null;
026        private double[] S1 = null;
027        
028        public MaxEntropyThresholder() {
029                super();
030        }
031        
032        public int getThreshold(int[] h) {
033                int K = h.length;       
034                double[] p = normalize(h);              // normalized histogram (probabilities)
035                makeTables(p);  // initialize S0, S1
036                
037                double P0 = 0, P1;
038                int qMax = -1;
039                double Hmax = Double.NEGATIVE_INFINITY;
040                
041                for (int q = 0; q <= K-1; q++) { // one more step for logging
042                        P0 = P0 + p[q]; 
043                        P1 = 1 - P0;    
044                        double H0 = (P0 > EPSILON) ? -S0[q]/P0 + Math.log(P0) : 0;                              
045                        double H1 = (P1 > EPSILON) ? -S1[q]/P1 + Math.log(P1) : 0;                      
046                        double H01 = H0 + H1;
047                        
048                        H0array[q] = H0;        // logging only
049                        H1array[q] = H1;        // logging only
050                        H01array[q] = H01;      // logging only
051                        
052                        if (H01 > Hmax) {
053                                Hmax = H01;
054                                qMax = q;
055                        }
056                }
057                return qMax;
058        }
059        
060        private void makeTables(double[] p) {
061                int K = p.length;
062
063                // make tables S0[], S1[]
064                S0 = new double[K];
065                S1 = new double[K];
066
067                double s0 = 0;
068                for (int i = 0; i < K; i++) {
069                        if (p[i] > EPSILON) {
070                                s0 = s0 + p[i] * Math.log(p[i]);
071                        }
072                        S0[i] = s0;
073                }
074
075                double s1 = 0;
076                for (int i = K-1; i >= 0; i--) {
077                        S1[i] = s1;
078                        if (p[i] > EPSILON) {
079                                s1 = s1 + p[i] * Math.log(p[i]);
080                        }
081                }
082        }
083}