画像のぼけの過程は巨大なスパース行列を使って表現できます．
入力画像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次元ベクトルに変換した後，（２）と同様にイテレータにてスパース行列とベクトルを掛け算します．

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

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

Xⁿ⁺¹ = Xⁿ – β(H^T(HXⁿ-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)比較用に一般的なアンシャープマスクフィルタにより鮮鋭化します．

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

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

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

実行結果例

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

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

入力画像

ぼけ画像

ぼけ除去した画像

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