ImageJ是一款非常优秀的、跨平台的开源软件，里面的源码是用Java编写，之所以用Java是因为Java可以跨平台。常用的图像处理是基于C++和OpenCV的，所以本人对ImageJ的自动二值化算法进行修改，使其能在C++上也能运行。
ImageJ的自动二值化算法大约有16种，之前有前辈修改之后在C#上运行。
读者也可以去下载本博客的完整测试代码:https://download.csdn.net/download/wfh2015/10424253
下面是ImageJ的原始Java代码，文件名为AutoThresholder.java，有兴趣的可以去下载：
package ij.process;
import ij.IJ;
import java.util.Arrays;
/** Autothresholding methods (limited to 256 bin histograms) from the Auto_Threshold plugin
(http://fiji.sc/Auto_Threshold) by G.Landini at bham dot ac dot uk). */
public class AutoThresholder {
private static String[] mStrings;
public enum Method {
Default,
Huang,
Intermodes,
IsoData,
IJ_IsoData,
Li,
MaxEntropy,
Mean,
MinError,
Minimum,
Moments,
Otsu,
Percentile,
RenyiEntropy,
Shanbhag,
Triangle,
Yen
};
public static String[] getMethods() {
if (mStrings==null) {
Method[] mVals = Method.values();
mStrings = new String[mVals.length];
for (int i=0; i<mVals.length; i++)
mStrings[i] = mVals[i].name();
}
return mStrings;
}
/** Calculates and returns a threshold using the specified
method and 256 bin histogram. */
public int getThreshold(Method method, int[] histogram) {
if (histogram==null)
throw new IllegalArgumentException("Histogram is null");
if (histogram.length!=256)
throw new IllegalArgumentException("Histogram length not 256");
int threshold = 0;
switch (method) {
case Default: threshold = defaultIsoData(histogram); break;
case IJ_IsoData: threshold = IJIsoData(histogram); break;
case Huang: threshold = Huang(histogram); break;
case Intermodes: threshold = Intermodes(histogram); break;
case IsoData: threshold = IsoData(histogram); break;
case Li: threshold = Li(histogram); break;
case MaxEntropy: threshold = MaxEntropy(histogram); break;
case Mean: threshold = Mean(histogram); break;
case MinError: threshold = MinErrorI(histogram); break;
case Minimum: threshold = Minimum(histogram); break;
case Moments: threshold = Moments(histogram); break;
case Otsu: threshold = Otsu(histogram); break;
case Percentile: threshold = Percentile(histogram); break;
case RenyiEntropy: threshold = RenyiEntropy(histogram); break;
case Shanbhag: threshold = Shanbhag(histogram); break;
case Triangle: threshold = Triangle(histogram); break;
case Yen: threshold = Yen(histogram); break;
}
if (threshold==-1) threshold = 0;
return threshold;
}
public int getThreshold(String mString, int[] histogram) {
// throws an exception if unknown argument
int index = mString.indexOf(" ");
if (index!=-1)
mString = mString.substring(0, index);
Method method = Method.valueOf(Method.class, mString);
return getThreshold(method, histogram);
}
int Huang(int [] data ) {
// 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 ( data[ih] != 0 ) {
first_bin = ih;
break;
}
}
/* Determine the last non-zero bin */
last_bin=255;
for (ih = 255; ih >= first_bin; ih-- ) {
if ( data[ih] != 0 ) {
last_bin = ih;
break;
}
}
term = 1.0 / ( double ) ( last_bin - first_bin );
double [] mu_0 = new double[256];
sum_pix = num_pix = 0;
for ( ih = first_bin; ih < 256; ih++ ){
sum_pix += (double)ih * data[ih];
num_pix += data[ih];
/* NUM_PIX cannot be zero ! */
mu_0[ih] = sum_pix / num_pix;
}
double [] mu_1 = new double[256];
sum_pix = num_pix = 0;
for ( ih = last_bin; ih > 0; ih-- ){
sum_pix += (double)ih * data[ih];
num_pix += data[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 = Double.MAX_VALUE;
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 * Math.abs ( ih - mu_0[it] ) );
if ( !((mu_x < 1e-06 ) || ( mu_x > 0.999999))) {
/* Equation (6) & (8) in Ref. 1 */
ent += data[ih] * ( -mu_x * Math.log ( mu_x ) - ( 1.0 - mu_x ) * Math.log ( 1.0 - mu_x ) );
}
}
for ( ih = it + 1; ih < 256; ih++ ) {
/* Equation (4) in Ref. 1 */
mu_x = 1.0 / ( 1.0 + term * Math.abs ( ih - mu_1[it] ) );
if ( !((mu_x < 1e-06 ) || ( mu_x > 0.999999))) {
/* Equation (6) & (8) in Ref. 1 */
ent += data[ih] * ( -mu_x * Math.log ( mu_x ) - ( 1.0 - mu_x ) * Math.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;
}
boolean bimodalTest(double [] y) {
int len=y.length;
boolean b = false;
int modes = 0;
for (int k=1;k<len-1;k++){
if (y[k-1] < y[k] && y[k+1] < y[k]) {
modes++;
if (modes>2)
return false;
}
}
if (modes == 2)
b = true;
return b;
}
int Intermodes(int[] data ) {
// J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
// Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
// 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.
//
// Assumes a bimodal histogram. The histogram needs is smoothed (using a
// running average of size 3, iteratively) until there are only two local maxima.
// j and k
// Threshold t is (j+k)/2.
// Images with histograms having extremely unequal peaks or a broad and
// flat valleys are unsuitable for this method.
int minbin=-1, maxbin=-1;
for (int i=0; i<data.length; i++)
if (data[i]>0) maxbin = i;
for (int i=data.length-1; i>=0; i--)
if (data[i]>0) minbin = i;
int length = (maxbin-minbin)+1;
double [] hist = new double[length];
for (int i=minbin; i<=maxbin; i++)
hist[i-minbin] = data[i];
int iter = 0;
int threshold=-1;
while (!bimodalTest(hist) ) {
//smooth with a 3 point running mean filter
double previous=0, current=0, next=hist[0];
for (int i=0; i<length-1; i++) {
previous = current;
current = next;
next = hist[i + 1];
hist[i] = (previous+current+next)/3;
}
hist[length-1] = (current+next)/3;
iter++;
if (iter>10000) {
threshold = -1;
IJ.log("Intermodes Threshold not found after 10000 iterations.");
return threshold;
}
}
// The threshold is the mean between the two peaks.
int tt=0;
for (int i=1; i<length - 1; i++) {
if (hist[i-1] < hist[i] && hist[i+1] < hist[i]){
tt += i;
//IJ.log("mode:" +i);
}
}
threshold = (int) Math.floor(tt/2.0);
return threshold+minbin;
}
int IsoData(int[] data ) {
// 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 (data[i] > 0){
g = i + 1;
break;
}
}
while (true){
l = 0;
totl = 0;
for (i = 0; i < g; i++) {
totl = totl + data[i];
l = l + (data[i] * i);
}
h = 0;
toth = 0;
for (i = g + 1; i < 256; i++){
toth += data[i];
h += ((double)data[i]*i);
}
if (totl > 0 && toth > 0){
l /= totl;
h /= toth;
if (g == (int) Math.round((l + h) / 2.0))
break;
}
g++;
if (g > 254)
return -1;
}
return g;
}
int defaultIsoData(int[] data) {
// This is the modified IsoData method used by the "Threshold" widget in "Default" mode
int n = data.length;
int[] data2 = new int[n];
int mode=0, maxCount=0;
for (int i=0; i<n; i++) {
int count = data[i];
data2[i] = data[i];
if (data2[i]>maxCount) {
maxCount = data2[i];
mode = i;
}
}
int maxCount2 = 0;
for (int i = 0; i<n; i++) {
if ((data2[i]>maxCount2) && (i!=mode))
maxCount2 = data2[i];
}
int hmax = maxCount;
if ((hmax>(maxCount2*2)) && (maxCount2!=0)) {
hmax = (int)(maxCount2 * 1.5);
data2[mode] = hmax;
}
return IJIsoData(data2);
}
int IJIsoData(int[] data) {
// This is the original ImageJ IsoData implementation, here for backward compatibility.
int level;
int maxValue = data.length - 1;
double result, sum1, sum2, sum3, sum4;
int count0 = data[0];
data[0] = 0; //set to zero so erased areas aren't included
int countMax = data[maxValue];
data[maxValue] = 0;
int min = 0;
while ((data[min]==0) && (min<maxValue))
min++;
int max = maxValue;
while ((data[max]==0) && (max>0))
max--;
if (min>=max) {
data[0]= count0; data[maxValue]=countMax;
level = data.length/2;
return level;
}
int movingIndex = min;
int inc = Math.max(max/40, 1);
do {
sum1=sum2=sum3=sum4=0.0;
for (int i=min; i<=movingIndex; i++) {
sum1 += (double)i*data[i];
sum2 += data[i];
}
for (int i=(movingIndex+1); i<=max; i++) {
sum3 += (double)i*data[i];
sum4 += data[i];
}
result = (sum1/sum2 + sum3/sum4)/2.0;
movingIndex++;
} while ((movingIndex+1)<=result && movingIndex<max-1);
data[0]= count0; data[maxValue]=countMax;
level = (int)Math.round(result);
return level;
}
int Li(int [] data ) {
// 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 += data[ih];
/* Calculate the mean gray-level */
mean = 0.0;
for (int ih = 0 + 1; ih < 256; ih++ ) //0 + 1?
mean += (double)ih * data[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 * data[ih];
num_back += data[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 * data[ih];
num_obj += data[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 ) / ( Math.log ( mean_back ) - Math.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 ( Math.abs ( new_thresh - old_thresh ) > tolerance );
return threshold;
}
int MaxEntropy(int [] data ) {
// 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 = new double[256]; /* normalized histogram */
double [] P1 = new double[256]; /* cumulative normalized histogram */
double [] P2 = new double[256];
double total =0;
for (ih = 0; ih < 256; ih++ )
total+=data[ih];
for (ih = 0; ih < 256; ih++ )
norm_histo[ih] = data[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 ( !(Math.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 ( !(Math.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 = Double.MIN_VALUE;
for ( it = first_bin; it <= last_bin; it++ ) {
/* Entropy of the background pixels */
ent_back = 0.0;
for ( ih = 0; ih <= it; ih++ ) {
if ( data[ih] !=0 ) {
ent_back -= ( norm_histo[ih] / P1[it] ) * Math.log ( norm_histo[ih] / P1[it] );
}
}
/* Entropy of the object pixels */
ent_obj = 0.0;
for ( ih = it + 1; ih < 256; ih++ ){
if (data[ih]!=0){
ent_obj -= ( norm_histo[ih] / P2[it] ) * Math.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;
}
int Mean(int [] data ) {
// 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+= data[i];
sum+=((double)i*data[i]);
}
threshold =(int) Math.floor(sum/tot);
return threshold;
}
int MinErrorI(int [] data ) {
// 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(data); //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(data, threshold)/A(data, threshold);
nu = (B(data, data.length - 1)-B(data, threshold))/(A(data, data.length - 1)-A(data, threshold));
p = A(data, threshold)/A(data, data.length - 1);
q = (A(data, data.length - 1)-A(data, threshold)) / A(data, data.length - 1);
sigma2 = C(data, threshold)/A(data, threshold)-(mu*mu);
tau2 = (C(data, data.length - 1)-C(data, threshold)) / (A(data, data.length - 1)-A(data, 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 + Math.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) {
IJ.log("MinError(I): not converging.");
return threshold;
}
//The updated threshold is the integer part of the solution of the quadratic equation.
Tprev = threshold;
temp = (w1+Math.sqrt(sqterm))/w0;
if (Double.isNaN(temp))
threshold = Tprev;
else
threshold =(int) Math.floor(temp);
}
return threshold;
}
private double A(int[] y, int j) {
if (j>=y.length) j=y.length-1;
double x = 0;
for (int i=0;i<=j;i++)
x+=y[i];
return x;
}
private double B(int[] y, int j) {
if (j>=y.length) j=y.length-1;
double x = 0;
for (int i=0;i<=j;i++)
x+=i*y[i];
return x;
}
private double C(int[] y, int j) {
if (j>=y.length) j=y.length-1;
double x = 0;
for (int i=0;i<=j;i++)
x+=i*i*y[i];
return x;
}
int Minimum(int [] data ) {
// J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
// Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
// 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.
//
// Assumes a bimodal histogram. The histogram needs is smoothed (using a
// running average of size 3, iteratively) until there are only two local maxima.
// Threshold t is such that yt-1 > yt <= yt+1.
// Images with histograms having extremely unequal peaks or a broad and
// flat valleys are unsuitable for this method.
int iter =0;
int threshold = -1;
double [] iHisto = new double [256];
for (int i=0; i<256; i++)
iHisto[i]=(double) data[i];
double [] tHisto = new double[iHisto.length] ;
while (!bimodalTest(iHisto) ) {
//smooth with a 3 point running mean filter
for (int i=1; i<255; i++)
tHisto[i]= (iHisto[i-1] + iHisto[i] +iHisto[i+1])/3;
tHisto[0] = (iHisto[0]+iHisto[1])/3; //0 outside
tHisto[255] = (iHisto[254]+iHisto[255])/3; //0 outside
System.arraycopy(tHisto, 0, iHisto, 0, iHisto.length) ;
iter++;
if (iter>10000) {
threshold = -1;
IJ.log("Minimum: threshold not found after 10000 iterations.");
return threshold;
}
}
// The threshold is the minimum between the two peaks.
for (int i=1; i<255; i++) {
if (iHisto[i-1] > iHisto[i] && iHisto[i+1] >= iHisto[i]) {
threshold = i;
break;
}
}
return threshold;
}
int Moments(int [] data ) {
// 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 = new double [256];
for (int i=0; i<256; i++)
total+=data[i];
for (int i=0; i<256; i++)
histo[i]=(double)(data[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 - Math.sqrt ( c1 * c1 - 4.0 * c0 ) );
z1 = 0.5 * ( -c1 + Math.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;
}
int Otsu(int [] data ) {
// 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 * data[k]; // Total histogram intensity
N += data[k]; // Total number of data points
}
Sk = 0;
N1 = data[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 * data[k];
N1 += data[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;
}
}
// kStar += 1; // Use QTI convention that intensity -> 1 if intensity >= k
// (the algorithm was developed for I-> 1 if I <= k.)
return kStar;
}
int Percentile(int [] data ) {
// 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 = new double [256];
for (int i=0; i<256; i++)
avec[i]=0.0;
double total =partialSum(data, 255);
double temp = 1.0;
for (int i=0; i<256; i++){
avec[i]=Math.abs((partialSum(data, i)/total)-ptile);
//IJ.log("Ptile["+i+"]:"+ avec[i]);
if (avec[i]<temp) {
temp = avec[i];
threshold = i;
}
}
return threshold;
}
double partialSum(int [] y, int j) {
double x = 0;
for (int i=0;i<=j;i++)
x+=y[i];
return x;
}
int RenyiEntropy(int [] data ) {
// 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 = new double[256]; /* normalized histogram */
double [] P1 = new double[256]; /* cumulative normalized histogram */
double [] P2 = new double[256];
double total =0;
for (ih = 0; ih < 256; ih++ )
total+=data[ih];
for (ih = 0; ih < 256; ih++ )
norm_histo[ih] = data[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 ( !(Math.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 ( !(Math.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 ( data[ih] !=0 ) {
ent_back -= ( norm_histo[ih] / P1[it] ) * Math.log ( norm_histo[ih] / P1[it] );
}
}
/* Entropy of the object pixels */
ent_obj = 0.0;
for ( ih = it + 1; ih < 256; ih++ ){
if (data[ih]!=0){
ent_obj -= ( norm_histo[ih] / P2[it] ) * Math.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 += Math.sqrt ( norm_histo[ih] / P1[it] );
/* Entropy of the object pixels */
ent_obj = 0.0;
for ( ih = it + 1; ih < 256; ih++ )
ent_obj += Math.sqrt ( norm_histo[ih] / P2[it] );
/* Total entropy */
tot_ent = term * ( ( ent_back * ent_obj ) > 0.0 ? Math.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 ? Math.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 ( Math.abs ( t_star1 - t_star2 ) <= 5 ) {
if ( Math.abs ( t_star2 - t_star3 ) <= 5 ) {
beta1 = 1;
beta2 = 2;
beta3 = 1;
}
else {
beta1 = 0;
beta2 = 1;
beta3 = 3;
}
}
else {
if ( Math.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 ));
return opt_threshold;
}
int Shanbhag(int [] data ) {
// 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 = new double[256]; /* normalized histogram */
double [] P1 = new double[256]; /* cumulative normalized histogram */
double [] P2 = new double[256];
double total =0;
for (ih = 0; ih < 256; ih++ )
total+=data[ih];
for (ih = 0; ih < 256; ih++ )
norm_histo[ih] = data[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 ( !(Math.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 ( !(Math.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 = Double.MAX_VALUE;
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] * Math.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] * Math.log ( 1.0 - term * P2[ih] );
}
ent_obj *= term;
/* Total entropy */
tot_ent = Math.abs ( ent_back - ent_obj );
if ( tot_ent < min_ent ) {
min_ent = tot_ent;
threshold = it;
}
}
return threshold;
}
int Triangle(int [] data ) {
// 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 < data.length; i++) {
if (data[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 (data[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 (data[i] >dmax) {
max=i;
dmax=data[i];
}
}
// find which is the furthest side
//IJ.log(""+min+" "+max+" "+min2);
boolean inverted = false;
if ((max-min)<(min2-max)){
// reverse the histogram
//IJ.log("Reversing 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 = data[left];
data[left] = data[right];
data[right] = temp;
// move the bounds toward the center
left++;
right--;
}
min=255-min2;
max=255-max;
}
if (min == max){
//IJ.log("Triangle: 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 = data[max]; //-min; // data[min]; // lowest value bmin = (p=0)% in the image
ny = min - max;
d = Math.sqrt(nx * nx + ny * ny);
nx /= d;
ny /= d;
d = nx * min + ny * data[min];
// find split point
int split = min;
double splitDistance = 0;
for (int i = min + 1; i <= max; i++) {
double newDistance = nx * i + ny * data[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 = data[left];
data[left] = data[right];
data[right] = temp;
left++;
right--;
}
return (255-split);
}
else
return split;
}
int Yen(int [] data ) {
// 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 = new double[256]; /* normalized histogram */
double [] P1 = new double[256]; /* cumulative normalized histogram */
double [] P1_sq = new double[256];
double [] P2_sq = new double[256];
double total =0;
for (ih = 0; ih < 256; ih++ )
total+=data[ih];
for (ih = 0; ih < 256; ih++ )
norm_histo[ih] = data[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 = Double.MIN_VALUE;
for ( it = 0; it < 256; it++ ) {
crit = -1.0 * (( P1_sq[it] * P2_sq[it] )> 0.0? Math.log( P1_sq[it] * P2_sq[it]):0.0) + 2 * ( ( P1[it] * ( 1.0 - P1[it] ) )>0.0? Math.log( P1[it] * ( 1.0 - P1[it] ) ): 0.0);
if ( crit > max_crit ) {
max_crit = crit;
threshold = it;
}
}
return threshold;
}
}
本人修改后的源码头文件，文件名称为auto_thresholder.h
/** @file
* @brief 自动计算单通道图像二值化阈值
* @author ImageJ fh
*/
#ifndef AUTO_THRESHOLDER
#define AUTO_THRESHOLDER
#include <vector>
#include "opencv2/core/core.hpp"
/** @brief 各种自动计算阈值方法
输入的图片必须是CV_8UC1类型cv::Mat
具体使用方法可参考下列代码:
@code
std::string path = "a.bmp";
cv::Mat img = cv::imread(path, cv::IMREAD_GRAYSCALE);
AutoThresholder::CalcHist(img, hist);
int threshValue = AutoThresholder::Triangle(hist);
cv::Mat bin;
cv::threshold(img, bin, threshValue, 255, cv::THRESH_BINARY);
@endcode
*/
class AutoThresholder
{
public:
/** @brief 统计每个灰度级及其对应的像素个数
@param src CV_8UC1类型的图像
@param data 所得到的结果, 索引为灰度级, 对应值为像素个数. 内有256个元素
*/
static void CalcHist(const cv::Mat& src, std::vector<int>& data);
public:
static int DefaultIsoData(const std::vector<int>& data);
/** @brief original ImageJ IsoData implementation
函数有可能会修改data
*/
static int IJIsoData(std::vector<int>& data);
/** @brief 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
*/
static int Huang(const std::vector<int>& data);
/** @brief Implements Intermodes thresholding method
J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
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.
Assumes a bimodal histogram. The histogram needs is smoothed (using a
running average of size 3, iteratively) until there are only two local maxima.
j and k
Threshold t is (j+k)/2.
Images with histograms having extremely unequal peaks or a broad and
flat valleys are unsuitable for this method.
*/
static int Intermodes(const std::vector<int>& data);
/** @brief Implements IsoData thresholding method
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 iteratingL = 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
*/
static int IsoData(const std::vector<int>& data);
/** @brief 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
*/
static int Li(const std::vector<int>& data);
/** @brief 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
*/
static int MaxEntropy(const std::vector<int>& data);
/** @brief Implements Mean thresholding method
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
*/
static int Mean(const std::vector<int>& data);
/** @brief Implements Minimum error thresholding thresholding method
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.
*/
static int MinErrorI(const std::vector<int>& data);
/**
J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
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.
Assumes a bimodal histogram. The histogram needs is smoothed (using a
running average of size 3, iteratively) until there are only two local maxima.
Threshold t is such that yt-1 > yt <= yt+1.
Images with histograms having extremely unequal peaks or a broad and
flat valleys are unsuitable for this method.
*/
static int Minimum(const std::vector<int>& data);
/**
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
*/
static int Moments(const std::vector<int>& data);
/**
Otsu's threshold algorithm
C++ code by Jordan Bevik <Jordan.Bevic@qtiworld.com>
ported to ImageJ plugin by G.Landini
*/
static int Otsu(const std::vector<int>& data);
/**
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.
*/
static int Percentile(const std::vector<int>& data);
/**
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
*/
static int RenyiEntropy(const std::vector<int>& data);
/**
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
*/
static int Shanbhag(const std::vector<int>& data);
/**
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
函数有可能会修改data
*/
static int Triangle(std::vector<int>& data);
/** @brief 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 routin
*/
static int Yen(const std::vector<int>& data);
};
#endif // !AUTO_THRESHOLDER
具体算法实现cpp为auto_thresholder.cpp
/** @file
* @brief 自动计算单通道图像二值化阈值
* @author ImageJ fh
* 由于VC编译器Max, Min函数可能会出现问题, 所以有关求最值函数全部自定义在cpp文件中
*/
#include "auto_thresholder.h"
/*****************************************常用函数*********************************************/
static int CalcMaxValue(int a, int b)
{
return (a > b) ? a : b;
}
static double CalcMaxValue(double a, double b)
{
return (a > b) ? a : b;
}
static int CalcMinValue(int a, int b)
{
return (a < b) ? a : b;
}
static double CalcMinValue(double a, double b)
{
return (a < b) ? a : b;
}
void AutoThresholder::CalcHist(const cv::Mat& src, std::vector<int>& data)
{
CV_Assert(CV_8UC1 == src.type());
data.clear();
// 初始化数组
const int GRAY_LEVEL = 256;
for (int i = 0; i < GRAY_LEVEL; ++i)
{
data.push_back(0);
}
// 统计每个灰度级及其对应的像素个数
for (int y = 0; y < src.rows; ++y)
{
for (int x = 0; x < src.cols; ++x)
{
int value = src.at<uchar>(y, x);
data.at(value) += 1;
}
}
}
int AutoThresholder::DefaultIsoData(const std::vector<int>& data)
{
// This is the modified IsoData method used by the "Threshold" widget in "Default" mode
int n = data.size();
std::vector<int> data2(n, 0);
int mode = 0, maxCount = 0;
for (int i = 0; i < n; i++)
{
int count = data[i];
data2[i] = data[i];
if (data2[i] > maxCount)
{
maxCount = data2[i];
mode = i;
}
}
int maxCount2 = 0;
for (int i = 0; i < n; i++)
{
if ((data2[i] > maxCount2) && (i != mode))
maxCount2 = data2[i];
}
int hmax = maxCount;
if ((hmax > (maxCount2 * 2)) && (maxCount2 != 0))
{
hmax = (int)(maxCount2 * 1.5);
data2[mode] = hmax;
}
return IJIsoData(data2);
}
int AutoThresholder::IJIsoData(std::vector<int>& data)
{
// This is the original ImageJ IsoData implementation, here for backward compatibility.
int level;
int maxValue = data.size() - 1;
double result, sum1, sum2, sum3, sum4;
int count0 = data[0];
data[0] = 0; //set to zero so erased areas aren't included
int countMax = data[maxValue];
data[maxValue] = 0;
int min = 0;
while ((data[min] == 0) && (min < maxValue))
min++;
int max = maxValue;
while ((data[max] == 0) && (max > 0))
max--;
if (min >= max)
{
data[0] = count0; data[maxValue] = countMax;
level = data.size() / 2;
return level;
}
int movingIndex = min;
int inc = CalcMaxValue(max / 40, 1);
do
{
sum1 = sum2 = sum3 = sum4 = 0.0;
for (int i = min; i <= movingIndex; i++)
{
sum1 += (double)i*data[i];
sum2 += data[i];
}
for (int i = (movingIndex + 1); i <= max; i++)
{
sum3 += (double)i*data[i];
sum4 += data[i];
}
result = (sum1 / sum2 + sum3 / sum4) / 2.0;
movingIndex++;
} while ((movingIndex + 1) <= result && movingIndex < max - 1);
data[0] = count0; data[maxValue] = countMax;
level = (int)std::round(result);
return level;
}
int AutoThresholder::Huang(const std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
int threshold = -1;
int ih = 0;
int it = 0;
int first_bin = 0;
int last_bin = 0;
double sum_pix = 0.0;
double num_pix = 0.0;
double term = 0.0;
double ent = 0.0; // entropy
double min_ent = 0.0; // min entropy
double mu_x = 0.0;
/* Determine the first non-zero bin */
first_bin = 0;
for (ih = 0; ih < GRAY_LEVELS; ih++)
{
if (data[ih] != 0)
{
first_bin = ih;
break;
}
}
/* Determine the last non-zero bin */
last_bin = 255;
for (ih = 255; ih >= first_bin; ih--)
{
if (data[ih] != 0)
{
last_bin = ih;
break;
}
}
term = 1.0 / (double)(last_bin - first_bin);
std::vector<double> mu_0(GRAY_LEVELS, 0.0);
sum_pix = num_pix = 0;
for (ih = first_bin; ih < GRAY_LEVELS; ih++)
{
sum_pix += (double)ih * data[ih];
num_pix += data[ih];
/* NUM_PIX cannot be zero ! */
mu_0[ih] = sum_pix / num_pix;
}
//double[] mu_1 = new double[256];
std::vector<double> mu_1(GRAY_LEVELS, 0.0);
sum_pix = num_pix = 0;
for (ih = last_bin; ih > 0; ih--)
{
sum_pix += (double)ih * data[ih];
num_pix += data[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 < GRAY_LEVELS; it++)
{
ent = 0.0;
for (ih = 0; ih <= it; ih++)
{
/* Equation (4) in Ref. 1 */
mu_x = 1.0 / (1.0 + term * std::abs(ih - mu_0[it]));
if (!((mu_x < 1e-06) || (mu_x > 0.999999)))
{
/* Equation (6) & (8) in Ref. 1 */
ent += data[ih] * (-mu_x * std::log(mu_x) - (1.0 - mu_x) * std::log(1.0 - mu_x));
}
}
for (ih = it + 1; ih < GRAY_LEVELS; ih++)
{
/* Equation (4) in Ref. 1 */
mu_x = 1.0 / (1.0 + term * std::abs(ih - mu_1[it]));
if (!((mu_x < 1e-06) || (mu_x > 0.999999)))
{
/* Equation (6) & (8) in Ref. 1 */
ent += data[ih] * (-mu_x * std::log(mu_x) - (1.0 - mu_x) * std::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 bool BimodalTest(std::vector<double>& y)
{
int len = y.size();
bool b = false;
int modes = 0;
for (int k = 1; k < len - 1; k++)
{
if (y[k - 1] < y[k] && y[k + 1] < y[k])
{
modes++;
if (modes > 2) return false;
}
}
if (modes == 2) b = true;
return b;
}
int AutoThresholder::Intermodes(const std::vector<int>& data)
{
int maxbin = -1;
for (int i = 0; i < data.size(); i++)
{
if (data[i] > 0) maxbin = i;
}
int minbin = -1;
for (int i = data.size() - 1; i >= 0; i--)
{
if (data[i] > 0) minbin = i;
}
int length = (maxbin - minbin) + 1;
std::vector<double> hist(data.size(), 0);
for (int i = minbin; i <= maxbin; i++)
{
hist[i - minbin] = data[i];
}
int iter = 0;
int threshold = -1;
while (!BimodalTest(hist))
{
//smooth with a 3 point running mean filter
double previous = 0;
double current = 0;
double next = hist[0];
for (int i = 0; i < length - 1; i++)
{
previous = current;
current = next;
next = hist[i + 1];
hist[i] = (previous + current + next) / 3;
}
hist[length - 1] = (current + next) / 3;
iter++;
if (iter > 10000)
{
threshold = -1;
return threshold;
}
}
// The threshold is the mean between the two peaks.
int tt = 0;
for (int i = 1; i < length - 1; i++)
{
if (hist[i - 1] < hist[i] && hist[i + 1] < hist[i])
{
tt += i;
}
}
threshold = (int)std::floor(tt / 2.0);
return threshold + minbin;
}
int AutoThresholder::IsoData(const std::vector<int>& data)
{
int i = 0;
int l = 0;
int totl = 0;
int g = 0;
double toth = 0.0;
double h = 0.0;
const int GRAY_LEVELS = 256;
for (i = 1; i < GRAY_LEVELS; i++)
{
if (data[i] > 0)
{
g = i + 1;
break;
}
}
while (true)
{
l = 0;
totl = 0;
for (i = 0; i < g; i++)
{
totl = totl + data[i];
l = l + (data[i] * i);
}
h = 0;
toth = 0;
for (i = g + 1; i < GRAY_LEVELS; i++)
{
toth += data[i];
h += ((double)data[i] * i);
}
if (totl > 0 && toth > 0)
{
l /= totl;
h /= toth;
if (g == (int)std::round((l + h) / 2.0))
break;
}
g++;
if (g > 254) return -1;
}
return g;
}
int AutoThresholder::Li(const std::vector<int>& data)
{
int threshold = 0;
double num_pixels = 0;
double sum_back = 0; /* sum of the background pixels at a given threshold */
double sum_obj = 0; /* sum of the object pixels at a given threshold */
double num_back = 0; /* number of background pixels at a given threshold */
double num_obj = 0; /* number of object pixels at a given threshold */
double old_thresh = 0;
double new_thresh = 0;
double mean_back = 0; /* mean of the background pixels at a given threshold */
double mean_obj = 0; /* mean of the object pixels at a given threshold */
double mean = 0; /* mean gray-level in the image */
double tolerance = 0; /* threshold tolerance */
double temp = 0;
tolerance = 0.5;
num_pixels = 0;
for (int ih = 0; ih < 256; ih++)
num_pixels += data[ih];
/* Calculate the mean gray-level */
mean = 0.0;
for (int ih = 0 + 1; ih < 256; ih++) //0 + 1?
mean += (double)ih * data[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 * data[ih];
num_back += data[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 * data[ih];
num_obj += data[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) / (std::log(mean_back) - std::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 (std::abs(new_thresh - old_thresh) > tolerance);
return threshold;
}
int AutoThresholder::MaxEntropy(const std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
int threshold = -1;
int ih = 0;
int it = 0;
int first_bin = 0;
int last_bin = 0;
double tot_ent = 0.0; /* total entropy */
double max_ent = 0.0; /* max entropy */
double ent_back = 0.0; /* entropy of the background pixels at a given threshold */
double ent_obj = 0.0; /* entropy of the object pixels at a given threshold */
std::vector<double> norm_histo(GRAY_LEVELS, 0.0);/* normalized histogram */
std::vector<double> P1(GRAY_LEVELS, 0.0);/* cumulative normalized histogram */
std::vector<double> P2(GRAY_LEVELS, 0.0);
double total = 0;
for (ih = 0; ih < GRAY_LEVELS; ih++)
total += data[ih];
for (ih = 0; ih < GRAY_LEVELS; ih++)
norm_histo[ih] = data[ih] / total;
P1[0] = norm_histo[0];
P2[0] = 1.0 - P1[0];
for (ih = 1; ih < GRAY_LEVELS; 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 < GRAY_LEVELS; ih++)
{
if (!(std::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 (!(std::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 (data[ih] != 0)
{
ent_back -= (norm_histo[ih] / P1[it]) * std::log(norm_histo[ih] / P1[it]);
}
}
/* Entropy of the object pixels */
ent_obj = 0.0;
for (ih = it + 1; ih < 256; ih++)
{
if (data[ih] != 0)
{
ent_obj -= (norm_histo[ih] / P2[it]) * std::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;
}
}
return threshold;
}
int AutoThresholder::Mean(const std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
int threshold = -1;
double tot = 0, sum = 0;
for (int i = 0; i < GRAY_LEVELS; i++)
{
tot += data[i];
sum += ((double)i*data[i]);
}
threshold = (int)std::floor(sum / tot);
return threshold;
}
/***************************************MinErrorI***************************************************/
static double A(const std::vector<int>& y, int j)
{
if (j >= y.size()) j = y.size() - 1;
double x = 0;
for (int i = 0; i <= j; i++)
x += y[i];
return x;
}
static double B(const std::vector<int>& y, int j)
{
if (j >= y.size()) j = y.size() - 1;
double x = 0;
for (int i = 0; i <= j; i++)
x += i*y[i];
return x;
}
static double C(const std::vector<int>& y, int j)
{
if (j >= y.size()) j = y.size() - 1;
double x = 0;
for (int i = 0; i <= j; i++)
x += i*i*y[i];
return x;
}
int AutoThresholder::MinErrorI(const std::vector<int>& data)
{
int threshold = Mean(data); //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;
while (threshold != Tprev)
{
//Calculate some statistics.
mu = B(data, threshold) / A(data, threshold);
nu = (B(data, data.size() - 1) - B(data, threshold)) / (A(data, data.size() - 1) - A(data, threshold));
p = A(data, threshold) / A(data, data.size() - 1);
q = (A(data, data.size() - 1) - A(data, threshold)) / A(data, data.size() - 1);
sigma2 = C(data, threshold) / A(data, threshold) - (mu*mu);
tau2 = (C(data, data.size() - 1) - C(data, threshold)) / (A(data, data.size() - 1) - A(data, 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 + std::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)
{
return threshold;
}
//The updated threshold is the integer part of the solution of the quadratic equation.
Tprev = threshold;
temp = (w1 + std::sqrt(sqterm)) / w0;
if (std::isnan(temp))
threshold = Tprev;
else
threshold = (int)std::floor(temp);
}
return threshold;
}
int AutoThresholder::Minimum(const std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
int iter = 0;
int threshold = -1;
std::vector<double> iHisto(GRAY_LEVELS, 0.0);
for (int i = 0; i < GRAY_LEVELS; i++)
{
iHisto[i] = (double)data[i];
}
std::vector<double> tHisto(iHisto.size(), 0.0);
while (!BimodalTest(iHisto))
{
//smooth with a 3 point running mean filter
for (int i = 1; i < 255; i++)
{
tHisto[i] = (iHisto[i - 1] + iHisto[i] + iHisto[i + 1]) / 3;
}
tHisto[0] = (iHisto[0] + iHisto[1]) / 3; //0 outside
tHisto[255] = (iHisto[254] + iHisto[255]) / 3; //0 outside
std::copy(tHisto.begin(), tHisto.end(), iHisto.begin());
iter++;
if (iter > 10000)
{
threshold = -1;
return threshold;
}
}
// The threshold is the minimum between the two peaks.
for (int i = 1; i < 255; i++)
{
if (iHisto[i - 1] > iHisto[i] && iHisto[i + 1] >= iHisto[i])
{
threshold = i;
break;
}
}
return threshold;
}
int AutoThresholder::Moments(const std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
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 = new double[256];
std::vector<double> histo(GRAY_LEVELS, 0.0);
for (int i = 0; i < GRAY_LEVELS; i++)
{
total += data[i];
}
for (int i = 0; i < GRAY_LEVELS; i++)
{
histo[i] = (double)(data[i] / total); //normalised histogram
}
/* Calculate the first, second, and third order moments */
for (int i = 0; i < GRAY_LEVELS; 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 - std::sqrt(c1 * c1 - 4.0 * c0));
z1 = 0.5 * (-c1 + std::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 < GRAY_LEVELS; i++)
{
sum += histo[i];
if (sum > p0)
{
threshold = i;
break;
}
}
return threshold;
}
int AutoThresholder::Otsu(const std::vector<int>& data)
{
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 * data[k]; // Total histogram intensity
N += data[k]; // Total number of data points
}
Sk = 0;
N1 = data[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 * data[k];
N1 += data[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;
}
}
// kStar += 1; // Use QTI convention that intensity -> 1 if intensity >= k
// (the algorithm was developed for I-> 1 if I <= k.)
return kStar;
}
static double PartialSum(const std::vector<int>& y, int j)
{
double x = 0;
for (int i = 0; i <= j; i++)
x += y[i];
return x;
}
int AutoThresholder::Percentile(const std::vector<int>& data)
{
int iter = 0;
int threshold = -1;
double ptile = 0.5; // default fraction of foreground pixels
std::vector<double> avec(256, 0.0);
for (int i = 0; i < 256; i++)
{
avec[i] = 0.0;
}
double total = PartialSum(data, 255);
double temp = 1.0;
for (int i = 0; i < 256; i++)
{
avec[i] = std::abs((PartialSum(data, i) / total) - ptile);
if (avec[i] < temp)
{
temp = avec[i];
threshold = i;
}
}
return threshold;
}
int AutoThresholder::RenyiEntropy(const std::vector<int>& data)
{
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;
std::vector<double> norm_histo(256, 0.0);/* normalized histogram */
std::vector<double> P1(256, 0.0);/* cumulative normalized histogram */
std::vector<double> P2(256, 0.0);
double total = 0;
for (ih = 0; ih < 256; ih++)
total += data[ih];
for (ih = 0; ih < 256; ih++)
norm_histo[ih] = data[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 (!(std::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 (!(std::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 (data[ih] != 0) {
ent_back -= (norm_histo[ih] / P1[it]) * std::log(norm_histo[ih] / P1[it]);
}
}
/* Entropy of the object pixels */
ent_obj = 0.0;
for (ih = it + 1; ih < 256; ih++)
{
if (data[ih] != 0) {
ent_obj -= (norm_histo[ih] / P2[it]) * std::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 += std::sqrt(norm_histo[ih] / P1[it]);
/* Entropy of the object pixels */
ent_obj = 0.0;
for (ih = it + 1; ih < 256; ih++)
ent_obj += std::sqrt(norm_histo[ih] / P2[it]);
/* Total entropy */
tot_ent = term * ((ent_back * ent_obj) > 0.0 ? std::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 ? std::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 (std::abs(t_star1 - t_star2) <= 5)
{
if (std::abs(t_star2 - t_star3) <= 5)
{
beta1 = 1;
beta2 = 2;
beta3 = 1;
}
else
{
beta1 = 0;
beta2 = 1;
beta3 = 3;
}
}
else
{
if (std::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;
}
int AutoThresholder::Shanbhag(const std::vector<int>& data)
{
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 */
std::vector<double> norm_histo(256, 0.0);/* normalized histogram */
std::vector<double> P1(256, 0.0);/* cumulative normalized histogram */
std::vector<double> P2(256, 0.0);
double total = 0;
for (ih = 0; ih < 256; ih++)
total += data[ih];
for (ih = 0; ih < 256; ih++)
norm_histo[ih] = data[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 (!(std::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 (!(std::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] * std::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] * std::log(1.0 - term * P2[ih]);
}
ent_obj *= term;
/* Total entropy */
tot_ent = std::abs(ent_back - ent_obj);
if (tot_ent < min_ent)
{
min_ent = tot_ent;
threshold = it;
}
}
return threshold;
}
int AutoThresholder::Triangle(std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
// find min and max
int min = 0, dmax = 0, max = 0, min2 = 0;
for (int i = 0; i < data.size(); i++)
{
if (data[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 (data[i] > 0)
{
min2 = i;
break;
}
}
if (min2 < 255) min2++; // line to the (p==0) point, not to data[min]
for (int i = 0; i < GRAY_LEVELS; i++)
{
if (data[i] > dmax)
{
max = i;
dmax = data[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 = data[left];
data[left] = data[right];
data[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 = data[max]; //-min; // data[min]; // lowest value bmin = (p=0)% in the image
ny = min - max;
d = std::sqrt(nx * nx + ny * ny);
nx /= d;
ny /= d;
d = nx * min + ny * data[min];
// find split point
int split = min;
double splitDistance = 0;
for (int i = min + 1; i <= max; i++)
{
double newDistance = nx * i + ny * data[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 = data[left];
data[left] = data[right];
data[right] = temp;
left++;
right--;
}
return (255 - split);
}
else
return split;
}
int AutoThresholder::Yen(const std::vector<int>& data)
{
const int GRAY_LEVELS = 256;
int threshold = 0;
int ih = 0;
int it=0;
double crit = 0.0;
double max_crit = 0.0;
std::vector<double> norm_histo(GRAY_LEVELS, 0.0);/* normalized histogram */
std::vector<double> P1(GRAY_LEVELS, 0.0);
std::vector<double> P1_sq(GRAY_LEVELS, 0.0);
std::vector<double> P2_sq(GRAY_LEVELS, 0.0);
double total = 0;
for (ih = 0; ih < GRAY_LEVELS; ih++)
{
total += data[ih];
}
for (ih = 0; ih < GRAY_LEVELS; ih++)
{
norm_histo[ih] = data[ih] / total;
}
P1[0] = norm_histo[0];
for (ih = 1; ih < GRAY_LEVELS; ih++)
{
P1[ih] = P1[ih - 1] + norm_histo[ih];
}
P1_sq[0] = norm_histo[0] * norm_histo[0];
for (ih = 1; ih < GRAY_LEVELS; 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 < GRAY_LEVELS; it++)
{
crit = -1.0 * ((P1_sq[it] * P2_sq[it]) > 0.0 ? std::log(P1_sq[it] * P2_sq[it]) : 0.0) + 2 * ((P1[it] * (1.0 - P1[it])) > 0.0 ? std::log(P1[it] * (1.0 - P1[it])) : 0.0);
if (crit > max_crit)
{
max_crit = crit;
threshold = it;
}
}
return threshold;
}
2656
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
