DSP移植-OTSU算法

有关OTSU算法的原理，网络上有很过资源，描述得也很清楚，这里引用已有代码，稍加修改后（原来的代码可能有溢出，无法找到正确的阈值），能够正确运行。

两段代码的效果类似，第一个函数的阈值是149， 第二个函数的阈值150，第一个函数的运算量要相对较小，优化后可以移植到到DSP。

源代码：

//http://blog.csdn.net/qq_26898461/article/details/46930183
//http://www.cnblogs.com/skyseraph/archive/2010/12/21/1913058.html
/*======================================================================*/
/* OTSU global thresholding routine */
/*======================================================================*/
int victOtsuThreshold(IplImage *image, int *threshold)
{
	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;
	double u1, u2;

	// 对直方图置零...
	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 = 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...
		*threshold =160;
		return 0;
	}
	// do the otsu global thresholding method
	fmax =-1.0;
	n1 = 0;
	n2 = 0;
	csum = 0.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;
		u1 = n1 / (double)n;
		u2 = n2 / (double)n;
		sb =  u1 * u2 *(m1 - m2) * (m1 - m2);
		/* bbg: note: can be optimized. */
		if (sb > fmax)
		{
			fmax = sb;
			thresholdValue = k;
		}
	}
	*threshold = thresholdValue;
	return 0;
}



int victOtsuThreshold_2(IplImage* src, int *threshold)
{    
	int height=src->height;    
	int width=src->width;        
	long size = height * width;   

	//histogram    
	float histogram[256] = {0};    
	for(int m=0; m < height; m++)  
	{    
		unsigned char* p=(unsigned char*)src->imageData + src->widthStep * m;    
		for(int n = 0; n < width; n++)   
		{    
			histogram[int(*p++)]++;    
		}    
	}    

	int threshold_temp;      
	long sum0 = 0, sum1 = 0; //存储前景的灰度总和和背景灰度总和  
	long cnt0 = 0, cnt1 = 0; //前景的总个数和背景的总个数  
	double w0 = 0, w1 = 0; //前景和背景所占整幅图像的比例  
	double u0 = 0, u1 = 0;  //前景和背景的平均灰度  
	double variance = 0; //最大类间方差  
	int i, j;  
	double u = 0;  
	double maxVariance = 0;  
	for(i = 1; i < 256; i++) //一次遍历每个像素  
	{    
		sum0 = 0;  
		sum1 = 0;   
		cnt0 = 0;  
		cnt1 = 0;  
		w0 = 0;  
		w1 = 0;  
		for(j = 0; j < i; j++)  
		{  
			cnt0 += histogram[j];  
			sum0 += j * histogram[j];  
		}  

		u0 = (double)sum0 /  cnt0;   
		w0 = (double)cnt0 / size;  

		for(j = i ; j <= 255; j++)  
		{  
			cnt1 += histogram[j];  
			sum1 += j * histogram[j];  
		}  

		u1 = (double)sum1 / cnt1;  
		w1 = 1 - w0; // (double)cnt1 / size;  

		u = u0 * w0 + u1 * w1; //图像的平均灰度  
		printf("u = %f\n", u);  
		//variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);  
		variance =  w0 * w1 *  (u0 - u1) * (u0 - u1);  
		if(variance > maxVariance)   
		{    
			maxVariance = variance;    
			threshold_temp = i;    
		}   
	}    

	printf("threshold = %d\n", threshold_temp); 
	*threshold = threshold_temp;
	return 0; 
}