積分画像(Integral Image)による、Sauvolaの手法の高速化

先日の記事にて、適応的二値化の一手法であるSauvola[1]の手法について紹介しました。

閾値の計算式は単純で、実装は割と簡単です。しかし1ピクセルごとに周囲の平均・標準偏差を求めるのは、結構コストがかかる処理です。これを積分画像(Integral Image)により高速に実装する方法を紹介します。

なお、今回はC++でコードを示します。ある程度わかりやすさ重視ですので、もちろん最速実装というわけではありません。

[1] J. Sauvola et. al., “Adaptive document image binarization,” Pattern Recognition 33(2), pp.225–236, 2000.
http://www.mediateam.oulu.fi/publications/pdf/24.pdf

Sauvolaの手法 おさらい

以下の式で、各ピクセルの閾値を求めています。
T(x,y)=m(x,y) \times \left[ 1+k \times \left( \frac{s(x,y)}{R}-1 \right) \right]
あるピクセル (x, y) の閾値 T(x, y) は、その周囲の平均値 m(x, y) と、標準偏差 s(x, y) によって決まります。kとRは適当に決める定数です。

(毎回ダシにして心苦しいですが)大津の手法による二値化とSauvolaの二値化を比べてみました。元の画像、大津、Sauvolaの順です。
f:id:Schima:20131025185346p:plain
f:id:Schima:20131025185350p:plain
f:id:Schima:20131025185349p:plain


素直な実装

式と見比べながら普通に書いてみます。

T(x,y)=m(x,y) \times \left[ 1+k \times \left( \frac{s(x,y)}{R}-1 \right) \right]

void sauvolaSimple(const cv::Mat &src, cv::Mat &dst, 
   int kernelSize, double k, double r)
{
    dst.create(src.size(), src.type());

    cv::Mat srcWithBorder;
    int borderSize = kernelSize / 2 + 1;
    // 画像端の処理を簡単にするため外枠を付ける
    cv::copyMakeBorder(src, srcWithBorder, 
        borderSize, borderSize, borderSize, borderSize, cv::BORDER_REPLICATE);

    for (int y = 0; y < src.rows; y++) 
    {
        for (int x = 0; x < src.cols; x++) 
        {			
            // (x,y)周囲の部分行列を取得
            cv::Mat kernel = srcWithBorder(cv::Rect(x, y, kernelSize, kernelSize));
            // 部分行列の平均・標準偏差を求め、式から閾値を求める
            cv::Scalar mean, stddev;            
            cv::meanStdDev(kernel, mean, stddev);
            double threshold = mean[0] * (1 + k * (stddev[0] / r - 1));
            // dstに白黒を設定
            if (src.at<uchar>(y, x) < threshold)
                dst.at<uchar>(y, x) = 0;
            else
                dst.at<uchar>(y, x) = 255;
        }
    }
}
cv::Mat src = cv::imread("joho.bmp", CV_LOAD_IMAGE_GRAYSCALE);
cv::Mat dst;
sauvolaSimple(src, dst, 11, 0.15, 50);

処理時間は、私の環境では 800x600 程度の画像でも1秒以上かかってしまいました。やはり、1ピクセルごとの平均・標準偏差計算が足を引っ張っています。

積分画像による実装

方法の概要

平均・標準偏差(分散)の計算には、画素値の和・二乗和が必要です。これの計算を効率よくできるのが積分画像(Integral Image)です。
3x4の画像での積分画像の例を、下に作ってみました。積分画像のそれぞれのピクセルの値は、左上の隅(0,0)からそのピクセルの位置までの画素値を総和したものになっています。例えば、積分画像の(行:2 列:2)にある値は 14 になっていて、これは元の画像の 1 + 2 + 5 + 6 です。

元の画像
1 2 3 4
5 6 7 8
9 10 11 12
積分画像
1 3 6 10
6 14 24 36
15 33 54 78

この積分画像で、「あるピクセルの周りの画素値の総和」が求められます。

下のがんばって作ったもののイマイチな図で、今求めたいのが右下の赤の部分の画素の総和とします。(x,y) からブロックサイズ分の周囲を広げているかたちです。
ここの値は直接積分画像には入っていませんが、A,B,C,Dの部分は求められています。これをつかうと、求めたい赤い部分Rは、以下の式で求められます。
 R=D-B-C+A
全体を含むDからBとCを引き、引きすぎた分Aを足しています。

f:id:Schima:20131025200041p:plain

実装

上の図的には、(x-w, y-w)から(x+w, y+w)の範囲の画素値総和を見ていくことになりますが、copyMakeBorderによりちょうどwの太さのボーダーを作っているため、実際にアクセスする部分がずれて記述は以下のようにシンプルになります。

エラー処理はありませんが、kernelSizeは奇数限定です。

void sauvolaFast(const cv::Mat &src, cv::Mat &dst, int kernelSize, double k, double r)
{
    dst.create(src.size(), src.type());

    cv::Mat srcWithBorder;
    int borderSize = kernelSize / 2 + 1;
    int kernelPixels = kernelSize * kernelSize;
    cv::copyMakeBorder(src, srcWithBorder, borderSize, borderSize, 
                       borderSize, borderSize, cv::BORDER_REPLICATE);

    cv::Mat sum, sqSum;
    cv::integral(srcWithBorder, sum, sqSum);	
    for(int y = 0; y < src.rows; y++)
    {
        for(int x = 0; x < src.cols; x++)
        {
            int kx = x + kernelSize;
            int ky = y + kernelSize;
            double sumVal = sum.at<int>(ky, kx)
                          - sum.at<int>(ky, x)
                          - sum.at<int>(y, kx)
                          + sum.at<int>(y, x);
            double sqSumVal = sqSum.at<double>(ky, kx)
                            - sqSum.at<double>(ky, x)
                            - sqSum.at<double>(y, kx)
                            + sqSum.at<double>(y, x);

            double mean = sumVal / kernelPixels;
            double var = (sqSumVal / kernelPixels) - (mean * mean);
            if (var < 0.0) 
                var = 0.0;
            double stddev = sqrt(var);
            double threshold = mean * (1 + k * (stddev / r - 1));

            if (src.at<uchar>(y, x) < threshold)
                dst.at<uchar>(y, x) = 0;
            else
                dst.at<uchar>(y, x) = 255;
        }
    }
}

実行速度

Core-i7-3520TM 2.90GHzのCPUで、OpenCV 2.4.5を使用、対象は上に例で示した画像(816x612)、Releaseビルドで実行した結果が以下です。

普通の実装
1.118秒
積分画像による実装
0.052秒

およそ20~30倍ほど高速になりました。効果てきめん。以下計測に使ったプログラムです。

int main(int argc, _TCHAR* argv[])
{
    const char *fileName = "C:\\joho.bmp";

    cv::Mat src = cv::imread(fileName, CV_LOAD_IMAGE_GRAYSCALE);
    cv::Mat dst1, dst2;

    cv::TickMeter meter;
    {
        meter.start();
        sauvolaSimple(src, dst1, 11, 0.15, 50);
        meter.stop();
    }
    std::cout << "time(simple): " << meter << std::endl;

    meter.reset();
    {
        meter.start();
        sauvolaFast(src, dst2, 11, 0.15, 50);
        meter.stop();
    }
    std::cout << "time(integral): " << meter << std::endl;

    // 同じ結果かどうか
    std::cout << "different pixels: " << cv::countNonZero(dst1 - dst2) << std::endl;

    cv::imshow("src", src);
    cv::imshow("sauvola simple", dst1);
    cv::imshow("sauvola integral", dst2);
    cv::waitKey();
    cv::destroyAllWindows();

    return 0;
}

OpenCvSharpでの実装

先日の記事にて説明した通り、Sauvolaの手法による二値化はOpenCvSharp.Extensions.Binarizerクラスから使うことができます。

最初に述べた普通の実装をしたものがBinarizer.Sauvolaメソッド、積分画像により実装したのがBinarizer.SauvolaFastメソッドです。普通はSauvolaFastを使えば良いはずです。