// otsu-image.c         Section: 3.7
// Threshold Selection by Maximizing Between Class Variance
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
//

#include "ImageIter.h"
#include "otsu-image.h"
#include "loop-ps.h"


///////////////////////////////////////////////////////////////////////
// The constructor for class: OtsuImage
// This class uses the loop_ps() constructor to initialize its
// point set.  Thus each point in its domain is equivalent to a
// possible partitioning of the normalized histogram.
//
OtsuImage::
OtsuImage(unsigned int tauCount,
	     IA_Image<IA_Point<int>,float> normalizedHistogram )
: IA_ClosureI<IA_Point<int>,float>(
//
///////////////////////////////////////////////////////////////////////
    loop_ps(tauCount,                              // Dimension
	    normalizedHistogram.domain().min()[0], // Start loop value
	    normalizedHistogram.domain().max()[0]) // End loop value
    ),
  d(tauCount),
  nh(normalizedHistogram),
  cnh(normalizedHistogram),    // We will correct the values shortly
  mu_cnh(normalizedHistogram), // ditto
  total_mu_squared(0),
  minNHpt(normalizedHistogram.domain().min()),
  maxNHpt(normalizedHistogram.domain().max())
{
    // Now correct value for cnh and mu_cnh
    IA_IPIter<IA_Point<int>,float>
	nhIter( nh );

    IA_Point<int>  ip; // Pixel iterator point temporary
    float  nhVal;    // Pixel iterator value temporary
    float cnhVal=0;  // For cumulative results
    while ( nhIter( ip, nhVal ) ) {
	cnhVal    += nhVal;
	cnh[ip]    = cnhVal;
	total_mu_squared += ip[0] * nhVal;  // Not complete until the end,
	mu_cnh[ip] = total_mu_squared;      // but it is the value here.
    }
    // Now that the loop is finished, total_mu_squared still must be ^2.
    total_mu_squared *= total_mu_squared;
};



///////////////////////////////////////////////////////////////////////
// OtsuImage::operator()
//    Define the image value as the objective function to be maximized.
//    NOTE:  tau_ip is a point who\'s co-ordinates represent a
//           possible partition of the image histogram.
//           tau_ip indices are one off from the book,
//           i.e. tau_ip[0] = tau   ...  tau_ip[i] = tau
//                               0                      i
float
OtsuImage::operator() (const IA_Point<int> &tau_ip) const {

    IA_Point<int>
	s = minNHpt[0], e = tau_ip[0];  // 1-D points

    double omega, mu, imValue = 0;

    // Calculate start & end values for class tau[j]
    for (int j=0; j<=this->d; j++) {
	s[0]  = (j==0      ) ? minNHpt[0] : tau_ip[j-1];
	e[0]  = (j==this->d) ? maxNHpt[0] : tau_ip[j];

	omega = (j==0) ? cnh(e)     : cnh(e) - cnh(s);

	if (omega != 0) {
	    mu = (j==0) ? mu_cnh(e) : mu_cnh(e) - mu_cnh(s);
	    imValue += mu * mu / omega;
	}
    }
    return imValue - total_mu_squared;
}


IA_ClosureI<IA_Point<int>,float> *
OtsuImage::clone_self() const {
    return new OtsuImage( *this );
}



///////////////////////////////////////////////////////////////////////
// OtsuImage::print_this()
//    For debugging purposes, we include a printing function
//    which prints a point value pair per output line.
ostream &
OtsuImage::print_this(ostream &os) const {

    IA_IPIter<IA_Point<int>,float>
	otsuIter( *this );

    IA_Point<int>  ip;   // Pixel iterator point temporary
    float        val;  // Pixel iterator value temporary
    while ( otsuIter( ip, val ) ) {
	os  << "(" << ip << "," << val << ")" << endl;
    }
    return os;
}
