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 imagingbook.common.histogram.HistogramUtils;
013
014/**
015 * <p>
016 * This is an implementation of the global thresholder proposed by Kapur et al. in [1]. See Sec. 9.1.5 (Alg. 9.5) of [2]
017 * for a detailed description.
018 * </p>
019 * <p>
020 * [1] J. N. Kapur, P. K. Sahoo, and A. K. C. Wong. A new method for gray-level picture thresholding using the entropy
021 * of the histogram. Computer Vision, Graphics, and Image Processing 29, 273–285 (1985). <br> [2] W. Burger, M.J. Burge,
022 * <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer (2022).
023 * </p>
024 *
025 * @author WB
026 * @version 2022/08/21
027 */
028public class MaxEntropyThresholder implements GlobalThresholder {
029        
030        static final double EPSILON = 1E-16;
031        
032        private double[] S0;
033        private double[] S1;
034        
035        public MaxEntropyThresholder() {
036                super();
037        }
038        
039        @Override
040        public float getThreshold(int[] h) {
041                int K = h.length;       
042                double[] p = HistogramUtils.pdf(h);                     // normalized histogram (to probabilities)      
043                makeTables(p);  // initialize S0, S1
044                
045                double P0 = 0;  // cumulative probability
046                
047                float qMax = NoThreshold;
048                double Hmax = Double.NEGATIVE_INFINITY;
049                
050                for (int q = 0; q <= K-2; q++) {                
051                        P0 = P0 + p[q]; 
052                        if (P0 < EPSILON) continue;             // empty histogram so far, nothing to do
053                                
054                        double P1 = 1 - P0;
055                        if (P1 < EPSILON) break;                // no more histogram entries, finished
056                        
057                        // P0, P1 are nonzero!                  
058                        double H0 = -S0[q]/P0 + Math.log(P0);
059                        double H1 = -S1[q]/P1 + Math.log(P1);
060                        double H01 = H0 + H1;           
061                        
062                        if (H01 > Hmax) {
063                                Hmax = H01;
064                                qMax = q;
065                        }               
066                }
067
068                return qMax;
069        }
070        
071        // pre-calculate tables S0 and S1
072        private void makeTables(double[] p) {
073                int K = p.length;               
074                S0 = new double[K];             
075                S1 = new double[K];             
076
077                // S0[q] = \sum_{i=0}^{q} p[i] * \log(p[i])
078                double s0 = 0;
079                for (int q = 0; q < K; q++) {
080                        if (p[q] > 0) {
081                                s0 = s0 + p[q] * Math.log(p[q]);
082                        }
083                        S0[q] = s0;
084                }
085                
086                // S1[q] = \sum_{i=q+1}^{K-1} p[i] * \log(p[i])
087                double s1 = 0;
088                for (int q = K-1; q >= 0; q--) {
089                        S1[q] = s1;
090                        if (p[q] > 0) {
091                                s1 = s1 + p[q] * Math.log(p[q]);
092                        }
093                }
094        }
095}