/*======================================================================*//* 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 countersint 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 everythingsum = 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 methodfmax = -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;}}}
