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

02 7月 2010 Under: opencv2.x-samples

Y=HX

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

# ``` #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); }   ```

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

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

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

（３）の逆変換を行います．ぼけがとれた理想画像とは，理想画像にぼけを乗じたものとぼけ画像のノルムを最小にする画像です．つまり|| 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)比較用に一般的なアンシャープマスクフィルタにより鮮鋭化します．

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

## 実行結果例

また，各出力と入力画像をL2ノルムで比較すると，ぼけ除去した画像が入力画像に近づいていることがわかります．ここではぼかしフィルタには7×7のボックスフィルタを用いています．

ぼけ画像：9554
ボケ除去画像：8569

ぼけ画像

ぼけ除去した画像

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

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

1. narutomi より:

はじめまして。

「(4)最急降下法によるボケ除去処理を行います．」の式で
Hの逆行列ではなく、Hの転置行列を掛けている理由を教えていただけないでしょうか。

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

• fukushima より:

福嶋です．

まず，Ｈのぼけ行列は大きな疎行列です．

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

もっと簡単には，ぼけ行列を平行移動の行列だと考えてみてください．（ぼけは一般的に平行移動の線形和です）そうすると１画素右に移動するための行列を転置すると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;

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 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 ]です．

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