一种局部二值化算法：Sauvola算法

之前接触过全局二值化（OTSU算法），还有OPENCV提供的自适应二值化，最近又了解到一种新的局部二值化算法，Sauvola算法。

转载自：http://www.dididongdong.com/archives/4048

值得注意的是，计算r×r邻域内像素值的时候，一种优化的策略是，使用OPENCV提供的积分图，计算整张图像的积分图，那么计算r×r区域内的均值可以在常数时间内实现。

CV_EXPORTS_W void integral( InputArray src, OutputArray sum, int sdepth = -1 );

我们常见的图像二值化算法大致可分为全局阈值方法与局部阈值方法这两种类型。其中OTSU算法是全局阈值的代表，而Sauvola算法则是局部阈值方法的标杆。Sauvola算法的输入是灰度图像，它以当前像素点为中心，根据当前像素点邻域内的灰度均值与标准方差来动态计算该像素点的阈值。

假定当前像素点的坐标为(x,y)，以该点为中心的领域为r*r，g(x,y)表示(x,y)处的灰度值，Sauvola算法的步骤为：

来自参考《基于光照不均匀图像的自适应二值化方法研究》郭佳

引用

在二值化的操作中，用的比较多的就是全局阈值话OTSU（大津法）和局部阈值NiBlack，Niblack方法是一种简单有效的动态阈值分割方法，修改得到最佳参数之后的效果比大津法要好，因为大津法是根据整个图像来确定一个阈值，而Niblack则是在不同的R*R区域会有不同的阈值。

Niblack的基本思想是：对于图像的每一个像素点，在rxr领域空问里，计算该像素点领域方位内其他像素点的均值和方差。然后利用公式（1）进行二值化。

其中，T（x,y）是阈值，k是预先设定的修正值,图像为f(x,y)，均值为m(x,y)，方差为s(x,y)。

使用Niblack法的优点在于：

对每一个像素点都独立的跟据此像素点的邻域的情况来计算门限，对于和邻域均值m(x，y)相近的像素点判断为背景而反之判断为前景；而具体相近到什么程度由标准差s(X’y)和修正系数k来决定，这保证了这种的方法的灵活性。

使用Niblack法的不足在于：

由于要利用域r×r模板遍历图像，导致边界区域（r-1）/2的像素范围内无法求取阈值；同时当进行图像遍历时，如果域r×r范围内都是背景，经NIBLACK计算后必有一部分被确定为目标，产生伪噪声。

总之，用Niblack方法进行图像分割时，选择的处理模板窗口R*R大小的选择很关键，选择的空间太小，则噪声抑制的效果不理想，目标主体不够突出，选择的空间太大，则目标的细节会被去除而丢失信息。

参考

1.第一次做MathOCR遇到的参考文献：
《图片中印刷体数学公式的自动识别》陈颂光
2.中文的的链接来自，追溯不到原文：

https://livezingy.com/derivations-of-sauvola-formula/
3.英文文献：值得参考的标准
Efficient implementation of local adaptive thresholding techniques using integral images
PDF链接为：https://pdfs.semanticscholar.org/8130/a9499715d22468492c3786c34ba1ba0b4ed3.pdf
4.Matlab代码参考
http://freesourcecode.net/matlabprojects/59687/sauvola-local-image-thresholding-in-matlab#.Wzsk2oq-vcs

5,https://www.cnblogs.com/guopengfei/p/4766526.html

//求区域内均值 integral即为积分图
float fastMean(cv::Mat& integral, int x, int y, int window)
{

    int min_y = std::max(0, y - window / 2);
    int max_y = std::min(integral.rows - 1, y + window / 2);
    int min_x = std::max(0, x - window / 2);
    int max_x = std::min(integral.cols - 1, x + window / 2);

    int topright = integral.at<int>(max_y, max_x);
    int botleft = integral.at<int>(min_y, min_x);
    int topleft = integral.at<int>(max_y, min_x);
    int botright = integral.at<int>(min_y, max_x);

    float res = (float)((topright + botleft - topleft - botright) / (float)((max_y - min_y) *(max_x - min_x)));

    return res;
}


cv::Mat& Sauvola(cv::Mat& inpImg, cv::Mat& resImg,  int window, float k)
{
    cv::Mat integral;
    int nYOffSet = 3;
    int nXOffSet = 3;
    cv::integral(inpImg, integral);  //计算积分图像
    for (int y = 0; y < inpImg.rows; y += nYOffSet)
    {
        for (int x = 0; x < inpImg.cols; x += nXOffSet)
        {

            float fmean = fastMean(integral, x, y, window);float fthreshold = (float)(fmean*(1.0 - k));  
            
            int nNextY = y + nYOffSet;
            int nNextX = x + nXOffSet;
            int nCurY = y;
            while (nCurY < nNextY && nCurY < inpImg.rows)
            {
                int nCurX = x;
                while (nCurX < nNextX && nCurX < inpImg.cols)
                {
                    uchar val = inpImg.at<uchar>(nCurY, nCurX) < fthreshold;
                    resImg.at<uchar>(nCurY, nCurX) = (val == 0 ? 0 : 255);
                    nCurX++;
                }
                nCurY++;
            }

        }
    }

    return resImg;
}

//************************************

// 函数名称: 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);

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

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

代码要注意下面几点:

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

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

 

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

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

以上转载自：

https://www.cnblogs.com/guopengfei/p/4766526.html