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;

}

}

}