スパース行列(SparseMat)の使い方:画像のぼかしとぼけ除去
画像のぼけの過程は巨大なスパース行列を使って表現できます.
入力画像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);
}
#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のボックスフィルタを用いています.