スパース行列(SparseMat)の使い方：画像のぼかしとぼけ除去
画像のぼけの過程は巨大なスパース行列を使って表現できます．
入力画像x(サイズは，S=width*height)，ぼけ画像yを1次元ベクトルX,Y (画像を縦横を区別せずならべたもの)としたものをぼけ過程を表す行列Hで表すと，
Y=HX
で表せます．この行列Hは画像サイズ×画像サイズ（SxS）の巨大な行列になり，Matクラスなどで確保すると あっという間にメモリが枯渇します．幸いなことに，この行列Hは要素がほとんどゼロなため，要素が非ゼロの 場所だけメモリを使用するスパース行列クラスSparseMatが使用可能です．ここでは， SparseMatの使い方を画像のぼかしとその逆過程であるぼけ除去を例に説明します．
C++
(1)ボケを表すスパース行列を生成します．
二次元のSxSの要素を持つスパース行列を確保し，必要な要素に値を入れていきます．データのアクセスには，ref，find，valueが使え，refはその要素がゼロなら初めてアクセスするときにメモリを確保します．findはゼロならNULLを返しますメモリの確保は行いません．valueは参照するだけで，メモリの確保は行いません．
(2)スパース行列を可視化します．
（1）で埋めたスパース行列を可視化します．イテレータで各行列の要素を対応するぼかしカーネルの値で埋めていきます．イテレータでは，非ゼロのの要素を持つ行列だけにアクセスします．
(3)画像を表すベクトルに（１）の行列をかけてぼかします．
入力画像を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)比較用に一般的なアンシャープマスクフィルタにより鮮鋭化します．
スパース行列によるぼけ除去と比較するために，一般的な鮮鋭化フィルタであるアンシャープマスクで画像を鮮鋭化します．（４）より非常に高速に動作します．
(6)それぞれの比較をします．
各出力画像との入力画像との違いを表示したり各結果の画像を出力します．何かキーを押すたびに，入力画像，ぼけ画像，ぼけ除去画像，鮮鋭化画像が出力されます．
実行結果例
下に，可視化したスパース行列と入力画像，出力画像を示します．スパース行列の要素がほとんどゼロであることや，
また，各出力と入力画像をL2ノルムで比較すると，ぼけ除去した画像が入力画像に近づいていることがわかります．ここではぼかしフィルタには7×7のボックスフィルタを用いています．
可視化したスパース行列．黒が０の領域．緑が非ゼロの領域
入力画像
ぼけ画像
ぼけ除去した画像
アンシャープマスクによる鮮鋭化した画像
はじめまして。
「(4)最急降下法によるボケ除去処理を行います．」の式で
Hの逆行列ではなく、Hの転置行列を掛けている理由を教えていただけないでしょうか。
Googleで検索してみたり文献などにも当たってみたのですが、今ひとつ理解できません。
よろしくおねがいします。
福嶋です．
まず，Ｈのぼけ行列は大きな疎行列です．
対応する１ｘ５の画像に５ｘ５の疎行列を書けるような行程を手書きで書いてみてください．
そうすると，転置がちょうど逆変換になるはずです．
もっと簡単には，ぼけ行列を平行移動の行列だと考えてみてください．（ぼけは一般的に平行移動の線形和です）そうすると１画素右に移動するための行列を転置すると1画素左に移動する行列になるかと思います．
この移動した値と元の値からの残差を小さくするように繰り返し処理を行っていきます．
#なお，逆行列が求まるなら，そのまま最急降下を使わなくても，そのまま掛けるだけで求まります．
なんとか理解することが出来ました。
初歩的なことにも関わらず、答えてくださり本当に助かりました。
ありがとうございます。
返信が遅くなってすいません。
（CUDA 4.1がリリースされたのでOpenCVをリビルドして環境を壊してしまって確認ができなかったもので。
まだ壊す前より実行速度が遅い気がします。戻すついでにTBBとEigenを入れたからか？）
無事ぼけ画像が表示されました。
有り難うございました。
文字化けの件ですが、超解像のソースは文字化けしてました。
ぼけ画像が出たのはReleaseビルドで、Debugビルドでは真っ黒になります。
Debugビルドでぼけ画像がでなかった件ですが、blurSparseMat32f()内のMat dstvecをゼロで初期化したやったら表示されるようになりました。
レポートありがとうございます．修正しておきます．
おそらく，現状の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; }新しいサンプルありがとうございます。
こちらをベースに色々試してみます。
はじめまして。
このサンプルを実行してみたのですが、うまく動作しませんでした。
512*512のlena画像を読込ませると例外が発生したので、サイズを256*256に変更したところ最後まで動作しました。しかし、入力画像と可視化したスパース行列しか表示されず、後の画像は真っ黒でした。
OpenCv 2.3.1では、動作しないのでしょうか？もしそうであれば修正箇所を教えていただけませんか？
よろしくお願いします。
動作環境
Windows7 X86
VS2008 SP1
OpenCV 2.3.1
OpenCV 2.3.1で実行してみましたが，同じように512×512の環境は実行できないようです．
（1年以上前に作成した時は，動作していたため恐らくどこか仕様変更したのでしょう．）
カーネルサイズなどのパラメータを変えてみたところ，小さなカーネルサイズなら実行可能でした．
恐らくOpenCVで実装されている疎行列のサイズに起因する問題だと推測されます．
疎行列処理をEigenという行列ライブラリで作り直したほうがよいかもしれません．
あと，ソースコード中でだとd1,s1などが表示されていますが，（私の環境だけかもしれませんが）
正しくはd[ c ],s[ c ]です．
文字化けしないように命名方法を変えたソースに更新しておきました．