动态阙值

int ihist[256] = {0}; // 图像直方图，256个点
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 i, j, k; // various counters
	int n, n1, n2, gmin, gmax;
	double m1, m2, sum, csum, fmax, sb;
		
	//初始化灰度直方图
//	for (k = 0; k <= 255; k++) 
//	{
//		ihist[k] = 0;
//	}
	//gmin=255; gmax=0;//图像阈值范围设置初始
	gmin=200; gmax=60;
	// 生成直方图
	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;
	   }
	}

	return(thresholdValue);
}