http://www.cnblogs.com/pixels/archive/2011/01/14/1935600.html
從ImageJ中挑選了三種閾值的計算方法,下圖是運行結果。可以看出均值的方法獲得的二值結果較差。
下面是源代碼(改自ImageJ中的Java代碼,其中給出的參考文獻,有興趣的可以看一下)
代碼
#include "stdafx.h"
#include "cv.h"
#include "cxcore.h"
#include "highgui.h"
#pragma comment(lib, "cv.lib")
#pragma comment(lib, "cxcore.lib")
#pragma comment(lib, "highgui.lib")
int hist_size = 256;
float range_0[] = {0, 256};
float* ranges[] = {range_0};
CvHistogram *hist;
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;
int sum_pix;
int 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[256]={0.0};
sum_pix = num_pix = 0;
for ( ih = first_bin; ih < 256; ih++ )
{
sum_pix += ih * data[ih];
num_pix += data[ih];
/* NUM_PIX cannot be zero ! */
mu_0[ih] = sum_pix / ( double ) num_pix;
}
double mu_1[256]={0.0};
sum_pix = num_pix = 0;
for ( ih = last_bin; ih > 0; ih-- )
{
sum_pix += 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 = 65535;
for ( it = 0; it < 256; it++ )
{
ent = 0.0;
for ( ih = 0; ih <= it; ih++ )
{
/* Equation (4) in Ref. 1 */
mu_x = 1.0 / ( 1.0 + term * abs ( ih - mu_0[it] ) );
if ( !((mu_x < 1e-06 ) || ( mu_x > 0.999999)))
{
/* Equation (6) & (8) in Ref. 1 */
ent += data[ih] * ( -mu_x * log ( mu_x ) - ( 1.0 - mu_x ) * log ( 1.0 - mu_x ) );
}
}
for ( ih = it + 1; ih < 256; ih++ )
{
/* Equation (4) in Ref. 1 */
mu_x = 1.0 / ( 1.0 + term * abs ( ih - mu_1[it] ) );
if ( !((mu_x < 1e-06 ) || ( mu_x > 0.999999)))
{
/* Equation (6) & (8) in Ref. 1 */
ent += data[ih] * ( -mu_x * log ( mu_x ) - ( 1.0 - mu_x ) * log ( 1.0 - mu_x ) );
}
}
/* No need to divide by NUM_ROWS * NUM_COLS * LOG(2) ! */
if ( ent < min_ent )
{
min_ent = ent;
threshold = it;
}
}
return threshold;
}
bool bimodalTest(double *y)
{
int len=256;//y.length;
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 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
// ?at valley are unsuitable for this method.
double iHisto[256]={0};
double tHisto[256]={0};
int iter =0;
int threshold=-1;
for (int i=0; i<256; i++)
{
iHisto[i]=(double) data[i];
tHisto[i] = (double)data[i];
}
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
//iHisto = tHisto;
for(int i=0; i<256; i++)
{
iHisto[i] = tHisto[i];
}
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<255; i++)
{
if (iHisto[i-1] < iHisto[i] && iHisto[i+1] < iHisto[i])
{
tt += i;
//IJ.log("mode:" +i);
}
}
threshold = (int) cvFloor(tt/2.0);
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+=(i*data[i]);
}
threshold =cvFloor(sum/tot);
return threshold;
}
int _tmain(int argc, _TCHAR* argv[])
{
IplImage *srcImg = cvLoadImage("D:\\luosi.jpg");
IplImage *grayImg = cvCreateImage(cvGetSize(srcImg), 8, 1);
cvCvtColor(srcImg, grayImg, CV_BGR2GRAY);
hist = cvCreateHist(1, &hist_size, CV_HIST_ARRAY, ranges, 1);
cvCalcHist(&grayImg, hist, 0, NULL);
int data[256];
for(int idx = 0; idx < 256; ++idx)
{
data[idx] = cvGet1D(hist->bins, idx).val[0];
}
//
int threshold = Huang(data);
IplImage *binaryImg = cvCreateImage(cvGetSize(grayImg), 8, 1);
cvThreshold(grayImg, binaryImg, threshold, 255, CV_THRESH_BINARY);
cvNamedWindow("Huang");
cvShowImage("Huang", binaryImg);
//
threshold = Intermodes(data);
cvThreshold(grayImg, binaryImg, threshold, 255, CV_THRESH_BINARY);
cvNamedWindow("Intermodes");
cvShowImage("Intermodes", binaryImg);
//
threshold = Mean(data);
cvThreshold(grayImg, binaryImg, threshold, 255, CV_THRESH_BINARY);
cvNamedWindow("Mean");
cvShowImage("Mean", binaryImg);
cvWaitKey(0);
return 0;
}