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= 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;k2)
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; i0) 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; i10000) {
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 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; imaxCount) {
maxCount = data2[i];
mode = i;
}
}
int maxCount2 = 0;
for (int i = 0; imaxCount2) && (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) && (min0))
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 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
// 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= 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]= 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;
}
}
Ads help maintain this website.