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;
}