package yi.stat.hist;

/**
 * Class OtsuThresholding implements the Otsu Method automatic thresholding
 * algorithm.
 * <p>
 * Linda Shapiro and G. Stockman, Computer Vision, p85.
 * </p>
 */
public class OtsuThresholding {
  /**
   * The histogram
   */
  int[] m_histogram;

  /**
   * The summary of variances for the lastest thresholding
   */
  double m_sumOfVariances;

  /**
   * Mean of the probabilistic histogram
   */
  double m_mean;

  /**
   * Variance of the probabilistic histogram
   */
  double m_variance;

  int m_threshold;

  int m_unit;

  double m_probability1, m_probability2, m_mean1, m_mean2;

  /**
   * When the histogram is not begin from 0, use this one.
   * @param hist
   * @param offset
   * @return
   */
  public int thresholding(int[] hist, int offset) {
    thresholding(hist);
    m_mean1 += offset;
    m_mean2 += offset;
    m_mean += offset;
    m_threshold += offset;
    return m_threshold;
  }

  private int thresholding(int[] hist) {
    m_histogram = hist;
    int hist_size = hist.length;

    int unit = (int)yi.stat.Util.sum(hist);
    m_unit = unit;
    int prob1 = 0, prob2 = unit;
    double mean1 = 0, mean = yi.stat.Util.meanOfProbHist(hist), mean2 = mean;
    double diff_mean = (mean1 - mean2);
    double max_between_group_var = prob1 * prob2 * (diff_mean * diff_mean);
    int threshold = 0;
    for (int t=0; t<hist_size; t++) {  // t is the threshold, [0,t],(t,size)
      int new_prob1 = prob1 + hist[t];
      if (new_prob1 !=0 ) mean1 = (prob1 * mean1 + t * hist[t]) / new_prob1;
      prob1 = new_prob1;
      prob2 = unit - prob1;
      if (prob2 != 0) mean2 = (mean * unit - mean1 * new_prob1) / prob2;
      diff_mean = mean1 - mean2;
      double between_group_var = prob1 * prob2 * (diff_mean * diff_mean);
      if (max_between_group_var < between_group_var) {
        max_between_group_var = between_group_var;
        threshold = t;
        m_probability1 = (double)prob1 / unit;
        m_probability2 = (double)prob2 / unit;
        m_mean1 = mean1; m_mean2 = mean2;
      }
    }
    m_mean = mean;
    m_variance = yi.stat.Util.varianceOfProbHist(hist, mean);
    m_sumOfVariances = m_variance - max_between_group_var/ ((double)unit * unit);

    m_threshold = threshold + 1; // [0, threshold), [threshold, size)

    return m_threshold;
  }

  public String toString() {
    String return_str = System.getProperty("line.separator");
    double var1 = 0, var2 = 0;
    int total1 = 0, total2 = 0;
    for (int i=0; i<m_histogram.length; i++) {
      if (i < m_threshold) {
        double diff = i - m_mean1;
        var1 += m_histogram[i] * diff * diff;
        total1 += m_histogram[i];
      }
      else {
        double diff = i - m_mean2;
        var2 += m_histogram[i] * diff * diff;
        total2 += m_histogram[i];
      }
    }
    var1 = var1 / total1; var2 = var2 / total2;
    return new String(
      "threshold = " + m_threshold + return_str +
      "mean, variance, probability of group 1 = " + m_mean1 + " , " +
      var1 + " , " + m_probability1 + return_str +
      "mean, variance, probability of group 2 = " + m_mean2 + " , " +
      var2 + " , " + m_probability2 + return_str +
      "summary of variance = " + m_sumOfVariances + return_str +
      "mean, variance before split = " + m_mean + " , " + m_variance + return_str +
      "normalized gain = " + (m_sumOfVariances/m_variance) + return_str);
  }

  /**
   * Return the summary of variances
   * @return the summary of variances
   */
  public double getSumOfVariance() {
    return m_sumOfVariances;
  }

  /**
   * Return the mean
   * @return the mean
   */
  public double getMean() {
    return m_mean;
  }

  /**
   * Return the mean1
   * @return the mean1
   */
  public double getMean1() {
    return m_mean1;
  }

  /**
   * Return the mean2
   * @return the mean2
   */
  public double getMean2() {
    return m_mean2;
  }

  /**
   * Return the variance
   * @return the variance
   */
  public double getVariance() {
    return m_variance;
  }

  /**
   * Return the threshold
   * @return the threshold
   */
  public int getThreshold() {
    return m_threshold;
  }

  /**
   * Return the total number of data
   * @return the total number of data
   */
  public int getUnit() {
    return m_unit;
  }

  public int[] getHistogram() {
    return m_histogram;
  }

  public int getNumber1() {
    return (int)Math.round(m_unit * m_probability1);
  }

  public int getNumber2() {
    return (int)Math.round(m_unit * m_probability2);
  }

  public double getScore() {
    return m_sumOfVariances / m_variance;
//    double var1 = 0, var2 = 0;
//    int total1 = 0, total2 = 0;
//    for (int i=0; i<m_histogram.length; i++) {
//      if (i < m_threshold) {
//        double diff = i - m_mean1;
//        var1 += m_histogram[i] * diff * diff;
//        total1 += m_histogram[i];
//      }
//      else {
//        double diff = i - m_mean2;
//        var2 += m_histogram[i] * diff * diff;
//        total2 += m_histogram[i];
//      }
//    }
//    var1 = var1 / total1; var2 = var2 / total2;
//    double score = (Math.sqrt(var1) + Math.sqrt(var2))/ (m_mean2 - m_mean1);
//    return score;
  }

  /**
   * Test method thresholding
   */
  public static class testThresholding{
    public static void main(String[] args) {
      int[] data = {0, 2, 5, 10, 5, 3, 8, 12, 8, 2};
      OtsuThresholding otsu = new OtsuThresholding();
      int t = otsu.thresholding(data);
      System.out.println("threshold = " + t);
      System.out.println("sum of variances = " + otsu.getSumOfVariance());
      System.out.println("variance = " + otsu.getVariance());
      System.out.println("normalized sum of variances = " + otsu.getSumOfVariance()/otsu.getVariance());
      System.out.println(otsu.toString());

      double sum = yi.stat.Util.sum(data);
      double m = yi.stat.Util.meanOfProbHist(data);
      for (int i=1; i<data.length; i++) {
        System.out.print("t = " + i);
        double[] hist = new double[data.length];
        double q1 = 0, q2 = 0, m1 = 0, m2 = 0;
        for (int j=0; j<data.length; j++) {
          hist[j] = data[j] / (double) sum;
          if (j<i) {
            q1 += hist[j];
            m1 += j*hist[j];
          }
          else {
            q2 += hist[j];
            m2 += j* hist[j];
          }
        }
        m1 /= q1; m2 /= q2;
        System.out.print(" q1 = " + q1);
        System.out.print(" q2 = " + q2);
        System.out.print(" m1 = " + m1);
        System.out.print(" m2 = " + m2);
        double v1 = 0, v2 = 0;
        for (int j=0; j<data.length; j++) {
          if (j<i) {
            v1 += (j - m1) * (j - m1) * hist[j];
          }
          else {
            v2 += (j - m2) * (j - m2) * hist[j];
          }
        }
        v1 /= q1; v2 /= q2;
        System.out.print(" v1 = " + v1);
        System.out.println(" v2 = " + v2);
        double sum_of_var = v1 * q1 + v2 * q2;
        System.out.println("summary of variances = " + sum_of_var);
        double between_group_var = q1*(m1-m)*(m1-m) + q2*(m2-m)*(m2-m);
        System.out.println("between group variance = " + between_group_var);
        System.out.println("variances = " + (sum_of_var + between_group_var));

      }
    }
  }


}