AutoThreshold (RenyiEntropy) C++

Java -> C++

 

Code

int RenyiEntropy(float* histo, int n_length) {
	// 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 = nullptr;
	double* P1 = nullptr; 
	double* P2 = nullptr;
	
	P1 = new double[n_length]; /* cumulative normalized histogram */
	P2 = new double[n_length];	
	norm_histo = new double[n_length]; /* normalized histogram */

	int total = 0;

	for (ih = 0; ih < n_length; ih++)
		total += histo[ih];

	for (ih = 0; ih < n_length; ih++)
		norm_histo[ih] = (double)histo[ih] / total;

	P1[0] = norm_histo[0];
	P2[0] = 1.0 - P1[0];

	for (ih = 1; ih < n_length; 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 < n_length; ih++) 
	{
		if (!(abs(P1[ih]) < 2.220446049250313E-16)) 
		{
			first_bin = ih;
			break;
		}
	}

	/* Determine the last non-zero bin */
	last_bin = n_length - 1;

	for (ih = n_length - 1; 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 (histo[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 < n_length; ih++) {
			if (histo[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;
		}
	}
	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 < n_length; 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 < n_length; 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;
		}
	}
	//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));

	if (norm_histo)
		delete[] norm_histo;

	if (P1)
		delete[] P1;

	if (P2)
		delete[] P2;

	return opt_threshold;
}

 

imageJ 관련 문서 : https://imagej.net/plugins/auto-threshold

github link : https://github.com/fiji/Auto_Threshold/blob/master/src/main/java/fiji/threshold/Auto_Threshold.java