sauvola二值化算法研究

sauvola二值化算法研究

 

sauvola是一种考虑局部均值亮度的图像二值化方法, 以局部均值为基准在根据标准差做些微调.算法实现上一般用积分图方法

来实现.这个方法能很好的解决全局阈值方法的短板—关照不均图像二值化不好的问题.先贴代码

//************************************ // 函数名称: sauvola // 函数说明: 局部均值二值化 // 参    数: //           const unsigned char * grayImage        [in]        输入图像数据 //           const unsigned char * biImage          [out]       输出图像数据      //           const int w                            [in]        输入输出图像数据宽 //           const int h                            [in]        输入输出图像数据高 //           const int k                            [in]        threshold = mean*(1 + k*((std / 128) - 1)) //           const int windowSize                   [in]        处理区域宽高 // 返 回 值: void //************************************ void sauvola(const unsigned char * grayImage, const unsigned char * biImage,     const int w, const int h, const int k, const int windowSize){     int whalf = windowSize >> 1;     int i, j;     int IMAGE_WIDTH = w;     int IMAGE_HEIGHT = h;     // create the integral image     unsigned long * integralImg = (unsigned long*)malloc(IMAGE_WIDTH*IMAGE_HEIGHT*sizeof(unsigned long*));     unsigned long * integralImgSqrt = (unsigned long*)malloc(IMAGE_WIDTH*IMAGE_HEIGHT*sizeof(unsigned long*));     int sum = 0;     int sqrtsum = 0;     int index;     //收集数据 integralImg像素和积分图 integralImgSqrt像素平方和积分图     for (i = 0; i < IMAGE_HEIGHT; i++){         // reset this column sum         sum = 0;         sqrtsum = 0;         for (j = 0; j < IMAGE_WIDTH; j++)         {             index = i*IMAGE_WIDTH + j;             sum += grayImage[index];             sqrtsum += grayImage[index] * grayImage[index];             if (i == 0){                 integralImg[index] = sum;                 integralImgSqrt[index] = sqrtsum;             }             else{                 integralImgSqrt[index] = integralImgSqrt[(i - 1)*IMAGE_WIDTH + j] + sqrtsum;                 integralImg[index] = integralImg[(i - 1)*IMAGE_WIDTH + j] + sum;             }         }     }     //Calculate the mean and standard deviation using the integral image     int xmin, ymin, xmax, ymax;     double mean, std, threshold;     double diagsum, idiagsum, diff, sqdiagsum, sqidiagsum, sqdiff, area;     for (i = 0; i < IMAGE_WIDTH; i++){         for (j = 0; j < IMAGE_HEIGHT; j++){             xmin = max(0, i - whalf);             ymin = max(0, j - whalf);             xmax = min(IMAGE_WIDTH - 1, i + whalf);             ymax = min(IMAGE_HEIGHT - 1, j + whalf);             area = (xmax - xmin + 1) * (ymax - ymin + 1);             if (area <= 0){                 biImage[i * IMAGE_WIDTH + j] = 255;                 continue;             }             if (xmin == 0 && ymin == 0){                 diff = integralImg[ymax * IMAGE_WIDTH + xmax];                 sqdiff = integralImgSqrt[ymax * IMAGE_WIDTH + xmax];             }             else if (xmin > 0 && ymin == 0){                 diff = integralImg[ymax * IMAGE_WIDTH + xmax] - integralImg[ymax * IMAGE_WIDTH + xmin - 1];                 sqdiff = integralImgSqrt[ymax * IMAGE_WIDTH + xmax] - integralImgSqrt[ymax * IMAGE_WIDTH + xmin - 1];             }             else if (xmin == 0 && ymin > 0){                 diff = integralImg[ymax * IMAGE_WIDTH + xmax] - integralImg[(ymin - 1) * IMAGE_WIDTH + xmax];                 sqdiff = integralImgSqrt[ymax * IMAGE_WIDTH + xmax] - integralImgSqrt[(ymin - 1) * IMAGE_WIDTH + xmax];;             }             else{                 diagsum = integralImg[ymax * IMAGE_WIDTH + xmax] + integralImg[(ymin - 1) * IMAGE_WIDTH + xmin - 1];                 idiagsum = integralImg[(ymin - 1) * IMAGE_WIDTH + xmax] + integralImg[ymax * IMAGE_WIDTH + xmin - 1];                 diff = diagsum - idiagsum;                 sqdiagsum = integralImgSqrt[ymax * IMAGE_WIDTH + xmax] + integralImgSqrt[(ymin - 1) * IMAGE_WIDTH + xmin - 1];                 sqidiagsum = integralImgSqrt[(ymin - 1) * IMAGE_WIDTH + xmax] + integralImgSqrt[ymax * IMAGE_WIDTH + xmin - 1];                 sqdiff = sqdiagsum - sqidiagsum;             }             mean = diff / area;             std = sqrt((sqdiff - diff*diff / area) / (area - 1));             threshold = mean*(1 + k*((std / 128) - 1));             if (grayImage[j*IMAGE_WIDTH + i] < threshold)                 biImage[j*IMAGE_WIDTH + i] = 0;             else                 biImage[j*IMAGE_WIDTH + i] = 255;         }     }     free(integralImg);     free(integralImgSqrt); }

 

代码要注意下面几点:

1 计算区域像素和,几乎使用积分图技术是必然的选择.

2 标准差的表示方法: std = sqrt((sqdiff - diff*diff / area) / (area - 1)) 终于感到高等代数没有白学,

可以看百度百科关于方差的说明

http://baike.baidu.com/link?url=uFltaqvwLYZHvCO4-IJipF89x8-EhuEfZW12lAnES_TGWyhG62ntmWKTcVs511PSfzE7nanQzgl37rKFMOwwYq

 

3 判定方程 threshold = mean*(1 + k*((std / 128) - 1)). 首先均值是基础, 如果标准差大写,阈值就会大些,标准差小些,阈值就会小些.

这个方法对一些不是光照不均的图片有时候效果不好,现在还在找较好的方法,初步打算先用全局均值做二值化,如何效果不好再用局部均值的方法.

----为什么我的博客没人看啊……………………………………

 

转载于:https://www.cnblogs.com/guopengfei/p/4766526.html