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 * This thresholder implements the algorithm proposed by Ridler and Calvard (1978),
016 * T.W. Ridler, S. Calvard, Picture thresholding using an iterative selection method,
017 * IEEE Trans. System, Man and Cybernetics, SMC-8 (August 1978) 630-632.
018 * described in Glasbey/Horgan: "Image Analysis for the Biological Sciences" (Ch. 4).
019 * 
020 * Fast version using tables of background and foreground means.
021 */
022public class IsodataThresholder extends GlobalThresholder {
023        
024        private int MAX_ITERATIONS = 100;
025
026        private double[] M0 = null;             // table of background means
027        private double[] M1 = null;             // table of foreground means
028        
029        public IsodataThresholder() {
030                super();
031        }
032
033        public int getThreshold(int[] h) {
034                makeMeanTables(h);
035                int K = h.length;
036                int q = (int) M0[K-1];  // start with total mean
037                int q_;
038                
039                int i = 0;
040                do {
041                        i++; 
042                        if (M0[q]<0 || M1[q]<0)  // background or foreground is empty
043                                return -1;
044                        q_ = q;                         
045                        q = (int) ((M0[q] + M1[q])/2);
046                } while (q != q_ && i < MAX_ITERATIONS);
047                
048                return q;
049        }
050        
051        private int makeMeanTables(int[] h) {
052                int K = h.length;
053                M0 = new double[K];
054                M1 = new double[K];
055                int n0 = 0;
056                long s0 = 0;
057                for (int q = 0; q < K; q++) {
058                        n0 = n0 + h[q];
059                        s0 = s0 + q * h[q];
060                        M0[q] = (n0 > 0) ? ((double) s0)/n0 : -1;
061                }
062                int N = n0;
063                
064                int n1 = 0;
065                long s1 = 0;
066                M1[K-1] = 0;
067                for (int q = h.length-2; q >= 0; q--) {
068                        n1 = n1 + h[q+1];
069                        s1 = s1 + (q+1) * h[q+1];
070                        M1[q] = (n1 > 0) ? ((double) s1)/n1 : -1;
071                }
072                
073                return N;
074        }
075        
076}