AutoThreshold (MinErrorI) C++

Java -> C++

 

Code

#include <cmath>

int Mean(float* histo, int n_length) 
{
	// 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;
	long tot = 0, sum = 0;
	
	for (int i = 0; i < n_length; i++) 
	{
		tot += histo[i];
		sum += ((long)i * (long)histo[i]);
	}

	threshold = (int)floor(sum / tot);

	return threshold;
}

inline double A(float* y, int j)
{
	double x = 0;

	for (int i = 0; i <= j; i++)
		x += y[i];

	return x;
}

inline double B(float* y, int j)
{
	double x = 0;

	for (int i = 0; i <= j; i++)
		x += i * y[i];

	return x;
}

inline double C(float* y, int j) 
{
	double x = 0;

	for (int i = 0; i <= j; i++)
		x += i * i * y[i];

	return x;
}


int MinErrorI(float* histo, int n_length) 
{
	// 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(histo, n_length); //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(histo, threshold) / A(histo, threshold);
		nu = (B(histo, n_length - 1) - B(histo, threshold)) / (A(histo, n_length - 1) - A(histo, threshold));
		p = A(histo, threshold) / A(histo, n_length - 1);
		q = (A(histo, n_length - 1) - A(histo, threshold)) / A(histo, n_length - 1);
		sigma2 = C(histo, threshold) / A(histo, threshold) - (mu * mu);
		tau2 = (C(histo, n_length - 1) - C(histo, threshold)) / (A(histo, n_length - 1) - A(histo, 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 + 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) {
			std::cout<< ("MinError(I): not converging. Try \'Ignore black/white\' options");
			return threshold;
		}

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

		if ((isnan(temp))) {
			std::cout<<("MinError(I): NaN, not converging. Try \'Ignore black/white\' options");
			threshold = Tprev;
		}
		else
			threshold = (int)floor(temp);
		//IJ.log("Iter: "+ counter+++"  t:"+threshold);
	}
	return 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