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)). 首先均值是基础, 如果标准差大写,阈值就会大些,标准差小些,阈值就会小些.

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

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