インテグラルイメージによるボックスフィルタの高速化

04 10月 2010 Under: opencv2.x-samples

フィルタリングは一般的に,カーネルサイズが大きくなればなるほど計算時間がかかるようになります.
ここでは,インテグラルイメージを使って,カーネルサイズに依存しない効率的なボックスフィルタを実装します.
ただし,OpenCVで実装されているボックスフィルタであるboxFilter関数のほうがより高速に動作します.
インテグラルイメージとは,原点(0,0)からピクセルp=(x,y)までの総和した画像
(integralimage(x,y)=Σ0xΣ0y im(p,q))
のことです.OpenCV(2.x)ではintegral関数で呼び出すことが可能です.

C++

#include <iostream>
#include <fstream>
#include <cv.h>
#include <highgui.h>
using namespace std;
using namespace cv;

void myboxBlur(Mat& src, Mat& dest, Size ksize);
void myintegral(Mat& src, Mat& dest);
void boxBlurIntegral(Mat& src, Mat& dest, Size ksize);

void speedTest(Mat& src);

int main(int argc, char** argv)
{
	Mat src = imread("garden.png");
	if(src.empty())cout<<"please input valid file"<<endl;

	Mat destMY;//my blur function
	Mat destII;//integral image
	Mat destCV;//OpenCV blur function
	Size kernel = Size(15,15);

	//(1) my basic implimentaion of box blur filtering
	int64 t = getTickCount();
	myboxBlur(src,destMY,kernel);
	cout<<"my box"<<(getTickCount()-t)*1000.0/getTickFrequency()<<"ms"<<endl;
	imshow("test",destMY);
	waitKey();

	//(3)integral image based implimentaion
	t = getTickCount();
	boxBlurIntegral(src,destII,kernel);
	cout<<"integral box"<<(getTickCount()-t)*1000.0/getTickFrequency()<<"ms"<<endl;
	imshow("test",destII);
	waitKey();

	t = getTickCount();
	//(5) OpenCV original implimentation
	boxFilter(src,destCV,-1,kernel);
	cout<<"cv box"<<(getTickCount()-t)*1000.0/getTickFrequency()<<"ms"<<endl;
	imshow("test",destCV);
	waitKey();

	//(6) test function of computational cost  
	speedTest(src);

	return 0;
}

void myboxBlur(Mat& src, Mat& dest, Size ksize)
{
	if(dest.empty())dest.create(src.size(),src.type());

	const int kh = ksize.height/2;
	const int kw = ksize.width/2;
	const double div  = 1.0/ksize.area();
	Mat im;
	//(2) make border image for outof range of kernel window
	copyMakeBorder(src,im,kh,kh,kw,kw,BORDER_REPLICATE);

#pragma omp parallel for
	for(int j=0;j<dest.rows;j++)
	{
		uchar* d = dest.ptr<uchar>(j);
		for(int i=0;i<dest.cols;i++)
		{
			int r=0;
			int g=0;
			int b=0;
			for(int l=0;l<ksize.height;l++)
			{
				uchar* s = im.ptr<uchar>(j+l);
				s+=3*i;
				for(int k=0;k<ksize.width;k++)
				{
					r += s[0];
					g += s[1];
					b += s[2];
					s+=3;
				}
			}
			d[0]=saturate_cast<uchar>(r*div);
			d[1]=saturate_cast<uchar>(g*div);
			d[2]=saturate_cast<uchar>(b*div);
			d+=3;
		}
	}
}

void boxBlurIntegral(Mat& src, Mat& dest, Size ksize)
{
	if(dest.empty())dest.create(src.size(),src.type());	
	const int kh = ksize.height/2;
	const int kw = ksize.width/2;
	const double div  = 1.0/ksize.area();
	Mat im;
	copyMakeBorder(src,im,kh,kh,kw,kw,BORDER_REPLICATE);

	Mat sum;
	integral(im,sum,CV_32SC3);
	const int step = 3*ksize.width;
#pragma omp parallel for
	for(int j=0;j<dest.rows;j++)
	{
		uchar* d = dest.ptr<uchar>(j);
		int* sp=sum.ptr<int>(j+ksize.height);
		int*  s=sum.ptr<int>(j);
		for(int i=0;i<dest.cols;i++)
		{
			//(4) get box blur value from integral image
			d[0]=saturate_cast<uchar>((s[0]-s[step+0]-sp[0]+sp[step+0])*div);
			d[1]=saturate_cast<uchar>((s[1]-s[step+1]-sp[1]+sp[step+1])*div);
			d[2]=saturate_cast<uchar>((s[2]-s[step+2]-sp[2]+sp[step+2])*div);
			d+=3;
			s+=3;
			sp+=3;
		}
	}
}

void speedTest(Mat& src)
{
	ofstream plot("plot");
	const int maxiter=10;

	Mat dest;

	for(int k=0;k<20;k++)
	{
		cout<<2*k+1<<endl;
		Size kernel = Size(2*k+1,2*k+1);

		double t1=0.0;
		double t2=0.0;
		double t3=0.0;
		for(int i=0;i<maxiter;i++)
		{
			int64 t = getTickCount();
			myboxBlur(src,dest,kernel);
			t1 += (getTickCount()-t)*1000.0/getTickFrequency();

			t = getTickCount();
			boxBlurIntegral(src,dest,kernel);
			t2 += (getTickCount()-t)*1000.0/getTickFrequency();

			t = getTickCount();
			boxFilter(src,dest,-1,kernel);
			t3 += (getTickCount()-t)*1000.0/getTickFrequency();

		}
		cout<<"my       box"<<t1/maxiter<<"ms"<<endl;
		cout<<"integral box"<<t2/maxiter<<"ms"<<endl;
		cout<<"cv       box"<<t3/maxiter<<"ms"<<endl;

		plot<<2*k+1<<" "<<t1/maxiter<<" "<<t2/maxiter<<" "<<t3/maxiter<<endl;
	}
}

(1)ボックスフィルタの古典的な実装

カーネルサイズに合わせて4重ループでフィルタを構成する,最も一般的な方法の実装です.
カーネルのウインドウが大きくなるほど計算コストは大きくなります.

(2)copyMakeBorder関数による縁のコピー

カーネル等を使うと画像の周辺領域で例外処理をする必要性が出てきます.この関数はその例外処理を避けるために,はみ出す予定の分だけ大きな画像を作り直します.
カーネルサイズが11×11だった場合,縦横に,10ピクセルづつ大きな画像が確保されます.

(3)インテグラルイメージによるボックスフィルタの実装

一定時間で動作するインテグラルイメージによるボックスフィルタの実装です.
まずはじめに,積分画像(インテグラルイメージ)を生成し,そのあとカーネルサイズに合わせて必要な値を取っていきます.
integral関数の詳細はreferenceを参照してください.

(4)インテグラルイメージからの値取り出し

カーネル中心から,カーネルウインドウサイズに応じた4つの角の値を加減算することで区間積分した値を取り出します.
copyMakeBorder関数によりオフセットがついたため,出力のカーネル中心と座標からカーネル半径だけ縦横にシフトした値からなる四角形の4つ角を使って計算します.

(5)OpenCVのオリジナルのボックスフィルタ

OpenCVのボックスフィルタの実装です.

(6)結果のグラフの計算コストプロット関数

実行結果内の計算速度のグラフを計算する関数です.

実行結果

入力画像,出力画像,またカーネルサイズに対する各手法の実行時間を以下に示します.(使用した計算機は,Core i7 920 3.2GHz)
また,OpenMPにより並列化した結果も同時に示します.
並列化なしでは,opencvの実装が10.8ms インテグラルイメージの実装が12.5msほどです.OpenMPによる並列化をするとインテグラルイメージの実装が7.8msほどになります.古典的な手法は,並列化により高速化されていますが,一定時間アルゴリズムよりも速くは動きません.

古典的な手法が,カーネルサイズが大きくなるほど計算時間がかかるようになっているのに対して,インテグラルイメージによる手法や,OpenCVの実装ではサイズに依存しない方法になっていることが分かります.

入力画像

出力画像

各手法の計算結果(並列化なし)


各手法の計算結果(並列化あり)