otsu二值化

/*
OTSU 算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。下面的代码最早由 Ryan Dibble提供，此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改，补正。

转自：http://forum.assuredigit.com/display_topic_threads.asp?ForumID=8&TopicID=3480

算法对输入的灰度图像的直方图进行分析，将直方图分成两个部分，使得两部分之间的距离最大。划分点就是求得的阈值。

parameter:   *image         --- buffer for image
              rows, cols     --- size of image
              x0, y0, dx, dy   --- region of vector used for computing threshold
              vvv             --- debug option, is 0, no debug information outputed
*/
/*======================================================================*/
/*   OTSU global thresholding routine                                 */
/*   takes a 2D unsigned char array pointer, number of rows, and     */
/*   number of cols in the array. returns the value of the threshold     */
/*======================================================================*/
int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
{

unsigned char *np;     // 图像指针
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 = y0 + 1; i < y0 + dy - 1; i++) {
  np = &image[i*cols+x0+1];
  for (j = x0 + 1; j < x0 + dx - 1; j++) {
    ihist[*np]++;
    if(*np > gmax) gmax=*np;
    if(*np < gmin) gmin=*np;
    np++; /* next pixel */
  }
}

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

for (k = 0; k <= 255; k++) {
  sum += (double) k * (double) 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");
  return (160);
}

// 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 += (double) k *ihist[k];
  m1 = csum / n1;
  m2 = (sum - csum) / n2;
  sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
  /* bbg: note: can be optimized. */
  if (sb > fmax) {
    fmax = sb;
    thresholdValue = k;
  }
}

// at this point we have our thresholding value

// debug code to display thresholding values
if ( vvv & 1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
  thresholdValue, gmin, gmax);

return(thresholdValue);
}