積分画像(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の手法 おさらい
以下の式で、各ピクセルの閾値を求めています。
あるピクセル (x, y) の閾値 T(x, y) は、その周囲の平均値 m(x, y) と、標準偏差 s(x, y) によって決まります。kとRは適当に決める定数です。
（毎回ダシにして心苦しいですが）大津の手法による二値化とSauvolaの二値化を比べてみました。元の画像、大津、Sauvolaの順です。
素直な実装
式と見比べながら普通に書いてみます。
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 |
実装
上の図的には、(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; }
