#include "eyemBin.h"

static int Huang(int hist[256])
{
	// 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 (hist[ih] != 0) {
			first_bin = ih;
			break;
		}
	}

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

	double mu_1[256];
	sum_pix = num_pix = 0;
	for (ih = last_bin; ih > 0; ih--) {
		sum_pix += (double)ih * hist[ih];
		num_pix += hist[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 = DBL_MAX;
	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 * abs(ih - mu_0[it]));
			if (!((mu_x < 1e-06) || (mu_x > 0.999999))) {
				/* Equation (6) & (8) in Ref. 1 */
				ent += hist[ih] * (-mu_x * log(mu_x) - (1.0 - mu_x) * log(1.0 - mu_x));
			}
		}

		for (ih = it + 1; ih < 256; ih++) {
			/* Equation (4) in Ref. 1 */
			mu_x = 1.0 / (1.0 + term * abs(ih - mu_1[it]));
			if (!((mu_x < 1e-06) || (mu_x > 0.999999))) {
				/* Equation (6) & (8) in Ref. 1 */
				ent += hist[ih] * (-mu_x * log(mu_x) - (1.0 - mu_x) * 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;
}

static int IsoData(int hist[256])
{
	// 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 (hist[i] > 0) {
			g = i + 1;
			break;
		}
	}
	while (true) {
		l = 0;
		totl = 0;
		for (i = 0; i < g; i++) {
			totl = totl + hist[i];
			l = l + (hist[i] * i);
		}
		h = 0;
		toth = 0;
		for (i = g + 1; i < 256; i++) {
			toth += hist[i];
			h += ((double)hist[i] * i);
		}
		if (totl > 0 && toth > 0) {
			l /= totl;
			h /= toth;
			if (g == (int)round((l + h) / 2.0))
				break;
		}
		g++;
		if (g > 254)
			return -1;
	}
	return g;
}

static int Li(int hist[256])
{
	// 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 += hist[ih];

	/* Calculate the mean gray-level */
	mean = 0.0;
	for (int ih = 0 + 1; ih < 256; ih++) //0 + 1?
		mean += (double)ih * hist[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 * hist[ih];
			num_back += hist[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 * hist[ih];
			num_obj += hist[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) / (log(mean_back) - 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 (abs(new_thresh - old_thresh) > tolerance);
	return threshold;
}

static int MaxEntropy(int hist[256])
{
	// 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[256]; /* normalized histogram */
	double P1[256]; /* cumulative normalized histogram */
	double P2[256];

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

	for (ih = 0; ih < 256; ih++)
		norm_histo[ih] = hist[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 (!(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 (!(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 = DBL_MIN;

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

		/* Entropy of the object pixels */
		ent_obj = 0.0;
		for (ih = it + 1; ih < 256; ih++) {
			if (hist[ih] != 0) {
				ent_obj -= (norm_histo[ih] / P2[it]) * 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;
}

static int Mean(int hist[256])
{
	// 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 += hist[i];
		sum += ((double)i*hist[i]);
	}
	threshold = (int)floor(sum / tot);
	return threshold;
}

static int Moments(int hist[256])
{
	//  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[256];

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

	for (int i = 0; i < 256; i++)
		histo[i] = (double)(hist[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 - sqrt(c1 * c1 - 4.0 * c0));
	z1 = 0.5 * (-c1 + 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;
}

static int Otsu(int hist[])
{
	// 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 * hist[k];	// Total histogram intensity
		N += hist[k];		// Total number of data points
	}

	Sk = 0;
	N1 = hist[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 * hist[k];
		N1 += hist[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;
		}
	}
	return kStar;
}

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

static int Percentile(int hist[])
{
	// 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[256];

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

	double total = partialSum(hist, 255);
	double temp = 1.0;
	for (int i = 0; i < 256; i++) {
		avec[i] = abs((partialSum(hist, i) / total) - ptile);
		if (avec[i] < temp) {
			temp = avec[i];
			threshold = i;
		}
	}
	return threshold;
}

static int RenyiEntropy(int hist[])
{
	// 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[256]; /* normalized histogram */
	double P1[256]; /* cumulative normalized histogram */
	double P2[256];

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

	for (ih = 0; ih < 256; ih++)
		norm_histo[ih] = hist[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 (!(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 (!(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 (hist[ih] != 0) {
				ent_back -= (norm_histo[ih] / P1[it]) * log(norm_histo[ih] / P1[it]);
			}
		}

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

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

		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 += sqrt(norm_histo[ih] / P1[it]);

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

		/* Total entropy */
		tot_ent = term * ((ent_back * ent_obj) > 0.0 ? 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 ? 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 (abs(t_star1 - t_star2) <= 5) {
		if (abs(t_star2 - t_star3) <= 5) {
			beta1 = 1;
			beta2 = 2;
			beta3 = 1;
		}
		else {
			beta1 = 0;
			beta2 = 1;
			beta3 = 3;
		}
	}
	else {
		if (abs(t_star2 - t_star3) <= 5) {
			beta1 = 3;
			beta2 = 1;
			beta3 = 0;
		}
		else {
			beta1 = 1;
			beta2 = 2;
			beta3 = 1;
		}
	}
	/* 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;
}

static int Shanbhag(int hist[])
{
	// 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[256]; /* normalized histogram */
	double P1[256]; /* cumulative normalized histogram */
	double P2[256];

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

	for (ih = 0; ih < 256; ih++)
		norm_histo[ih] = hist[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 (!(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 (!(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 = DBL_MAX;

	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] * 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] * log(1.0 - term * P2[ih]);
		}
		ent_obj *= term;

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

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

static int Triangle(int hist[]) {
	//  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 < 256; i++) {
		if (hist[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 (hist[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 (hist[i] > dmax) {
			max = i;
			dmax = hist[i];
		}
	}
	// find which is the furthest side
	//IJ.log(""+min+" "+max+" "+min2);
	bool inverted = false;
	if ((max - min) < (min2 - max)) {
		// reverse the 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 = hist[left];
			hist[left] = hist[right];
			hist[right] = temp;
			// move the bounds toward the center
			left++;
			right--;
		}
		min = 255 - min2;
		max = 255 - max;
	}

	if (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 = hist[max];   //-min; // data[min]; //  lowest value bmin = (p=0)% in the image
	ny = min - max;
	d = sqrt(nx * nx + ny * ny);
	nx /= d;
	ny /= d;
	d = nx * min + ny * hist[min];

	// find split point
	int split = min;
	double splitDistance = 0;
	for (int i = min + 1; i <= max; i++) {
		double newDistance = nx * i + ny * hist[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 = hist[left];
			hist[left] = hist[right];
			hist[right] = temp;
			left++;
			right--;
		}
		return (255 - split);
	}
	else
		return split;
}

static int Yen(int hist[])
{
	// 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[256]; /* normalized histogram */
	double P1[256]; /* cumulative normalized histogram */
	double P1_sq[256];
	double P2_sq[256];

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

	for (ih = 0; ih < 256; ih++)
		norm_histo[ih] = hist[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 = DBL_MIN;
	for (it = 0; it < 256; it++) {
		crit = -1.0 * ((P1_sq[it] * P2_sq[it]) > 0.0 ? log(P1_sq[it] * P2_sq[it]) : 0.0) + 2 * ((P1[it] * (1.0 - P1[it])) > 0.0 ? log(P1[it] * (1.0 - P1[it])) : 0.0);
		if (crit > max_crit) {
			max_crit = crit;
			threshold = it;
		}
	}
	return threshold;
}

static double getThreshVal_Otsu_8u(const cv::Mat& _src)
{
	cv::Size size = _src.size();
	int step = (int)_src.step;
	if (_src.isContinuous())
	{
		size.width *= size.height;
		size.height = 1;
		step = size.width;
	}

#ifdef HAVE_IPP
	unsigned char thresh = 0;
	CV_IPP_RUN_FAST(ipp_getThreshVal_Otsu_8u(_src.ptr(), step, size, thresh), thresh);
#endif

	const int N = 256;
	int i, j, h[N] = { 0 };
#if CV_ENABLE_UNROLLED
	int h_unrolled[3][N] = {};
#endif
	for (i = 0; i < size.height; i++)
	{
		const uchar* src = _src.ptr() + step*i;
		j = 0;
#if CV_ENABLE_UNROLLED
		for (; j <= size.width - 4; j += 4)
		{
			int v0 = src[j], v1 = src[j + 1];
			h[v0]++; h_unrolled[0][v1]++;
			v0 = src[j + 2]; v1 = src[j + 3];
			h_unrolled[1][v0]++; h_unrolled[2][v1]++;
		}
#endif
		for (; j < size.width; j++)
			h[src[j]]++;
	}

	double mu = 0, scale = 1. / (size.width*size.height);
	for (i = 0; i < N; i++)
	{
#if CV_ENABLE_UNROLLED
		h[i] += h_unrolled[0][i] + h_unrolled[1][i] + h_unrolled[2][i];
#endif
		mu += i*(double)h[i];
	}

	mu *= scale;
	double mu1 = 0, q1 = 0;
	double max_sigma = 0, max_val = 0;

	for (i = 0; i < N; i++)
	{
		double p_i, q2, mu2, sigma;

		p_i = h[i] * scale;
		mu1 *= q1;
		q1 += p_i;
		q2 = 1. - q1;

		if (std::min(q1, q2) < FLT_EPSILON || std::max(q1, q2) > 1. - FLT_EPSILON)
			continue;

		mu1 = (mu1 + i*p_i) / q1;
		mu2 = (mu - q1*mu1) / q2;
		sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
		if (sigma > max_sigma)
		{
			max_sigma = sigma;
			max_val = i;
		}
	}

	return max_val;
}

int eyemBinThreshold(EyemImage tpSrcImg, int iLightDark, double dThresh, double dMaxVal, EyemImage *tpDstImg)
{
	cv::Mat image = cv::Mat(tpSrcImg.iHeight, tpSrcImg.iWidth, MAKETYPE(tpSrcImg.iDepth, tpSrcImg.iChannels), tpSrcImg.vpImage).clone();
	if (image.empty()) {
		return FUNC_IMAGE_NOT_EXIST;
	}
	cv::Mat binary;
	//执行二值化操作
	cv::threshold(image, binary, dThresh, dMaxVal, iLightDark);
	tpDstImg->iWidth = binary.cols; tpDstImg->iHeight = binary.rows; tpDstImg->iDepth = binary.depth(); tpDstImg->iChannels = binary.channels();

	//内存尺寸
	int _Size = tpDstImg->iWidth*tpDstImg->iHeight*tpDstImg->iChannels * sizeof(uint8_t);
	//分配初始化内存
	tpDstImg->vpImage = (uint8_t *)malloc(_Size);
	if (NULL == tpDstImg->vpImage)
		return FUNC_NOT_ENOUGH_MEM;
	memset(tpDstImg->vpImage, 0, _Size);
	//拷贝数据
	memcpy(tpDstImg->vpImage, binary.data, _Size);

	return FUNC_OK;
}

int	eyemBinNiBlack(EyemImage tpSrcImg, int iType, int iWinSize, double dK, int binMethod, double dR, EyemImage *tpDstImg)
{
	cv::Mat src = cv::Mat(tpSrcImg.iHeight, tpSrcImg.iWidth, MAKETYPE(tpSrcImg.iDepth, tpSrcImg.iChannels), tpSrcImg.vpImage).clone();

	if (src.empty()) {
		return FUNC_IMAGE_NOT_EXIST;
	}

	CV_Assert(iWinSize % 2 == 1 && iWinSize > 1);
	if (binMethod == BINARIZATION_SAUVOLA) {
		CV_Assert(src.depth() == CV_8U);
		CV_Assert(dR != 0);
	}
	iType &= cv::THRESH_MASK;

	cv::Mat thresh;
	{
		cv::Mat mean, sqmean, variance, stddev, sqrtVarianceMeanSum;
		double srcMin, stddevMax;
		boxFilter(src, mean, CV_32F, cv::Size(iWinSize, iWinSize),
			cv::Point(-1, -1), true, cv::BORDER_REPLICATE);
		sqrBoxFilter(src, sqmean, CV_32F, cv::Size(iWinSize, iWinSize),
			cv::Point(-1, -1), true, cv::BORDER_REPLICATE);
		variance = sqmean - mean.mul(mean);
		sqrt(variance, stddev);
		switch (binMethod)
		{
		case BINARIZATION_NIBLACK:
			thresh = mean + stddev * static_cast<float>(dK);
			break;
		case BINARIZATION_SAUVOLA:
			thresh = mean.mul(1. + static_cast<float>(dK) * (stddev / dR - 1.));
			break;
		case BINARIZATION_WOLF:
			minMaxIdx(src, &srcMin);
			minMaxIdx(stddev, NULL, &stddevMax);
			thresh = mean - static_cast<float>(dK) * (mean - srcMin - stddev.mul(mean - srcMin) / stddevMax);
			break;
		case BINARIZATION_NICK:
			sqrt(variance + sqmean, sqrtVarianceMeanSum);
			thresh = mean + static_cast<float>(dK) * sqrtVarianceMeanSum;
			break;
		default:
			break;
		}
		thresh.convertTo(thresh, src.depth());
	}

	cv::Mat dst(src.size(), CV_8U);

	cv::Mat mask;
	switch (iType)
	{
	case cv::THRESH_BINARY:
	case cv::THRESH_BINARY_INV:
		compare(src, thresh, mask, (iType == cv::THRESH_BINARY ? cv::CMP_GT : cv::CMP_LE));
		dst.setTo(0);
		dst.setTo(255, mask);
		break;
	case cv::THRESH_TRUNC:
		compare(src, thresh, mask, cv::CMP_GT);
		src.copyTo(dst);
		thresh.copyTo(dst, mask);
		break;
	case cv::THRESH_TOZERO:
	case cv::THRESH_TOZERO_INV:
		compare(src, thresh, mask, (iType == cv::THRESH_TOZERO ? cv::CMP_GT : cv::CMP_LE));
		dst.setTo(0);
		src.copyTo(dst, mask);
		break;
	default:
		break;
	}

	//输出图像
	{
		tpDstImg->iWidth = dst.cols; tpDstImg->iHeight = dst.rows; tpDstImg->iDepth = dst.depth(); tpDstImg->iChannels = dst.channels();

		//内存尺寸
		int _Size = tpDstImg->iWidth*tpDstImg->iHeight*tpDstImg->iChannels * sizeof(uint8_t);

		//分配内存
		tpDstImg->vpImage = (uint8_t *)malloc(_Size);
		if (NULL == tpDstImg->vpImage)
			return FUNC_NOT_ENOUGH_MEM;
		memset(tpDstImg->vpImage, 0, _Size);

		//拷贝数据
		memcpy(tpDstImg->vpImage, dst.data, _Size);
	}

	return FUNC_OK;
}

int eyemBinDynThreshold(EyemImage tpSrcImg, EyemImage tpPreImg, double dOffset, int iType, EyemImage *tpDstImg)
{
	CV_Assert(MAKETYPE(tpSrcImg.iDepth, tpSrcImg.iChannels) == MAKETYPE(tpPreImg.iDepth, tpPreImg.iChannels));

	cv::Mat image = cv::Mat(tpSrcImg.iHeight, tpSrcImg.iWidth, MAKETYPE(tpSrcImg.iDepth, tpSrcImg.iChannels), tpSrcImg.vpImage).clone();

	if (image.empty()) {
		return FUNC_IMAGE_NOT_EXIST;
	}

	cv::Mat imagePre, variance;
	cv::Mat thresh;
	{
		imagePre = cv::Mat(tpPreImg.iHeight, tpPreImg.iWidth, MAKETYPE(tpPreImg.iDepth, tpPreImg.iChannels), tpPreImg.vpImage).clone();
		switch (iType)
		{
		case LIGHT:
			variance = image - imagePre;
			break;
		case DARK:
			variance = imagePre - image;
			break;
		case EQUAL:
			variance = abs(image - imagePre);
			break;
		case NOT_EQUAL:
			variance = abs(imagePre - image);
			break;
		default:
			break;
		}
	}
	cv::Mat binary, showMat;
	cv::compare(variance, cv::Mat::ones(imagePre.size(), imagePre.type())*dOffset, binary, cv::CMP_GT);

	//输出结果图像
	{
		if (NULL != tpDstImg->vpImage) {
			tpDstImg->iWidth = tpDstImg->iHeight = tpDstImg->iDepth = tpDstImg->iChannels = 0;
			//释放
			free(tpDstImg->vpImage);
			tpDstImg->vpImage = NULL;
		}

		tpDstImg->iWidth = binary.cols; tpDstImg->iHeight = binary.rows; tpDstImg->iDepth = binary.depth(); tpDstImg->iChannels = binary.channels();

		//内存尺寸
		int _Size = tpDstImg->iWidth*tpDstImg->iHeight*tpDstImg->iChannels * sizeof(uint8_t);

		//分配初始化内存
		tpDstImg->vpImage = (uint8_t *)malloc(_Size);
		if (NULL == tpDstImg->vpImage)
			return FUNC_NOT_ENOUGH_MEM;
		memset(tpDstImg->vpImage, 0, _Size);

		//拷贝数据
		memcpy(tpDstImg->vpImage, binary.data, _Size);
	}
	return FUNC_OK;
}

int eyemBinAutoThreshold(EyemImage tpImage, double dSigma, int iLightDark, int binMethod, EyemImage *tpDstImg)
{
	cv::Mat image = cv::Mat(tpImage.iHeight, tpImage.iWidth, MAKETYPE(tpImage.iDepth, tpImage.iChannels), tpImage.vpImage).clone();

	if (image.empty()) {
		return FUNC_IMAGE_NOT_EXIST;
	}
	int(*calc_threshold_param) (int *) = 0;

	int threshold = 0;
	switch (binMethod)
	{
	case HUANG:
		calc_threshold_param = Huang;
		break;
	case ISODATA:
		calc_threshold_param = IsoData;
		break;
	case LI:
		calc_threshold_param = Li;
		break;
	case MAXENTROPY:
		calc_threshold_param = MaxEntropy;
		break;
	case MEAN:
		calc_threshold_param = Mean;
		break;
	case MOMENTS:
		calc_threshold_param = Moments;
		break;
	case OTSU:
		calc_threshold_param = Otsu;
		break;
	default:
		calc_threshold_param = Otsu;
		break;
	}
	int hist[256];
	std::memset(hist, 0, sizeof(hist));

	//计算直方图
	for (int Y = 0; Y < 256; Y++) hist[Y] = 0;
	for (int Y = 0; Y < image.rows; Y++)
	{
		uchar *uPtr = image.data + Y * image.cols;
		for (int X = 0; X < image.cols; X++, uPtr++) hist[*uPtr]++;
	}

	threshold = calc_threshold_param(hist);

	cv::Mat binary;
	cv::threshold(image, binary, threshold, 255, iLightDark);

	//输出结果图像
	{
		if (NULL != tpDstImg->vpImage) {
			tpDstImg->iWidth = tpDstImg->iHeight = tpDstImg->iDepth = tpDstImg->iChannels = 0;
			//释放
			free(tpDstImg->vpImage);
			tpDstImg->vpImage = NULL;
		}

		tpDstImg->iWidth = binary.cols; tpDstImg->iHeight = binary.rows; tpDstImg->iDepth = binary.depth(); tpDstImg->iChannels = binary.channels();

		//内存尺寸
		int _Size = tpDstImg->iWidth*tpDstImg->iHeight*tpDstImg->iChannels * sizeof(uint8_t);

		//分配初始化内存
		tpDstImg->vpImage = (uint8_t *)malloc(_Size);
		if (NULL == tpDstImg->vpImage)
			return FUNC_NOT_ENOUGH_MEM;
		memset(tpDstImg->vpImage, 0, _Size);

		//拷贝数据
		memcpy(tpDstImg->vpImage, binary.data, _Size);
	}

	return FUNC_OK;
}

int eyemBinThresholdC(EyemImage tpImage, EyemHSVModel tpHSVModel, EyemImage *tpDstImg)
{
	cv::Mat image = cv::Mat(tpImage.iHeight, tpImage.iWidth, MAKETYPE(tpImage.iDepth, tpImage.iChannels), tpImage.vpImage).clone();

	if (image.empty())
		return FUNC_IMAGE_NOT_EXIST;
	//图像尺寸
	const int X = image.cols, Y = image.rows;
	//非彩色图像处理
	int incn = image.channels();
	if (incn > 3) {
		cv::cvtColor(image, image, cv::COLOR_BGRA2BGR);
	}
	else if (incn == 1) {
		cv::cvtColor(image, image, cv::COLOR_GRAY2BGR);
	}
	//转hsv空间
	cv::Mat imghsv;
	cv::cvtColor(image, imghsv, cv::COLOR_BGR2HSV);

	//红色比较特殊，分两个区间
	cv::Mat mask1, mask2(cv::Size(X, Y), CV_8UC1, cv::Scalar(0));
	cv::inRange(imghsv, cv::Scalar(tpHSVModel.dpRangeL[0], tpHSVModel.dpRangeL[1], tpHSVModel.dpRangeL[2]),
		cv::Scalar(tpHSVModel.dpRangeU[0], tpHSVModel.dpRangeU[1], tpHSVModel.dpRangeU[2]), mask1);
	//多个分割阈值
	if ((tpHSVModel.dpRangeLExt[0] + tpHSVModel.dpRangeLExt[1] + tpHSVModel.dpRangeLExt[2]) != 0 ||
		(tpHSVModel.dpRangeUExt[0] + tpHSVModel.dpRangeUExt[1] + tpHSVModel.dpRangeUExt[2]) != 0) {
		cv::inRange(imghsv, cv::Scalar(tpHSVModel.dpRangeLExt[0], tpHSVModel.dpRangeLExt[1], tpHSVModel.dpRangeLExt[2]),
			cv::Scalar(tpHSVModel.dpRangeUExt[0], tpHSVModel.dpRangeUExt[1], tpHSVModel.dpRangeUExt[2]), mask2);
	}
	//合并
	cv::Mat maskj;
	cv::bitwise_or(mask1, mask2, maskj);

	//输出结果图像
	if (NULL != tpDstImg->vpImage) {
		tpDstImg->iWidth = tpDstImg->iHeight = tpDstImg->iDepth = tpDstImg->iChannels = 0;
		//释放
		free(tpDstImg->vpImage);
		tpDstImg->vpImage = NULL;
	}

	tpDstImg->iWidth = maskj.cols; tpDstImg->iHeight = maskj.rows; tpDstImg->iDepth = maskj.depth(); tpDstImg->iChannels = maskj.channels();

	//内存尺寸
	int _Size = tpDstImg->iWidth*tpDstImg->iHeight*tpDstImg->iChannels * sizeof(uint8_t);

	//分配初始化内存
	tpDstImg->vpImage = (uint8_t *)malloc(_Size);
	if (NULL == tpDstImg->vpImage)
		return FUNC_NOT_ENOUGH_MEM;
	memset(tpDstImg->vpImage, 0, _Size);

	//拷贝数据
	memcpy(tpDstImg->vpImage, maskj.data, _Size);

	return FUNC_OK;
}

int	eyemBinDilation(EyemImage tpSrcImg, int iBinLevel, int iNum, EyemImage *tpDstImg)
{

	return FUNC_OK;
}

int	eyemBinErosion(EyemImage tpSrcImg, int iBinLevel, int iNum, EyemImage *tpDstImg)
{
	return FUNC_OK;
}

int	eyemBinOpening(EyemImage tpSrcImg, int iBinLevel, int iNum, EyemImage *tpDstImg)
{
	return FUNC_OK;
}

int	eyemBinClosing(EyemImage tpSrcImg, int iBinLevel, int iNum, EyemImage *tpDstImg)
{
	return FUNC_OK;
}

int eyemBinBlob(EyemImage tpImage, IntPtr *hObject, int iAreaThrs, EyemBinBlob **tpResult, int *ipNum, EyemImage *tpDstImg)
{
	cv::Mat image = cv::Mat(tpImage.iHeight, tpImage.iWidth, MAKETYPE(tpImage.iDepth, tpImage.iChannels), tpImage.vpImage).clone();
	if (image.empty()) {
		return FUNC_IMAGE_NOT_EXIST;
	}
	//判断图像
	if (image.type() != CV_8UC1 || image.channels() != 1) {
		return FUNC_CANNOT_CALC;
	}

	cv::threshold(image, image, 0, 255, cv::THRESH_BINARY_INV | cv::THRESH_OTSU);

	//显示图像
	cv::Mat showMat;
	cv::cvtColor(image, showMat, cv::COLOR_GRAY2RGB);
	//图像尺寸
	const int X = image.cols, Y = image.rows;

	//斑点大小限制
	bool filterByArea = true;
	int minArea = 25, maxArea = 25000;

	//斑点圆度限制
	bool filterByCircularity = false;
	float minCircularity = 0.8f, maxCircularity = std::numeric_limits<float>::max();

	//斑点的惯性率限制
	bool filterByInertia = false;
	float minInertiaRatio = 0.1f, maxInertiaRatio = std::numeric_limits<float>::max();

	//斑点凸度限制
	bool filterByConvexity = false;
	float minConvexity = 0.8f, maxConvexity = std::numeric_limits<float>::max();

	//斑点检测
	cv::Mat labels, stats, centroids;
	int nccomps = cv::connectedComponentsWithStats(image, labels, stats, centroids);

	std::vector<uchar> colors(nccomps + 1, 0);
	//按面积过滤
	if (filterByArea) {
		//过滤连通域面积
		for (int i = 0; i < nccomps; i++) {
			colors[i] = 255;
			double dArea = stats.ptr<int>(i)[cv::CC_STAT_AREA];
			if (!(dArea >= minArea&&dArea <= maxArea)) {
				colors[i] = 0;
			}
		}
	}

	//斑点方位
	cv::Mat mOrientation(cv::Size(1, nccomps + 1), CV_64FC1, cv::Scalar(0));
	//根据轮廓属性过滤
	std::vector<std::vector<cv::Point>> contours;
	cv::findContours(image, contours, cv::RETR_LIST, cv::CHAIN_APPROX_NONE);
	for (auto&contour : contours)
	{
		int label = labels.at<int>(contour[0]);
		//计算轮廓矩
		cv::Moments moms = cv::moments(contour);
		//主要方向计算
		cv::Mat pts((int)contour.size(), 2, CV_64FC1);
		for (int i = 0; i < pts.rows; i++)
		{
			pts.ptr<double>(i)[0] = contour[i].x;
			pts.ptr<double>(i)[1] = contour[i].y;
		}

		if (pts.rows > 2) {
			cv::PCA pca(pts, cv::Mat(), cv::PCA::DATA_AS_ROW);

			cv::Point pt((int)pca.mean.at<double>(0, 0), (int)pca.mean.at<double>(0, 1));

			//特征值和特征向量
			std::vector<cv::Point2d> eigvec(2);
			std::vector<double> eigval(2);
			for (int i = 0; i < 2; ++i)
			{
				eigvec[i] = cv::Point2d(pca.eigenvectors.at<double>(i, 0), pca.eigenvectors.at<double>(i, 1));
				eigval[i] = pca.eigenvalues.at<double>(i, 0);
			}

			mOrientation.ptr<float>(label)[0] = (float)atan2(eigvec[0].y, eigvec[0].x);
		}

		//按圆度过滤
		if (filterByCircularity) {
			double perimeter = cv::arcLength(contour, true);
			double ratio = 4 * CV_PI * moms.m00 / (perimeter * perimeter);
			if (ratio < minCircularity || ratio >= maxCircularity)
				colors[label] = 0;
		}
		//按惯性率过滤
		if (filterByInertia) {
			double denominator = std::sqrt(std::pow(2 * moms.mu11, 2) + std::pow(moms.mu20 - moms.mu02, 2));
			const double eps = 1e-2;
			double ratio;
			if (denominator > eps) {
				double cosmin = (moms.mu20 - moms.mu02) / denominator;
				double sinmin = 2 * moms.mu11 / denominator;
				double cosmax = -cosmin;
				double sinmax = -sinmin;

				double imin = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmin - moms.mu11 * sinmin;
				double imax = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmax - moms.mu11 * sinmax;
				ratio = imin / imax;
			}
			else {
				ratio = 1;
			}
			if (ratio < minInertiaRatio || ratio >= maxInertiaRatio)
				colors[label] = 0;
		}
		//按凸度过滤
		if (filterByConvexity) {
			std::vector <cv::Point> hull;
			cv::convexHull(contour, hull);
			double area = cv::contourArea(contour);
			double hullArea = contourArea(hull);
			if (fabs(hullArea) < DBL_EPSILON)
				colors[label] = 0;
			double ratio = area / hullArea;
			if (ratio < minConvexity || ratio >= maxConvexity)
				colors[label] = 0;
		}
	}

	Palete pal;
	unsigned int colorCount = 0;
	for (int i = 1; i < nccomps; i++)
	{
		CvLabel _label = i;
		double r, g, b;
		_HSV2RGB_((double)((colorCount * 77) % 360), .5, 1., r, g, b);
		colorCount++;
		pal[_label] = CV_RGB(r, g, b);
	}

	//过滤
	cv::parallel_for_(cv::Range(0, Y), [&](const cv::Range& range)->void {
		for (int y = range.start; y < range.end; y++) {
			uint8_t *ptrRow = image.ptr<uint8_t>(y);
			for (int x = 0; x < X; x++) {
				int label = labels.ptr<int>(y)[x];
				CV_Assert(0 <= label && label <= nccomps);
				ptrRow[x] = colors[label];
				if (colors[label]) {
					showMat.ptr<cv::Vec3b>(y)[x] = cv::Vec3b((uchar)pal[label].val[0], (uchar)pal[label].val[1], (uchar)pal[label].val[2]);
				}
			}
		}
	});

	EyemBinBlob blob;
	std::vector<EyemBinBlob> * tpResults = new std::vector<EyemBinBlob>();
	for (int i = 1; i < nccomps; i++) {
		if (colors[i]) {
			/*cv::rectangle(showMat, cv::Rect(stats.ptr<int>(i)[cv::CC_STAT_LEFT], stats.ptr<int>(i)[cv::CC_STAT_TOP],
				stats.ptr<int>(i)[cv::CC_STAT_WIDTH], stats.ptr<int>(i)[cv::CC_STAT_HEIGHT]), cv::Scalar(0, 0, 255));*/
			cv::drawMarker(showMat, cv::Point((int)centroids.ptr<double>(i)[0], (int)centroids.ptr<double>(i)[1]), cv::Scalar(255, 0, 0), cv::MARKER_CROSS, 6);

			double x1, y1, x2, y2;
			double lengthLine = MAX(stats.ptr<int>(i)[cv::CC_STAT_WIDTH], stats.ptr<int>(i)[cv::CC_STAT_HEIGHT]) / 2.;

			x1 = centroids.ptr<double>(i)[0] - lengthLine*cos(mOrientation.ptr<float>(i)[0]);
			y1 = centroids.ptr<double>(i)[1] - lengthLine*sin(mOrientation.ptr<float>(i)[0]);
			x2 = centroids.ptr<double>(i)[0] + lengthLine*cos(mOrientation.ptr<float>(i)[0]);
			y2 = centroids.ptr<double>(i)[1] + lengthLine*sin(mOrientation.ptr<float>(i)[0]);
			cv::line(showMat, cv::Point(int(x1), int(y1)), cv::Point(int(x2), int(y2)), cv::Scalar(0, 255, 0));
			blob.iXs = stats.ptr<int>(i)[cv::CC_STAT_LEFT];
			blob.iYs = stats.ptr<int>(i)[cv::CC_STAT_TOP];
			blob.iXe = blob.iYs + stats.ptr<int>(i)[cv::CC_STAT_WIDTH];
			blob.iYe = blob.iYs + stats.ptr<int>(i)[cv::CC_STAT_HEIGHT];
			blob.iWidth = stats.ptr<int>(i)[cv::CC_STAT_WIDTH];
			blob.iHeight = stats.ptr<int>(i)[cv::CC_STAT_HEIGHT];
			blob.dCenterX = centroids.ptr<double>(i)[0];
			blob.dCenterY = centroids.ptr<double>(i)[1];
			blob.iArea = stats.ptr<int>(i)[cv::CC_STAT_AREA];
			blob.dTheta = mOrientation.ptr<float>(i)[0];
			tpResults->push_back(blob);
		}
	}
	//<输出结果图像
	tpDstImg->iWidth = showMat.cols; tpDstImg->iHeight = showMat.rows; tpDstImg->iDepth = showMat.depth(); tpDstImg->iChannels = showMat.channels();
	//内存尺寸
	int _Size = tpDstImg->iWidth*tpDstImg->iHeight*tpDstImg->iChannels * sizeof(uint8_t);
	//分配初始化内存
	tpDstImg->vpImage = (uint8_t *)malloc(_Size);
	if (NULL == tpDstImg->vpImage)
		return FUNC_NOT_ENOUGH_MEM;
	memset(tpDstImg->vpImage, 0, _Size);
	//拷贝数据
	memcpy(tpDstImg->vpImage, showMat.data, _Size);
	//输出结果
	*ipNum = static_cast<int>(tpResults->size());
	*hObject = reinterpret_cast<IntPtr>(tpResults);
	*tpResult = tpResults->data();
	return FUNC_OK;
}

bool eyemBinFree(IntPtr hObject)
{
	std::vector<EyemBinBlob>  *tpResult = reinterpret_cast<std::vector<EyemBinBlob>*>(hObject);
	delete tpResult;
	tpResult = NULL;
	return true;
}