大津法otsu的自动阈值分割源代码

/*======================================================================*/
/* OTSU global thresholding routine */
/*======================================================================*/
void otsu (IplImage *image)
{
int w = image->width;
int h = image->height;

unsigned char *np; // 图像指针
unsigned char pixel;
int thresholdValue=1; // 阈值
int ihist[256]; // 图像直方图，256个点

int i, j, k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;

// 对直方图置零...
memset(ihist, 0, sizeof(ihist));

gmin=255; gmax=0;
// 生成直方图
for (i = 0; i < h; i++) 
{
np = (unsigned char*)(image->imageData + image->widthStep*i);
for (j = 0; j < w; j++) 
{
pixel = np[j];
ihist[ pixel]++;
if(pixel > gmax) gmax= pixel;
if(pixel < gmin) gmin= pixel;
}
}

// set up everything
sum = csum = 0.0;
n = 0;

for (k = 0; k <= 255; k++) 
{
sum += k * ihist[k]; /* x*f(x) 质量矩*/
n += ihist[k]; /* f(x) 质量 */
}

if (!n) 
{
// if n has no value, there is problems...
//fprintf (stderr, "NOT NORMAL thresholdValue = 160/n");
thresholdValue = 160;
goto L;
}

// do the otsu global thresholding method
fmax = -1.0;
n1 = 0;
for (k = 0; k < 255; k++) 
{
n1 += ihist[k];
if (!n1) { continue; }
n2 = n - n1;
if (n2 == 0) { break; }
csum += k *ihist[k];
m1 = csum / n1;
m2 = (sum - csum) / n2;
sb = n1 * n2 *(m1 - m2) * (m1 - m2);
/* bbg: note: can be optimized. */
if (sb > fmax)
{
fmax = sb;
thresholdValue = k;
}
}

L:
for (i = 0; i < h; i++) 
{
np = (unsigned char*)(image->imageData + image->widthStep*i);
for (j = 0; j < w; j++) 
{
if(np[j] >= thresholdValue)
np[j] = 255;
else np[j] = 0;
}
}
}