package ij.process;
import ij.IJ;
import java.util.Arrays;

/** Autothresholding methods (limited to 256 bin histograms) from the Auto_Threshold plugin 
    (http://fiji.sc/Auto_Threshold) by G.Landini at bham dot ac dot uk). */
public class AutoThresholder {
    private static String[] mStrings;
            
    public enum Method {
        Default, 
        Huang,
        Intermodes,
        IsoData, 
        IJ_IsoData,
        Li, 
        MaxEntropy, 
        Mean, 
        MinError, 
        Minimum,
        Moments, 
        Otsu, 
        Percentile, 
        RenyiEntropy, 
        Shanbhag, 
        Triangle, 
        Yen
    };

    public static String[] getMethods() {
        if (mStrings==null) {
            Method[] mVals = Method.values();
            mStrings = new String[mVals.length];
            for (int i=0; i<mVals.length; i++)
                mStrings[i] = mVals[i].name();
        }
        return mStrings;
    }
    
    /** Calculates and returns a threshold using the specified
        method and 256 bin histogram. */
    public int getThreshold(Method method, int[] histogram) {
        if (histogram==null)
            throw new IllegalArgumentException("Histogram is null");
        if (histogram.length!=256)
            throw new IllegalArgumentException("Histogram length not 256");
        int threshold = 0;
        switch (method) {
            case Default: threshold =  defaultIsoData(histogram); break;
            case IJ_IsoData: threshold =  IJIsoData(histogram); break;
            case Huang: threshold = Huang(histogram); break;
            case Intermodes: threshold = Intermodes(histogram); break;
            case IsoData: threshold = IsoData(histogram); break;
            case Li: threshold = Li(histogram); break;
            case MaxEntropy: threshold = MaxEntropy(histogram); break;
            case Mean: threshold = Mean(histogram); break;
            case MinError: threshold = MinErrorI(histogram); break;
            case Minimum: threshold = Minimum(histogram); break;
            case Moments: threshold = Moments(histogram); break;
            case Otsu: threshold = Otsu(histogram); break;
            case Percentile: threshold = Percentile(histogram); break;
            case RenyiEntropy: threshold = RenyiEntropy(histogram); break;
            case Shanbhag: threshold = Shanbhag(histogram); break;
            case Triangle: threshold = Triangle(histogram); break;
            case Yen: threshold = Yen(histogram); break;
        }
        if (threshold==-1) threshold = 0;
        return threshold;
    }

    public int getThreshold(String mString, int[] histogram) {
        // throws an exception if unknown argument
        int index = mString.indexOf(" ");
        if (index!=-1)
            mString = mString.substring(0, index);
        Method method = Method.valueOf(Method.class, mString); 
        return getThreshold(method, histogram);
    }
    
    int Huang(int [] data ) {
        // Implements Huang's fuzzy thresholding method 
        // Uses Shannon's entropy function (one can also use Yager's entropy function) 
        // Huang L.-K. and Wang M.-J.J. (1995) "Image Thresholding by Minimizing  
        // the Measures of Fuzziness" Pattern Recognition, 28(1): 41-51
        // M. Emre Celebi  06.15.2007
        // Ported to ImageJ plugin by G. Landini from E Celebi's fourier_0.8 routines
        int threshold=-1;
        int ih, it;
        int first_bin;
        int last_bin;
        double sum_pix;
        double num_pix;
        double term;
        double ent;  // entropy 
        double min_ent; // min entropy 
        double mu_x;

        /* Determine the first non-zero bin */
        first_bin=0;
        for (ih = 0; ih < 256; ih++ ) {
            if ( data[ih] != 0 ) {
                first_bin = ih;
                break;
            }
        }

        /* Determine the last non-zero bin */
        last_bin=255;
        for (ih = 255; ih >= first_bin; ih-- ) {
            if ( data[ih] != 0 ) {
                last_bin = ih;
                break;
            }
        }
        term = 1.0 / ( double ) ( last_bin - first_bin );
        double [] mu_0 = new double[256];
        sum_pix = num_pix = 0;
        for ( ih = first_bin; ih < 256; ih++ ){
            sum_pix += (double)ih * data[ih];
            num_pix += data[ih];
            /* NUM_PIX cannot be zero ! */
            mu_0[ih] = sum_pix / num_pix;
        }

        double [] mu_1 = new double[256];
        sum_pix = num_pix = 0;
        for ( ih = last_bin; ih > 0; ih-- ){
            sum_pix += (double)ih * data[ih];
            num_pix += data[ih];
            /* NUM_PIX cannot be zero ! */
            mu_1[ih - 1] = sum_pix / ( double ) num_pix;
        }

        /* Determine the threshold that minimizes the fuzzy entropy */
        threshold = -1;
        min_ent = Double.MAX_VALUE;
        for ( it = 0; it < 256; it++ ){
            ent = 0.0;
            for ( ih = 0; ih <= it; ih++ ) {
                /* Equation (4) in Ref. 1 */
                mu_x = 1.0 / ( 1.0 + term * Math.abs ( ih - mu_0[it] ) );
                if ( !((mu_x  < 1e-06 ) || ( mu_x > 0.999999))) {
                    /* Equation (6) & (8) in Ref. 1 */
                    ent += data[ih] * ( -mu_x * Math.log ( mu_x ) - ( 1.0 - mu_x ) * Math.log ( 1.0 - mu_x ) );
                }
            }

            for ( ih = it + 1; ih < 256; ih++ ) {
                /* Equation (4) in Ref. 1 */
                mu_x = 1.0 / ( 1.0 + term * Math.abs ( ih - mu_1[it] ) );
                if ( !((mu_x  < 1e-06 ) || ( mu_x > 0.999999))) {
                    /* Equation (6) & (8) in Ref. 1 */
                    ent += data[ih] * ( -mu_x * Math.log ( mu_x ) - ( 1.0 - mu_x ) * Math.log ( 1.0 - mu_x ) );
                }
            }
            /* No need to divide by NUM_ROWS * NUM_COLS * LOG(2) ! */
            if ( ent < min_ent ) {
                min_ent = ent;
                threshold = it;
            }
        }
        return threshold;
    }

    boolean bimodalTest(double [] y) {
        int len=y.length;
        boolean b = false;
        int modes = 0;
 
        for (int k=1;k<len-1;k++){
            if (y[k-1] < y[k] && y[k+1] < y[k]) {
                modes++;
                if (modes>2)
                    return false;
            }
        }
        if (modes == 2)
            b = true;
        return b;
    }

    int Intermodes(int[] data ) {
        // J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
        // Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
        // ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
        // Original Matlab code Copyright (C) 2004 Antti Niemisto
        // See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
        // and the original Matlab code.
        //
        // Assumes a bimodal histogram. The histogram needs is smoothed (using a
        // running average of size 3, iteratively) until there are only two local maxima.
        // j and k
        // Threshold t is (j+k)/2.
        // Images with histograms having extremely unequal peaks or a broad and
        // flat valleys are unsuitable for this method.
        
        int minbin=-1, maxbin=-1;
        for (int i=0; i<data.length; i++)
            if (data[i]>0) maxbin = i;
        for (int i=data.length-1; i>=0; i--)
            if (data[i]>0) minbin = i;
        int length = (maxbin-minbin)+1;
        double [] hist = new double[length];
        for (int i=minbin; i<=maxbin; i++)
            hist[i-minbin] = data[i];
            
        int iter = 0;
        int threshold=-1;
        while (!bimodalTest(hist) ) {
             //smooth with a 3 point running mean filter
            double previous=0, current=0, next=hist[0];
            for (int i=0; i<length-1; i++) {
                previous = current;
                current = next;
                next = hist[i + 1];
                hist[i] = (previous+current+next)/3;
            }
            hist[length-1] = (current+next)/3;
            iter++;
            if (iter>10000) {
                threshold = -1;
                IJ.log("Intermodes Threshold not found after 10000 iterations.");
                return threshold;
            }
        }

        // The threshold is the mean between the two peaks.
        int tt=0;
        for (int i=1; i<length - 1; i++) {
            if (hist[i-1] < hist[i] && hist[i+1] < hist[i]){
                tt += i;
                //IJ.log("mode:" +i);
            }
        }
        threshold = (int) Math.floor(tt/2.0);
        return threshold+minbin;
    }

    int IsoData(int[] data ) {
        // Also called intermeans
        // Iterative procedure based on the isodata algorithm [T.W. Ridler, S. Calvard, Picture 
        // thresholding using an iterative selection method, IEEE Trans. System, Man and 
        // Cybernetics, SMC-8 (1978) 630-632.] 
        // The procedure divides the image into objects and background by taking an initial threshold,
        // then the averages of the pixels at or below the threshold and pixels above are computed. 
        // The averages of those two values are computed, the threshold is incremented and the 
        // process is repeated until the threshold is larger than the composite average. That is,
        //  threshold = (average background + average objects)/2
        // The code in ImageJ that implements this function is the getAutoThreshold() method in the ImageProcessor class. 
        //
        // From: Tim Morris (dtm@ap.co.umist.ac.uk)
        // Subject: Re: Thresholding method?
        // posted to sci.image.processing on 1996/06/24
        // The algorithm implemented in NIH Image sets the threshold as that grey
        // value, G, for which the average of the averages of the grey values
        // below and above G is equal to G. It does this by initialising G to the
        // lowest sensible value and iterating:

        // L = the average grey value of pixels with intensities < G
        // H = the average grey value of pixels with intensities > G
        // is G = (L + H)/2?
        // yes => exit
        // no => increment G and repeat
        //
        int i, l, totl, g=0;
        double toth, h;
        for (i = 1; i < 256; i++) {
            if (data[i] > 0){
                g = i + 1;
                break;
            }
        }
        while (true){
            l = 0;
            totl = 0;
            for (i = 0; i < g; i++) {
                 totl = totl + data[i];
                 l = l + (data[i] * i);
            }
            h = 0;
            toth = 0;
            for (i = g + 1; i < 256; i++){
                toth += data[i];
                h += ((double)data[i]*i);
            }
            if (totl > 0 && toth > 0){
                l /= totl;
                h /= toth;
                if (g == (int) Math.round((l + h) / 2.0))
                    break;
            }
            g++;
            if (g > 254)
                return -1;
        }
        return g;
    }
    
    int defaultIsoData(int[] data) {
        // This is the modified IsoData method used by the "Threshold" widget in "Default" mode
        int n = data.length;
        int[] data2 = new int[n];
        int mode=0, maxCount=0;
        for (int i=0; i<n; i++) {
            int count = data[i];
            data2[i] = data[i];
            if (data2[i]>maxCount) {
                maxCount = data2[i];
                mode = i;
            }
        }
        int maxCount2 = 0;
        for (int i = 0; i<n; i++) {
            if ((data2[i]>maxCount2) && (i!=mode))
                maxCount2 = data2[i];
        }
        int hmax = maxCount;
        if ((hmax>(maxCount2*2)) && (maxCount2!=0)) {
            hmax = (int)(maxCount2 * 1.5);
            data2[mode] = hmax;
        }
        return IJIsoData(data2);
    }

    int IJIsoData(int[] data) {
        // This is the original ImageJ IsoData implementation, here for backward compatibility.
        int level;
        int maxValue = data.length - 1;
        double result, sum1, sum2, sum3, sum4;
        int count0 = data[0];
        data[0] = 0; //set to zero so erased areas aren't included
        int countMax = data[maxValue];
        data[maxValue] = 0;
        int min = 0;
        while ((data[min]==0) && (min<maxValue))
            min++;
        int max = maxValue;
        while ((data[max]==0) && (max>0))
            max--;
        if (min>=max) {
            data[0]= count0; data[maxValue]=countMax;
            level = data.length/2;
            return level;
        }
        int movingIndex = min;
        int inc = Math.max(max/40, 1);
        do {
            sum1=sum2=sum3=sum4=0.0;
            for (int i=min; i<=movingIndex; i++) {
                sum1 += (double)i*data[i];
                sum2 += data[i];
            }
            for (int i=(movingIndex+1); i<=max; i++) {
                sum3 += (double)i*data[i];
                sum4 += data[i];
            }           
            result = (sum1/sum2 + sum3/sum4)/2.0;
            movingIndex++;
        } while ((movingIndex+1)<=result && movingIndex<max-1);
        data[0]= count0; data[maxValue]=countMax;
        level = (int)Math.round(result);
        return level;
    }


    int Li(int [] data ) {
        // Implements Li's Minimum Cross Entropy thresholding method
        // This implementation is based on the iterative version (Ref. 2) of the algorithm.
        // 1) Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" 
        //    Pattern Recognition, 26(4): 617-625
        // 2) Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum 
        //    Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776
        // 3) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 
        //    Techniques and Quantitative Performance Evaluation" Journal of 
        //    Electronic Imaging, 13(1): 146-165 
        //    http://citeseer.ist.psu.edu/sezgin04survey.html
        // Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
        int threshold;
        double num_pixels;
        double sum_back; /* sum of the background pixels at a given threshold */
        double sum_obj;  /* sum of the object pixels at a given threshold */
        double num_back; /* number of background pixels at a given threshold */
        double num_obj;  /* number of object pixels at a given threshold */
        double old_thresh;
        double new_thresh;
        double mean_back; /* mean of the background pixels at a given threshold */
        double mean_obj;  /* mean of the object pixels at a given threshold */
        double mean;  /* mean gray-level in the image */
        double tolerance; /* threshold tolerance */
        double temp;

        tolerance=0.5;
        num_pixels = 0;
        for (int ih = 0; ih < 256; ih++ ) 
            num_pixels += data[ih];

        /* Calculate the mean gray-level */
        mean = 0.0;
        for (int ih = 0 + 1; ih < 256; ih++ ) //0 + 1?
            mean += (double)ih * data[ih];
        mean /= num_pixels;
        /* Initial estimate */
        new_thresh = mean;

        do {
            old_thresh = new_thresh;
            threshold = (int) (old_thresh + 0.5);   /* range */
            /* Calculate the means of background and object pixels */
            /* Background */
            sum_back = 0;
            num_back = 0;
            for (int ih = 0; ih <= threshold; ih++ ) {
                sum_back += (double)ih * data[ih];
                num_back += data[ih];
            }
            mean_back = ( num_back == 0 ? 0.0 : ( sum_back / ( double ) num_back ) );
            /* Object */
            sum_obj = 0;
            num_obj = 0;
            for (int ih = threshold + 1; ih < 256; ih++ ) {
                sum_obj += (double)ih * data[ih];
                num_obj += data[ih];
            }
            mean_obj = ( num_obj == 0 ? 0.0 : ( sum_obj / ( double ) num_obj ) );

            /* Calculate the new threshold: Equation (7) in Ref. 2 */
            //new_thresh = simple_round ( ( mean_back - mean_obj ) / ( Math.log ( mean_back ) - Math.log ( mean_obj ) ) );
            //simple_round ( double x ) {
            // return ( int ) ( IS_NEG ( x ) ? x - .5 : x + .5 );
            //}
            //
            //#define IS_NEG( x ) ( ( x ) < -DBL_EPSILON ) 
            //DBL_EPSILON = 2.220446049250313E-16
            temp = ( mean_back - mean_obj ) / ( Math.log ( mean_back ) - Math.log ( mean_obj ) );

            if (temp < -2.220446049250313E-16)
                new_thresh = (int) (temp - 0.5);
            else
                new_thresh = (int) (temp + 0.5);
            /*  Stop the iterations when the difference between the
            new and old threshold values is less than the tolerance */
        }
        while ( Math.abs ( new_thresh - old_thresh ) > tolerance );
        return threshold;
    }

    int MaxEntropy(int [] data ) {
        // Implements Kapur-Sahoo-Wong (Maximum Entropy) thresholding method
        // Kapur J.N., Sahoo P.K., and Wong A.K.C. (1985) "A New Method for
        // Gray-Level Picture Thresholding Using the Entropy of the Histogram"
        // Graphical Models and Image Processing, 29(3): 273-285
        // M. Emre Celebi
        // 06.15.2007
        // Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
        int threshold=-1;
        int ih, it;
        int first_bin;
        int last_bin;
        double tot_ent;  /* total entropy */
        double max_ent;  /* max entropy */
        double ent_back; /* entropy of the background pixels at a given threshold */
        double ent_obj;  /* entropy of the object pixels at a given threshold */
        double [] norm_histo = new double[256]; /* normalized histogram */
        double [] P1 = new double[256]; /* cumulative normalized histogram */
        double [] P2 = new double[256]; 

        double total =0;
        for (ih = 0; ih < 256; ih++ ) 
            total+=data[ih];

        for (ih = 0; ih < 256; ih++ )
            norm_histo[ih] = data[ih]/total;

        P1[0]=norm_histo[0];
        P2[0]=1.0-P1[0];
        for (ih = 1; ih < 256; ih++ ){
            P1[ih]= P1[ih-1] + norm_histo[ih];
            P2[ih]= 1.0 - P1[ih];
        }

        /* Determine the first non-zero bin */
        first_bin=0;
        for (ih = 0; ih < 256; ih++ ) {
            if ( !(Math.abs(P1[ih])<2.220446049250313E-16)) {
                first_bin = ih;
                break;
            }
        }

        /* Determine the last non-zero bin */
        last_bin=255;
        for (ih = 255; ih >= first_bin; ih-- ) {
            if ( !(Math.abs(P2[ih])<2.220446049250313E-16)) {
                last_bin = ih;
                break;
            }
        }

        // Calculate the total entropy each gray-level
        // and find the threshold that maximizes it 
        max_ent = Double.MIN_VALUE;

        for ( it = first_bin; it <= last_bin; it++ ) {
            /* Entropy of the background pixels */
            ent_back = 0.0;
            for ( ih = 0; ih <= it; ih++ )  {
                if ( data[ih] !=0 ) {
                    ent_back -= ( norm_histo[ih] / P1[it] ) * Math.log ( norm_histo[ih] / P1[it] );
                }
            }

            /* Entropy of the object pixels */
            ent_obj = 0.0;
            for ( ih = it + 1; ih < 256; ih++ ){
                if (data[ih]!=0){
                ent_obj -= ( norm_histo[ih] / P2[it] ) * Math.log ( norm_histo[ih] / P2[it] );
                }
            }

            /* Total entropy */
            tot_ent = ent_back + ent_obj;

            // IJ.log(""+max_ent+"  "+tot_ent);
            if ( max_ent < tot_ent ) {
                max_ent = tot_ent;
                threshold = it;
            }
        }
        return threshold;
    }

    int Mean(int [] data ) {
        // C. A. Glasbey, "An analysis of histogram-based thresholding algorithms,"
        // CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
        //
        // The threshold is the mean of the greyscale data
        int threshold = -1;
        double tot=0, sum=0;
        for (int i=0; i<256; i++){
            tot+= data[i];
            sum+=((double)i*data[i]);
        }
        threshold =(int) Math.floor(sum/tot);
        return threshold;
    }

    int MinErrorI(int [] data ) {
        // Kittler and J. Illingworth, "Minimum error thresholding," Pattern Recognition, vol. 19, pp. 41-47, 1986.
        // C. A. Glasbey, "An analysis of histogram-based thresholding algorithms," CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
        // Ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
        // Original Matlab code Copyright (C) 2004 Antti Niemisto
        // See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
        // and the original Matlab code.

        int threshold = Mean(data); //Initial estimate for the threshold is found with the MEAN algorithm.
        int Tprev =-2;
        double mu, nu, p, q, sigma2, tau2, w0, w1, w2, sqterm, temp;
        //int counter=1;
        while (threshold!=Tprev){
            //Calculate some statistics.
            mu = B(data, threshold)/A(data, threshold);
            nu = (B(data, data.length - 1)-B(data, threshold))/(A(data, data.length - 1)-A(data, threshold));
            p = A(data, threshold)/A(data, data.length - 1);
            q = (A(data, data.length - 1)-A(data, threshold)) / A(data, data.length - 1);
            sigma2 = C(data, threshold)/A(data, threshold)-(mu*mu);
            tau2 = (C(data, data.length - 1)-C(data, threshold)) / (A(data, data.length - 1)-A(data, threshold)) - (nu*nu);

            //The terms of the quadratic equation to be solved.
            w0 = 1.0/sigma2-1.0/tau2;
            w1 = mu/sigma2-nu/tau2;
            w2 = (mu*mu)/sigma2 - (nu*nu)/tau2 + Math.log10((sigma2*(q*q))/(tau2*(p*p)));

            //If the next threshold would be imaginary, return with the current one.
            sqterm = (w1*w1)-w0*w2;
            if (sqterm < 0) {
                IJ.log("MinError(I): not converging.");
                return threshold;
            }

            //The updated threshold is the integer part of the solution of the quadratic equation.
            Tprev = threshold;
            temp = (w1+Math.sqrt(sqterm))/w0;

            if (Double.isNaN(temp)) {
                IJ.log ("MinError(I): NaN, not converging.");
                threshold = Tprev;
            } else
                threshold =(int) Math.floor(temp);
        }
        return threshold;
    }

    private double A(int[] y, int j) {
        if (j>=y.length) j=y.length-1;
        double x = 0;
        for (int i=0;i<=j;i++)
            x+=y[i];
        return x;
    }

    private double B(int[] y, int j) {
        if (j>=y.length) j=y.length-1;
        double x = 0;
        for (int i=0;i<=j;i++)
            x+=i*y[i];
        return x;
    }

    private double C(int[] y, int j) {
        if (j>=y.length) j=y.length-1;
        double x = 0;
        for (int i=0;i<=j;i++)
            x+=i*i*y[i];
        return x;
    }

    int Minimum(int [] data ) {
        // J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
        // Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
        // ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
        // Original Matlab code Copyright (C) 2004 Antti Niemisto
        // See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
        // and the original Matlab code.
        //
        // Assumes a bimodal histogram. The histogram needs is smoothed (using a
        // running average of size 3, iteratively) until there are only two local maxima.
        // Threshold t is such that yt-1 > yt <= yt+1.
        // Images with histograms having extremely unequal peaks or a broad and
        // flat valleys are unsuitable for this method.
        int iter =0;
        int threshold = -1;
        double [] iHisto = new double [256];
        for (int i=0; i<256; i++)
            iHisto[i]=(double) data[i];
        double [] tHisto = new double[iHisto.length] ;

        while (!bimodalTest(iHisto) ) {
             //smooth with a 3 point running mean filter
            for (int i=1; i<255; i++)
                tHisto[i]= (iHisto[i-1] + iHisto[i] +iHisto[i+1])/3;
            tHisto[0] = (iHisto[0]+iHisto[1])/3; //0 outside
            tHisto[255] = (iHisto[254]+iHisto[255])/3; //0 outside
            System.arraycopy(tHisto, 0, iHisto, 0, iHisto.length) ;
            iter++;
            if (iter>10000) {
                threshold = -1;
                IJ.log("Minimum: threshold not found after 10000 iterations.");
                return threshold;
            }
        }
        // The threshold is the minimum between the two peaks.
        for (int i=1; i<255; i++) {
            //IJ.log(" "+i+"  "+iHisto[i]);
            if (iHisto[i-1] > iHisto[i] && iHisto[i+1] >= iHisto[i]) {
                threshold = i;
                break;
            }
        }
        return threshold;
    }

    int Moments(int [] data ) {
        //  W. Tsai, "Moment-preserving thresholding: a new approach," Computer Vision,
        // Graphics, and Image Processing, vol. 29, pp. 377-393, 1985.
        // Ported to ImageJ plugin by G.Landini from the the open source project FOURIER 0.8
        // by  M. Emre Celebi , Department of Computer Science,  Louisiana State University in Shreveport
        // Shreveport, LA 71115, USA
        //  http://sourceforge.net/projects/fourier-ipal
        //  http://www.lsus.edu/faculty/~ecelebi/fourier.htm
        double total =0;
        double m0=1.0, m1=0.0, m2 =0.0, m3 =0.0, sum =0.0, p0=0.0;
        double cd, c0, c1, z0, z1;  /* auxiliary variables */
        int threshold = -1;

        double [] histo = new  double [256];

        for (int i=0; i<256; i++)
            total+=data[i];

        for (int i=0; i<256; i++)
            histo[i]=(double)(data[i]/total); //normalised histogram

        /* Calculate the first, second, and third order moments */
        for ( int i = 0; i < 256; i++ ) {
            double di = i;
            m1 += di * histo[i];
            m2 += di * di * histo[i];
            m3 += di * di * di * histo[i];
        }
        /* 
        First 4 moments of the gray-level image should match the first 4 moments
        of the target binary image. This leads to 4 equalities whose solutions 
        are given in the Appendix of Ref. 1 
        */
        cd = m0 * m2 - m1 * m1;
        c0 = ( -m2 * m2 + m1 * m3 ) / cd;
        c1 = ( m0 * -m3 + m2 * m1 ) / cd;
        z0 = 0.5 * ( -c1 - Math.sqrt ( c1 * c1 - 4.0 * c0 ) );
        z1 = 0.5 * ( -c1 + Math.sqrt ( c1 * c1 - 4.0 * c0 ) );
        p0 = ( z1 - m1 ) / ( z1 - z0 );  /* Fraction of the object pixels in the target binary image */

        // The threshold is the gray-level closest  
        // to the p0-tile of the normalized histogram 
        sum=0;
        for (int i=0; i<256; i++){
            sum+=histo[i];
            if (sum>p0) {
                threshold = i;
                break;
            }
        }
        return threshold;
    }

    int Otsu(int [] data ) {
        // Otsu's threshold algorithm
        // C++ code by Jordan Bevik <Jordan.Bevic@qtiworld.com>
        // ported to ImageJ plugin by G.Landini
        int k,kStar;  // k = the current threshold; kStar = optimal threshold
        double N1, N;    // N1 = # points with intensity <=k; N = total number of points
        double BCV, BCVmax; // The current Between Class Variance and maximum BCV
        double num, denom;  // temporary bookeeping
        double Sk;  // The total intensity for all histogram points <=k
        double S, L=256; // The total intensity of the image

        // Initialize values:
        S = N = 0;
        for (k=0; k<L; k++){
            S += (double)k * data[k];   // Total histogram intensity
            N += data[k];       // Total number of data points
        }

        Sk = 0;
        N1 = data[0]; // The entry for zero intensity
        BCV = 0;
        BCVmax=0;
        kStar = 0;

        // Look at each possible threshold value,
        // calculate the between-class variance, and decide if it's a max
        for (k=1; k<L-1; k++) { // No need to check endpoints k = 0 or k = L-1
            Sk += (double)k * data[k];
            N1 += data[k];

            // The float casting here is to avoid compiler warning about loss of precision and
            // will prevent overflow in the case of large saturated images
            denom = (double)( N1) * (N - N1); // Maximum value of denom is (N^2)/4 =  approx. 3E10

            if (denom != 0 ){
                // Float here is to avoid loss of precision when dividing
                num = ( (double)N1 / N ) * S - Sk;  // Maximum value of num =  255*N = approx 8E7
                BCV = (num * num) / denom;
            }
            else
                BCV = 0;

            if (BCV >= BCVmax){ // Assign the best threshold found so far
                BCVmax = BCV;
                kStar = k;
            }
        }
        // kStar += 1;  // Use QTI convention that intensity -> 1 if intensity >= k
        // (the algorithm was developed for I-> 1 if I <= k.)
        return kStar;
    }


    int Percentile(int [] data ) {
        // W. Doyle, "Operation useful for similarity-invariant pattern recognition,"
        // Journal of the Association for Computing Machinery, vol. 9,pp. 259-267, 1962.
        // ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
        // Original Matlab code Copyright (C) 2004 Antti Niemisto
        // See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
        // and the original Matlab code.

        int iter =0;
        int threshold = -1;
        double ptile= 0.5; // default fraction of foreground pixels
        double [] avec = new double [256];

        for (int i=0; i<256; i++)
            avec[i]=0.0;

        double total =partialSum(data, 255);
        double temp = 1.0;
        for (int i=0; i<256; i++){
            avec[i]=Math.abs((partialSum(data, i)/total)-ptile);
            //IJ.log("Ptile["+i+"]:"+ avec[i]);
            if (avec[i]<temp) {
                temp = avec[i];
                threshold = i;
            }
        }
        return threshold;
    }


    double partialSum(int [] y, int j) {
        double x = 0;
        for (int i=0;i<=j;i++)
            x+=y[i];
        return x;
    }


    int RenyiEntropy(int [] data ) {
        // Kapur J.N., Sahoo P.K., and Wong A.K.C. (1985) "A New Method for
        // Gray-Level Picture Thresholding Using the Entropy of the Histogram"
        // Graphical Models and Image Processing, 29(3): 273-285
        // M. Emre Celebi
        // 06.15.2007
        // Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines

        int threshold; 
        int opt_threshold;

        int ih, it;
        int first_bin;
        int last_bin;
        int tmp_var;
        int t_star1, t_star2, t_star3;
        int beta1, beta2, beta3;
        double alpha;/* alpha parameter of the method */
        double term;
        double tot_ent;  /* total entropy */
        double max_ent;  /* max entropy */
        double ent_back; /* entropy of the background pixels at a given threshold */
        double ent_obj;  /* entropy of the object pixels at a given threshold */
        double omega;
        double [] norm_histo = new double[256]; /* normalized histogram */
        double [] P1 = new double[256]; /* cumulative normalized histogram */
        double [] P2 = new double[256]; 

        double total =0;
        for (ih = 0; ih < 256; ih++ ) 
            total+=data[ih];

        for (ih = 0; ih < 256; ih++ )
            norm_histo[ih] = data[ih]/total;

        P1[0]=norm_histo[0];
        P2[0]=1.0-P1[0];
        for (ih = 1; ih < 256; ih++ ){
            P1[ih]= P1[ih-1] + norm_histo[ih];
            P2[ih]= 1.0 - P1[ih];
        }

        /* Determine the first non-zero bin */
        first_bin=0;
        for (ih = 0; ih < 256; ih++ ) {
            if ( !(Math.abs(P1[ih])<2.220446049250313E-16)) {
                first_bin = ih;
                break;
            }
        }

        /* Determine the last non-zero bin */
        last_bin=255;
        for (ih = 255; ih >= first_bin; ih-- ) {
            if ( !(Math.abs(P2[ih])<2.220446049250313E-16)) {
                last_bin = ih;
                break;
            }
        }

        /* Maximum Entropy Thresholding - BEGIN */
        /* ALPHA = 1.0 */
        /* Calculate the total entropy each gray-level
        and find the threshold that maximizes it 
        */
        threshold =0; // was MIN_INT in original code, but if an empty image is processed it gives an error later on.
        max_ent = 0.0;

        for ( it = first_bin; it <= last_bin; it++ ) {
            /* Entropy of the background pixels */
            ent_back = 0.0;
            for ( ih = 0; ih <= it; ih++ )  {
                if ( data[ih] !=0 ) {
                    ent_back -= ( norm_histo[ih] / P1[it] ) * Math.log ( norm_histo[ih] / P1[it] );
                }
            }

            /* Entropy of the object pixels */
            ent_obj = 0.0;
            for ( ih = it + 1; ih < 256; ih++ ){
                if (data[ih]!=0){
                ent_obj -= ( norm_histo[ih] / P2[it] ) * Math.log ( norm_histo[ih] / P2[it] );
                }
            }

            /* Total entropy */
            tot_ent = ent_back + ent_obj;

            // IJ.log(""+max_ent+"  "+tot_ent);

            if ( max_ent < tot_ent ) {
                max_ent = tot_ent;
                threshold = it;
            }
        }
        t_star2 = threshold;

        /* Maximum Entropy Thresholding - END */
        threshold =0; //was MIN_INT in original code, but if an empty image is processed it gives an error later on.
        max_ent = 0.0;
        alpha = 0.5;
        term = 1.0 / ( 1.0 - alpha );
        for ( it = first_bin; it <= last_bin; it++ ) {
            /* Entropy of the background pixels */
            ent_back = 0.0;
            for ( ih = 0; ih <= it; ih++ )
                ent_back += Math.sqrt ( norm_histo[ih] / P1[it] );

            /* Entropy of the object pixels */
            ent_obj = 0.0;
            for ( ih = it + 1; ih < 256; ih++ )
                ent_obj += Math.sqrt ( norm_histo[ih] / P2[it] );

            /* Total entropy */
            tot_ent = term * ( ( ent_back * ent_obj ) > 0.0 ? Math.log ( ent_back * ent_obj ) : 0.0);

            if ( tot_ent > max_ent ){
                max_ent = tot_ent;
                threshold = it;
            }
        }

        t_star1 = threshold;

        threshold = 0; //was MIN_INT in original code, but if an empty image is processed it gives an error later on.
        max_ent = 0.0;
        alpha = 2.0;
        term = 1.0 / ( 1.0 - alpha );
        for ( it = first_bin; it <= last_bin; it++ ) {
            /* Entropy of the background pixels */
            ent_back = 0.0;
            for ( ih = 0; ih <= it; ih++ )
                ent_back += ( norm_histo[ih] * norm_histo[ih] ) / ( P1[it] * P1[it] );

            /* Entropy of the object pixels */
            ent_obj = 0.0;
            for ( ih = it + 1; ih < 256; ih++ )
                ent_obj += ( norm_histo[ih] * norm_histo[ih] ) / ( P2[it] * P2[it] );

            /* Total entropy */
            tot_ent = term *( ( ent_back * ent_obj ) > 0.0 ? Math.log(ent_back * ent_obj ): 0.0 );

            if ( tot_ent > max_ent ){
                max_ent = tot_ent;
                threshold = it;
            }
        }

        t_star3 = threshold;

        /* Sort t_star values */
        if ( t_star2 < t_star1 ){
            tmp_var = t_star1;
            t_star1 = t_star2;
            t_star2 = tmp_var;
        }
        if ( t_star3 < t_star2 ){
            tmp_var = t_star2;
            t_star2 = t_star3;
            t_star3 = tmp_var;
        }
        if ( t_star2 < t_star1 ) {
            tmp_var = t_star1;
            t_star1 = t_star2;
            t_star2 = tmp_var;
        }

        /* Adjust beta values */
        if ( Math.abs ( t_star1 - t_star2 ) <= 5 )  {
            if ( Math.abs ( t_star2 - t_star3 ) <= 5 ) {
                beta1 = 1;
                beta2 = 2;
                beta3 = 1;
            }
            else {
                beta1 = 0;
                beta2 = 1;
                beta3 = 3;
            }
        }
        else {
            if ( Math.abs ( t_star2 - t_star3 ) <= 5 ) {
                beta1 = 3;
                beta2 = 1;
                beta3 = 0;
            }
            else {
                beta1 = 1;
                beta2 = 2;
                beta3 = 1;
            }
        }
        //IJ.log(""+t_star1+" "+t_star2+" "+t_star3);
        /* Determine the optimal threshold value */
        omega = P1[t_star3] - P1[t_star1];
        opt_threshold = (int) (t_star1 * ( P1[t_star1] + 0.25 * omega * beta1 ) + 0.25 * t_star2 * omega * beta2  + t_star3 * ( P2[t_star3] + 0.25 * omega * beta3 ));

        return opt_threshold;
    }


    int Shanbhag(int [] data ) {
        // Shanhbag A.G. (1994) "Utilization of Information Measure as a Means of
        //  Image Thresholding" Graphical Models and Image Processing, 56(5): 414-419
        // Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
        int threshold;
        int ih, it;
        int first_bin;
        int last_bin;
        double term;
        double tot_ent;  /* total entropy */
        double min_ent;  /* max entropy */
        double ent_back; /* entropy of the background pixels at a given threshold */
        double ent_obj;  /* entropy of the object pixels at a given threshold */
        double [] norm_histo = new double[256]; /* normalized histogram */
        double [] P1 = new double[256]; /* cumulative normalized histogram */
        double [] P2 = new double[256]; 

        double total =0;
        for (ih = 0; ih < 256; ih++ ) 
            total+=data[ih];

        for (ih = 0; ih < 256; ih++ )
            norm_histo[ih] = data[ih]/total;

        P1[0]=norm_histo[0];
        P2[0]=1.0-P1[0];
        for (ih = 1; ih < 256; ih++ ){
            P1[ih]= P1[ih-1] + norm_histo[ih];
            P2[ih]= 1.0 - P1[ih];
        }

        /* Determine the first non-zero bin */
        first_bin=0;
        for (ih = 0; ih < 256; ih++ ) {
            if ( !(Math.abs(P1[ih])<2.220446049250313E-16)) {
                first_bin = ih;
                break;
            }
        }

        /* Determine the last non-zero bin */
        last_bin=255;
        for (ih = 255; ih >= first_bin; ih-- ) {
            if ( !(Math.abs(P2[ih])<2.220446049250313E-16)) {
                last_bin = ih;
                break;
            }
        }

        // Calculate the total entropy each gray-level
        // and find the threshold that maximizes it 
        threshold =-1;
        min_ent = Double.MAX_VALUE;

        for ( it = first_bin; it <= last_bin; it++ ) {
            /* Entropy of the background pixels */
            ent_back = 0.0;
            term = 0.5 / P1[it];
            for ( ih = 1; ih <= it; ih++ )  { //0+1?
                ent_back -= norm_histo[ih] * Math.log ( 1.0 - term * P1[ih - 1] );
            }
            ent_back *= term;

            /* Entropy of the object pixels */
            ent_obj = 0.0;
            term = 0.5 / P2[it];
            for ( ih = it + 1; ih < 256; ih++ ){
                ent_obj -= norm_histo[ih] * Math.log ( 1.0 - term * P2[ih] );
            }
            ent_obj *= term;

            /* Total entropy */
            tot_ent = Math.abs ( ent_back - ent_obj );

            if ( tot_ent < min_ent ) {
                min_ent = tot_ent;
                threshold = it;
            }
        }
        return threshold;
    }


    int Triangle(int [] data ) {
        //  Zack, G. W., Rogers, W. E. and Latt, S. A., 1977,
        //  Automatic Measurement of Sister Chromatid Exchange Frequency,
        // Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753
        //
        //  modified from Johannes Schindelin plugin
        // 
        // find min and max
        int min = 0, dmax=0, max = 0, min2=0;
        for (int i = 0; i < data.length; i++) {
            if (data[i]>0){
                min=i;
                break;
            }
        }
        if (min>0) min--; // line to the (p==0) point, not to data[min]

        // The Triangle algorithm cannot tell whether the data is skewed to one side or another.
        // This causes a problem as there are 2 possible thresholds between the max and the 2 extremes
        // of the histogram.
        // Here I propose to find out to which side of the max point the data is furthest, and use that as
        //  the other extreme.
        for (int i = 255; i >0; i-- ) {
            if (data[i]>0){
                min2=i;
                break;
            }
        }
        if (min2<255) min2++; // line to the (p==0) point, not to data[min]

        for (int i =0; i < 256; i++) {
            if (data[i] >dmax) {
                max=i;
                dmax=data[i];
            }
        }
        // find which is the furthest side
        //IJ.log(""+min+" "+max+" "+min2);
        boolean inverted = false;
        if ((max-min)<(min2-max)){
            // reverse the histogram
            //IJ.log("Reversing histogram.");
            inverted = true;
            int left  = 0;          // index of leftmost element
            int right = 255; // index of rightmost element
            while (left < right) {
                // exchange the left and right elements
                int temp = data[left]; 
                data[left]  = data[right]; 
                data[right] = temp;
                // move the bounds toward the center
                left++;
                right--;
            }
            min=255-min2;
            max=255-max;
        }

        if (min == max){
            //IJ.log("Triangle:  min == max.");
            return min;
        }

        // describe line by nx * x + ny * y - d = 0
        double nx, ny, d;
        // nx is just the max frequency as the other point has freq=0
        nx = data[max];   //-min; // data[min]; //  lowest value bmin = (p=0)% in the image
        ny = min - max;
        d = Math.sqrt(nx * nx + ny * ny);
        nx /= d;
        ny /= d;
        d = nx * min + ny * data[min];

        // find split point
        int split = min;
        double splitDistance = 0;
        for (int i = min + 1; i <= max; i++) {
            double newDistance = nx * i + ny * data[i] - d;
            if (newDistance > splitDistance) {
                split = i;
                splitDistance = newDistance;
            }
        }
        split--;

        if (inverted) {
            // The histogram might be used for something else, so let's reverse it back
            int left  = 0; 
            int right = 255;
            while (left < right) {
                int temp = data[left]; 
                data[left]  = data[right]; 
                data[right] = temp;
                left++;
                right--;
            }
            return (255-split);
        }
        else
            return split;
    }


    int Yen(int [] data ) {
        // Implements Yen  thresholding method
        // 1) Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion 
        //    for Automatic Multilevel Thresholding" IEEE Trans. on Image 
        //    Processing, 4(3): 370-378
        // 2) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 
        //    Techniques and Quantitative Performance Evaluation" Journal of 
        //    Electronic Imaging, 13(1): 146-165
        //    http://citeseer.ist.psu.edu/sezgin04survey.html
        //
        // M. Emre Celebi
        // 06.15.2007
        // Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
        int threshold;
        int ih, it;
        double crit;
        double max_crit;
        double [] norm_histo = new double[256]; /* normalized histogram */
        double [] P1 = new double[256]; /* cumulative normalized histogram */
        double [] P1_sq = new double[256]; 
        double [] P2_sq = new double[256]; 

        double total =0;
        for (ih = 0; ih < 256; ih++ ) 
            total+=data[ih];

        for (ih = 0; ih < 256; ih++ )
            norm_histo[ih] = data[ih]/total;

        P1[0]=norm_histo[0];
        for (ih = 1; ih < 256; ih++ )
            P1[ih]= P1[ih-1] + norm_histo[ih];

        P1_sq[0]=norm_histo[0]*norm_histo[0];
        for (ih = 1; ih < 256; ih++ )
            P1_sq[ih]= P1_sq[ih-1] + norm_histo[ih] * norm_histo[ih];

        P2_sq[255] = 0.0;
        for ( ih = 254; ih >= 0; ih-- )
            P2_sq[ih] = P2_sq[ih + 1] + norm_histo[ih + 1] * norm_histo[ih + 1];

        /* Find the threshold that maximizes the criterion */
        threshold = -1;
        max_crit = Double.MIN_VALUE;
        for ( it = 0; it < 256; it++ ) {
            crit = -1.0 * (( P1_sq[it] * P2_sq[it] )> 0.0? Math.log( P1_sq[it] * P2_sq[it]):0.0) +  2 * ( ( P1[it] * ( 1.0 - P1[it] ) )>0.0? Math.log(  P1[it] * ( 1.0 - P1[it] ) ): 0.0);
            if ( crit > max_crit ) {
                max_crit = crit;
                threshold = it;
            }
        }
        return threshold;
    }

}