カメラ画像のボケぐあいであるPSF(Point Spread Function、点拡がり関数)を求める１つの手法をまとめる。 Gnuplotによる区分関数の最小二乗法 を用いるとプログラミングの必要なく非常に簡単にPSFを求めることができる。 ここで仮定するPSFは、次の式で表現されるピルボックス関数（Pillbox function）である。
PSFとしてピルボックス関数を仮定する
\[ R(x, y, r)= \left\{ \begin{array}{cl} \frac{1}{\pi r^2} & x^2+y^2 \leq r^2\\ 0 & {\rm otherwise} \end{array} \right. \]
階段関数のテストパターン
テストパターンとして次の図のような白黒2色のパターンを撮影する。 これは水平方向に見れば階段関数（step function, staircase function）と見なすことができる。
\[ S(x,y)= \left\{ \begin{array}{cl} 1 & (x \leq 0)\\ 0 & (x < 0) \end{array} \right. \]
観測パターンとそのモデル
上記のパターンがピンぼけて下図（原画の4倍スケール）のようなパターンが観測される。
ピルボックス関数とテストパターンの畳み込み積分を考えると、 ピルボックス関数の一部がテストパターンの黒い領域に重なり、その残りの部分が白い領域に重なる瞬間がある。 この２つ重なり領域のピルボックス関数の体積比が画像の輝度として観測されていると考えることができる。 こう考えると観測画像の水平方向の輝度は下式で表されることが知られている。 $-1 \leq t \leq 1$の範囲がピルボックス関数とテストパターンの黒白境界が重なる範囲である。
\[ f(x)= \left\{ \begin{array}{ll} I_{\rm min} & t<-1 \\ g(t)(I_{\rm max} - I_{\rm min}) + I_{\rm min} & -1 \leq t \leq 1 \\ I_{\rm max} & t>1 \end{array} \right. \] \[ g(t) = \frac{1}{2} + \frac{1}{\pi}(t \sqrt{1-t^2}+\sin^{-1}t) \] \[ t = \frac{x - c}{r} \]PSFのパラメータ推定
上記関数のパラメータは何か？ パラメータは白い領域の輝度$I_{\rm max}$、黒い領域の輝度$I_{\rm min}$、中心位置$c$、ボケ幅$r$だけど、 PSFとして意味のあるパラメータは$r$だけである。 4つのパラメータを求めてそのうちボケ幅$r$だけがピルボックス関数の形を決めるパラメータとなる。
一方データとしては、画像の水平方向に走査したときの輝度値を用いる。 垂直方向に輝度値を平均すれば画像のノイズの影響を軽減することができる。 このデータに対して先ほどの関数をGnuplotを用いてフィッティングする。 フィッティングされた関数とデータの例を下に示す。
reset # 関数の定義 t(x,c,r)=(x-c)/r g(x,c,r)=0.5+(t(x,c,r)*sqrt(1-t(x,c,r)**2)+asin(t(x,c,r)))/pi f(x,c,r,imax,imin)=(t(x,c,r)<-1) ? imin :((t(x,c,r)>1) ? imax : g(x,c,r)*(imax-imin)+imin) # 初期値の設定 c=50 r=25 imax=95 imin=5 # フィッティング fit f(x,c,r,imax,imin) 'sample.txt' via c,r,imax,imin # データと関数描画 set xr[-5:105] plot 'sample.txt' with points pt 3 replot f(x,c,r,imax,imin) linewidth 3
Gnuplotの出力は以下の通り。 以下の結果からボケ幅は$r=21.298$画素であることが分かる。
Iteration 0
WSSR : 587.799 delta(WSSR)/WSSR : 0
delta(WSSR) : 0 limit for stopping : 1e-005
lambda : 0.859572
initial set of free parameter values
c = 50
r = 25
imax = 95
imin = 5
/
Iteration 1
WSSR : 68.9751 delta(WSSR)/WSSR : -7.5219
delta(WSSR) : -518.824 limit for stopping : 1e-005
lambda : 0.0859572
resultant parameter values
c = 49.3057
r = 21.7139
imax = 94.3184
imin = 6.75013
/
：
（省略）
：
Iteration 5
WSSR : 63.5581 delta(WSSR)/WSSR : -8.02681e-007
delta(WSSR) : -5.10168e-005 limit for stopping : 1e-005
lambda : 8.59572e-006
resultant parameter values
c = 49.278
r = 21.298
imax = 94.2626
imin = 6.80794
After 5 iterations the fit converged.
final sum of squares of residuals : 63.5581
rel. change during last iteration : -8.02681e-007
degrees of freedom (FIT_NDF) : 96
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 0.813673
variance of residuals (reduced chisquare) = WSSR/ndf : 0.662063
Final set of parameters Asymptotic Standard Error
======================= ==========================
c = 49.278 +/- 0.07501 (0.1522%)
r = 21.298 +/- 0.1477 (0.6935%)
imax = 94.2626 +/- 0.1459 (0.1548%)
imin = 6.80794 +/- 0.1517 (2.228%)
correlation matrix of the fit parameters:
c r imax imin
c 1.000
r -0.013 1.000
imax 0.432 0.325 1.000
imin 0.449 -0.338 -0.013 1.000
参考
- Gnuplotによる区分関数の最小二乗法
http://buaiso.blogspot.jp/2013/11/gnuplot.html - B. Okumura, M. Kanbara, and N. Yokoya: "Augmented reality based on estimation of defocusing and motion blurring from captured images", Proc. IEEE and ACM Int. Sympo. on Mixed Augmented Reality (ISMAR 06), pp. 219-225, 2006.
