DSP

DSP移植-OTSU算法

2019-07-13 10:47发布

站内文章 / DSP
16169 0 1161
有关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 ", 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 ", threshold_temp); *threshold = threshold_temp; return 0; }

热门文章