一種局部二值化算法：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