スパース行列(SparseMat)の使い方:画像のぼかしとぼけ除去

02 7月 2010 Under: opencv2.x-samples

画像のぼけの過程は巨大なスパース行列を使って表現できます.
入力画像x(サイズは,S=width*height),ぼけ画像yを1次元ベクトルX,Y (画像を縦横を区別せずならべたもの)としたものをぼけ過程を表す行列Hで表すと,

Y=HX

で表せます.この行列Hは画像サイズ×画像サイズ(SxS)の巨大な行列になり,Matクラスなどで確保すると あっという間にメモリが枯渇します.幸いなことに,この行列Hは要素がほとんどゼロなため,要素が非ゼロの 場所だけメモリを使用するスパース行列クラスSparseMatが使用可能です.ここでは, SparseMatの使い方を画像のぼかしとその逆過程であるぼけ除去を例に説明します.

C++


#include <cv.h>
#include <highgui.h>
using namespace cv;

#include <iostream>
using namespace std;

SparseMat createblurBoxSparseMat32f(Mat& src, Size window);
Mat visualizeSparseMat(SparseMat& smat, Size out_imsize);

void blurSparseMat32f(Mat& src, Mat& dest, SparseMat& H);
void deblurSparseMat32f(Mat& src, Mat& dest, SparseMat& H, int iteration, float beta);
void mulSparseMat32f(SparseMat& smat, Mat& src, Mat& dest, bool isTranspose = false);

void unsharpMaskBox(Mat& src, Mat& dest,Size window);

int main(int argc, char** argv)
{
	Mat image = imread("lenna.png");//input image
	Mat bimage(image.size(),CV_8UC3);// blured image
	Mat dbimage(image.size(),CV_8UC3);// deblured blured image
	Mat dbimage2(image.size(),CV_8UC3);// deblured blured image for unsharp mask

	Size kernelwin = Size(7,7);//blur kernel is 7x7 box matrix
	int64 time;
	double f = 1000.0/getTickFrequency();

	//(1) generate sparse matrix of point spread function
	time = getTickCount();
	SparseMat psf=createblurBoxSparseMat32f(image, kernelwin);
	cout<<"generating point spread fuction time:"<<(getTickCount()-time)*f<<" ms"<<endl;

	//(2) visualize the sparse matrix
	time = getTickCount();
	Mat smatim = visualizeSparseMat(psf, Size(256,256));
	cout<<"visualizing sparse matrix time:"<<(getTickCount()-time)*f<<" ms"<<endl;
	imshow("Sparse Mat",smatim);
	waitKey(10);

	//(3) blur a input image  with the sparse matrix
	time = getTickCount();
	blurSparseMat32f(image,bimage,psf);
	cout<<"bluring time:"<<(getTickCount()-time)*f<<" ms"<<endl;

	//(4) deblur the degraded image with a steepest descent method for L2 norm minimization
	time = getTickCount();
	deblurSparseMat32f(bimage,dbimage,psf,10, 0.2f);//iteration:10, beta = 0.2
	cout<<"debluring time:"<<(getTickCount()-time)*f<<" ms"<<endl;

	//(5) sharp the degraded image with common sharpness filter for mutch up
	time = getTickCount();
	unsharpMaskBox(bimage,dbimage2,kernelwin);
	cout<<"unsharpMask time:"<<(getTickCount()-time)*f<<" ms"<<endl;

	imshow("test",image); //input image
	waitKey();
	imshow("test",bimage);//blured image
	waitKey();
	imshow("test",dbimage);//deblured image
	waitKey();
	imshow("test",dbimage2);//unsharpmasked image
	waitKey();

	//(6)evaluation by L2 norm
	cout<<"norm: ideal <-> blured                "<<norm(image,bimage)<<endl;
	cout<<"norm: ideal <-> deblured              "<<norm(image,dbimage)<<endl;
	cout<<"norm: ideal <-> unsharpMask "<<norm(image,dbimage2)<<endl;
	return 0;
}

Mat visualizeSparseMat(SparseMat& smat, Size out_imsize)
{
	Mat data = Mat::zeros(out_imsize,CV_8U);
	double inv_size = 1.0/smat.size(0); 

	SparseMatIterator it = smat.begin(),it_end = smat.end();
	for(;it!=it_end;++it)
	{
		int j = (int)(((double)it.node()->idx[0]*inv_size*out_imsize.height));
		int i = (int)(((double)it.node()->idx[1]*inv_size*out_imsize.width));
		data.at<uchar>(j,i)++;
	}

	double minv,maxv;
	minMaxLoc(data,&minv,&maxv);
	data.convertTo(data,-1,255.0/maxv);

	Mat zeromat = Mat::zeros(out_imsize,CV_8U);
	vector<Mat> image;
	image.push_back(zeromat);
	image.push_back(data);
	image.push_back(zeromat);

	Mat ret;
	cv::merge(image,ret);

	cout<<"number of non zero elements: "<<smat.nzcount()<<endl;	
	return ret;
}

SparseMat createblurBoxSparseMat32f(Mat& src, Size window)
{
	int matsize =src.cols*src.rows ;
	int size[2]={matsize,matsize};
	SparseMat H(2,size,CV_32FC1);//sparse mat: image_size x image_size

	int bsy=window.height/2;// offset for kernel window
	int bsx=window.width/2;// offset for kernel window

	Mat kernel(window, CV_32FC1);
	kernel.setTo(1);
	normalize(kernel,kernel,1,0,NORM_L1);
	for(int j=0;j<src.rows;j++)
	{
		for(int i=0;i<src.cols;i++)
		{
			int y = src.cols*j+i;
			if(i>=bsx &&i<src.cols-bsx &&j>=bsy &&j<src.rows-bsy)
			{
				for(int l=0;l<window.height;l++)
				{
					for(int k=0;k<window.height;k++)
					{
						// higher-level element access functions:
						// ref<_Tp>(i0,...[,hashval]) - equivalent to *(_Tp*)ptr(i0,...true[,hashval]).
						//    always return valid reference to the element.
						//    If it's did not exist, it is created.
						// find<_Tp>(i0,...[,hashval]) - equivalent to (_const Tp*)ptr(i0,...false[,hashval]).
						//    return pointer to the element or NULL pointer if the element is not there.
						// value<_Tp>(i0,...[,hashval]) - equivalent to
						//    { const _Tp* p = find<_Tp>(i0,...[,hashval]); return p ? *p : _Tp(); }
						//    that is, 0 is returned when the element is not there.
						// note that _Tp must match the actual matrix type -
						// the functions do not do any on-fly type conversion
						H.ref<float>(y,y+src.cols*(l-bsy)+(k-bsx))=kernel.at<float>(l,k);
					}
				}
			}
			else //exception handling
			{
				H.ref<float>(y,y)=1.0;
			}
		}
	}
	return H;
}


void mulSparseMat32f(SparseMat& smat, Mat& src, Mat& dest, bool isTranspose)
{
	dest.setTo(0);
	SparseMatIterator it = smat.begin(),it_end = smat.end();
	if(!isTranspose)
	{
		for(;it!=it_end;++it)
		{	
			float* d = dest.ptr<float>(it.node()->idx[0]);
			float* s = src.ptr<float>(it.node()->idx[1]);
			for(int color=0;color<3;color++)
			{
				d[color]+= it.value<float>() * s[color];
			}
		}
	}
	else // for transpose matrix
	{
		for(;it!=it_end;++it)
		{	
			float* d = dest.ptr<float>(it.node()->idx[1]);
			float* s = src.ptr<float>(it.node()->idx[0]);
             for(int color=0;color<3;color++)
			{
				d[color]+= it.value<float>() * s[color];
			}
		}
	}
}

void deblurSparseMat32f(Mat& src, Mat& dest, SparseMat& H, int iteration, float beta)
{
	const int matsize = src.cols*src.rows;
	Mat dstvec(matsize,1,CV_32FC3);
	Mat svec; 
	//convert Mat image structure to 1D vecor structure
	src.reshape(3,matsize).convertTo(svec,CV_32FC3);

	//steepest descent method for L2 norm minimization
	for(int i=0;i<iteration;i++)
	{
		//blur input vector
		mulSparseMat32f(H,svec,dstvec);
		//subtract b blured image from non-blured image
		Mat temp = dstvec - svec;
		//blur the subtructed vector with transposed matrix
		mulSparseMat32f(H,temp,dstvec,true);
		//creep ideal image, beta is parameter of the creeping speed.
		svec -= (beta*dstvec);
	}
	//re-convert  1D vecor structure to Mat image structure
	svec.reshape(3,src.rows).convertTo(dest,CV_8UC3);
}

void blurSparseMat32f(Mat& src, Mat& dest, SparseMat& H)
{	
	const int matsize = src.cols*src.rows;
	Mat dstvec(matsize,1,CV_32FC3);
	//convert Mat image structure to 1D vecor structure
	Mat svec = src.reshape(3,matsize);

	//bluring operation
	SparseMatIterator it = H.begin(),it_end = H.end();
	for(;it!=it_end;++it)
	{	
		float* d = dstvec.ptr<float>(it.node()->idx[0]);
		uchar* s = svec.ptr<uchar>(it.node()->idx[1]);
        
        for(int color=0;color<3;color++)
		{
			d[color]+= it.value<float>() * s[color];
		}
	}

	//re-convert  1D vecor structure to Mat image structure
	dstvec.reshape(3,src.rows).convertTo(dest,CV_8UC3);
}

//unsharp mask for Mat image; it is not matrix domain representation, just image filtering.
void unsharpMaskBox(Mat& src, Mat& dest,Size window)
{
	Mat src_f;
	Mat dest_f;
	Mat sub_f;

	src.convertTo(src_f,CV_32F);

	blur(src_f,dest_f, window);
	subtract(src_f,dest_f,sub_f);
	add(src_f,sub_f,dest_f);

	dest_f.convertTo(dest,CV_8U);
}

 

(1)ボケを表すスパース行列を生成します.

二次元のSxSの要素を持つスパース行列を確保し,必要な要素に値を入れていきます.データのアクセスには,ref,find,valueが使え,refはその要素がゼロなら初めてアクセスするときにメモリを確保します.findはゼロならNULLを返しますメモリの確保は行いません.valueは参照するだけで,メモリの確保は行いません.

(2)スパース行列を可視化します.

(1)で埋めたスパース行列を可視化します.イテレータで各行列の要素を対応するぼかしカーネルの値で埋めていきます.イテレータでは,非ゼロのの要素を持つ行列だけにアクセスします.

(3)画像を表すベクトルに(1)の行列をかけてぼかします.

入力画像を1次元ベクトルに変換した後,(2)と同様にイテレータにてスパース行列とベクトルを掛け算します.

(4)最急降下法によるボケ除去処理を行います.

(3)の逆変換を行います.ぼけがとれた理想画像とは,理想画像にぼけを乗じたものとぼけ画像のノルムを最小にする画像です.つまり|| HX-Y||を最小にするXを求めることです.内部では,ノルムをL2として最急降下法による繰り返し演算により最小化します.具体的には,

Xn+1 = Xn – β(HT(HXn-Y))

の式を繰り返し行います.βは漸近パラメータを(小さいほどステップが小さい),nが何度目の繰り返しかを表しています.ノイズ除去の実際の用途には,L1ノルムの最小化やTV正則化などが必要になってきます.こちらを参考にして下さい.
Farsiu, S.,Robinson, D., Elad, M., Milanfar, P.”Fast and robust multiframe super resolution,” IEEETrans.ImageProcessing 13 (2004)1327–1344.

(5)比較用に一般的なアンシャープマスクフィルタにより鮮鋭化します.

スパース行列によるぼけ除去と比較するために,一般的な鮮鋭化フィルタであるアンシャープマスクで画像を鮮鋭化します.(4)より非常に高速に動作します.

(6)それぞれの比較をします.

各出力画像との入力画像との違いを表示したり各結果の画像を出力します.何かキーを押すたびに,入力画像,ぼけ画像,ぼけ除去画像,鮮鋭化画像が出力されます.

実行結果例

下に,可視化したスパース行列と入力画像,出力画像を示します.スパース行列の要素がほとんどゼロであることや,
また,各出力と入力画像をL2ノルムで比較すると,ぼけ除去した画像が入力画像に近づいていることがわかります.ここではぼかしフィルタには7×7のボックスフィルタを用いています.

ぼけ画像:9554
ボケ除去画像:8569
鮮鋭化画像:8953

可視化したスパース行列.黒が0の領域.緑が非ゼロの領域

入力画像

ぼけ画像

ぼけ除去した画像

アンシャープマスクによる鮮鋭化した画像

“スパース行列(SparseMat)の使い方:画像のぼかしとぼけ除去” への10件のコメント

  1. narutomi より:

    はじめまして。

    「(4)最急降下法によるボケ除去処理を行います.」の式で
    Hの逆行列ではなく、Hの転置行列を掛けている理由を教えていただけないでしょうか。
    Googleで検索してみたり文献などにも当たってみたのですが、今ひとつ理解できません。

    よろしくおねがいします。

    • fukushima より:

      福嶋です.

      まず,Hのぼけ行列は大きな疎行列です.

      対応する1x5の画像に5x5の疎行列を書けるような行程を手書きで書いてみてください.
      そうすると,転置がちょうど逆変換になるはずです.

      もっと簡単には,ぼけ行列を平行移動の行列だと考えてみてください.(ぼけは一般的に平行移動の線形和です)そうすると1画素右に移動するための行列を転置すると1画素左に移動する行列になるかと思います.

      この移動した値と元の値からの残差を小さくするように繰り返し処理を行っていきます.

      #なお,逆行列が求まるなら,そのまま最急降下を使わなくても,そのまま掛けるだけで求まります.

      • narutomi より:

        なんとか理解することが出来ました。
        初歩的なことにも関わらず、答えてくださり本当に助かりました。
        ありがとうございます。

  2. Kirin55 より:

    返信が遅くなってすいません。
    (CUDA 4.1がリリースされたのでOpenCVをリビルドして環境を壊してしまって確認ができなかったもので。
    まだ壊す前より実行速度が遅い気がします。戻すついでにTBBとEigenを入れたからか?)

    無事ぼけ画像が表示されました。
    有り難うございました。

    文字化けの件ですが、超解像のソースは文字化けしてました。

    • Kirin55 より:

      ぼけ画像が出たのはReleaseビルドで、Debugビルドでは真っ黒になります。

      • Kirin55 より:

        Debugビルドでぼけ画像がでなかった件ですが、blurSparseMat32f()内のMat dstvecをゼロで初期化したやったら表示されるようになりました。

        • fukushima より:

          レポートありがとうございます.修正しておきます.

          おそらく,現状のcv::SparseMatでは,解像度を変えることが苦しいと思われます.代わりに疎行列を使わないバージョンのコード(Iterative Back Projectionによる実装)を張り付けておきます.行っていることは,このページのコードと等価です.
          ガウシアンブラーのみの対応ですが,コードをすこし修正すれば,いろんなボケに対応可能なはずです.

          Michal Irani, Shmuel Peleg, “Improving resolution by image registration,” CVGIP: Graphical Models and Image Processing, Volume 53, Issue 3, Pages 231–239, May 1991.
          http://pdn.sciencedirect.com/science?_ob=MiamiImageURL&_cid=272335&_user=130270&_pii=104996529190045L&_check=y&_origin=article&_zone=toolbar&_coverDate=31-May-1991&view=c&originContentFamily=serial&wchp=dGLbVBA-zSkzV&md5=cc1930b10a723161d12a0ea010716333/1-s2.0-104996529190045L-main.pdf

          //コードここから

          #include <opencv2/opencv.hpp>
          
          #ifdef _DEBUG
          #pragma comment(lib, "opencv_imgproc232d.lib")
          #pragma comment(lib, "opencv_highgui232d.lib")
          #pragma comment(lib, "opencv_core232d.lib")
          #else
          #pragma comment(lib, "opencv_imgproc232.lib")
          #pragma comment(lib, "opencv_highgui232.lib")
          #pragma comment(lib, "opencv_core232.lib")
          #endif
          
          using namespace cv;
          using namespace std;
          
          double getPSNR(const Mat& src, const Mat& src2, const Mat& mask=Mat())
          {
          	Mat src1Y,src2Y;
          	if(src.channels()==3)cvtColor(src,src1Y,CV_BGR2GRAY);
          	else src.copyTo(src1Y);
          	if(src2.channels()==3)cvtColor(src2,src2Y,CV_BGR2GRAY);
          	else src2.copyTo(src2Y);
          
          	Mat msk;
          	if(mask.empty())msk = Mat::ones(src.size(),CV_8U);
          	else mask.copyTo(msk);
          
          	Mat sub = Mat::zeros(src.size(),CV_32F);
          
          	subtract(src1Y,src2Y,sub,msk,CV_32F);
          	multiply(sub,sub,sub);
          	const Scalar sse = sum(sub);
          	const int count = countNonZero(msk);
          
          	if(sse[0] == 0.0 || count==0)
          	{
          		return 0.0;
          	}
          	else
          	{
          		double mse =sse[0] /(double)(count);
          		return 10.0*log10((255.0*255.0)/mse);
          	}
          }
          
          void iterativeBackProjection(const Mat& src, Mat& dest, const Size ksize, const double sigma, const double lambda, const int iteration, bool isWrite=true)
          {
          	Mat srcf;
          	Mat destf;
          	Mat subf;
          	src.convertTo(srcf,CV_32FC3);
          	src.convertTo(destf,CV_32FC3);
          	Mat bdest;
          
          	double maxe = DBL_MAX;
          
          	for(int i=0;i<iteration;i++)
          	{
          //for gaussian blur. If you change this kernel, this function can deblur various kernel
          		GaussianBlur(destf,bdest,ksize,sigma);
          		subf = bdest-srcf;
          		double e = norm(subf);
          //for gaussian blur. If you change this kernel, this function can deblur various kernel
          		GaussianBlur(subf,subf,ksize,sigma);
          
          		destf -= lambda*subf;
          
          		cout<<e<<endl;
          		if(i!=0)
          		{
          			if(maxe>e)
          			{
          				maxe=e;
          			}
          			else break;
          		}
          		if(isWrite)
          		{
          		destf.convertTo(dest,CV_8UC3);
          		imwrite(format("B%03d.png",i),dest);
          		}
          	}
          	destf.convertTo(dest,CV_8UC3);
          }
          
          int main(int argc, char** argv)
          {
          	Mat src = imread("lenna.png");
          	if(src.empty())cout<<"please input valid file name\n";
          
          	Mat bsrc;
          	Size ksize = Size(7,7);
          	double sigma = 5.0;
          	GaussianBlur(src,bsrc,ksize,sigma);
          
          	cout<<"PSNR before"<<getPSNR(src,bsrc)<<endl;
          
          	Mat dest;
          
          	iterativeBackProjection(bsrc,dest,ksize,sigma,0.55,500);
          	cout<<"PSNR after "<<": "<<getPSNR(src,dest)<<endl;
          }
          
          • Kirin55 より:

            新しいサンプルありがとうございます。

            こちらをベースに色々試してみます。

  3. Kirin55 より:

    はじめまして。
    このサンプルを実行してみたのですが、うまく動作しませんでした。
    512*512のlena画像を読込ませると例外が発生したので、サイズを256*256に変更したところ最後まで動作しました。しかし、入力画像と可視化したスパース行列しか表示されず、後の画像は真っ黒でした。

    OpenCv 2.3.1では、動作しないのでしょうか?もしそうであれば修正箇所を教えていただけませんか?
    よろしくお願いします。

    動作環境
    Windows7 X86
    VS2008 SP1
    OpenCV 2.3.1

    • fukushima より:

      OpenCV 2.3.1で実行してみましたが,同じように512×512の環境は実行できないようです.
      (1年以上前に作成した時は,動作していたため恐らくどこか仕様変更したのでしょう.)

      カーネルサイズなどのパラメータを変えてみたところ,小さなカーネルサイズなら実行可能でした.
      恐らくOpenCVで実装されている疎行列のサイズに起因する問題だと推測されます.
      疎行列処理をEigenという行列ライブラリで作り直したほうがよいかもしれません.

      あと,ソースコード中でだとd1,s1などが表示されていますが,(私の環境だけかもしれませんが)
      正しくはd[ c ],s[ c ]です.

      文字化けしないように命名方法を変えたソースに更新しておきました.

コメントをどうぞ