はじめに
先日，SG法というデータ平滑化の手法を紹介いたしました．
SG法をPythonで実装したらえらい簡単だった．
SG法は，等間隔で並んだデータ点の各点に対し，点の前後 $N$ 点（合計 $2N+1$ 点）を最小二乗多項式フィッティングし，その中心点を平滑化後の値とする手法です．．
この手法の最大の特徴はデータの等間隔を仮定することにより，非常に単純な計算でこれを実現できる点です．
より具体的には，カーネルサイズ $N$ に対応した重みを最初の1度のみ計算して，あとは同じ重みでデータ全体に畳み込み演算を施すだけです．
個別のサブデータセットに対して最小二乗多項式フィッティングするなど複雑な演算に感じますが，
結局のところ，ただの重み付き移動平均の一種として扱えるわけです．
ところでニューラルネットワークの流行で，畳み込みといえば2次元のイメージがどうしても付きまといますよね．
…ん？これって画像の平滑化フィルタとしても使えるのでは？
2次元のSG法を考える
点が画素値，メッシュがフィッティングした曲面↑．
例として2次の2変数関数による最小二乗法フィッティングを考えます．
関数
\begin{align}
f(x, y) &= a_0 + a_1 x + a_2 y + a_3 xy + a_4 x^2 + a_5 y^2\\
&=
\left(
\begin{array}\\
1& x& y& xy& x^2& y^2
\end{array}
\right)
\left(
\begin{array}\\
a_0\\
a_1\\
a_2\\
a_3\\
a_4\\
a_5\\
\end{array}
\right)
\end{align}
で $(2N+1)\times (2N+1)$ の格子点
\begin{align}
\left[
\begin{array}\\
(-N, -N)& (-N+1, -N)& \cdots & (0, -N) &\cdots & (N-1, -N) & (N, -N)\\
(-N, -N+1)& (-N+1, -N+1)& \cdots & (0, -N+1) &\cdots & (N-1, -N+1) & (N, -N+1)\\
\vdots&\vdots&&\vdots&&\vdots&\vdots\\
(-N, 0)&(-N+1, 0)&\cdots&(0, 0)&\cdots&(N-1, 0)&(N, 0)\\
\vdots&\vdots&&\vdots&&\vdots&\vdots\\
(-N, N-1)& (-N+1, N-1)& \cdots & (0, N-1) &\cdots & (N-1, N-1) & (N, N-1)\\
(-N, N)& (-N+1, N)& \cdots & (0, N) &\cdots & (N-1, N)& (N, N)\\
\end{array}
\right]
\end{align}
の値 $f_{x, y}$ をフィッティングします．
中心座標を $(0, 0)$ に設定する点がポイントです．
これには1次元も2次元も関係ありません．各点における $f_{x, y}$ を $(2N+1)^2\times 6$ の行列に均した以下の連立方程式を解けば良いので，
\begin{align}
\left(
\begin{array}\\
1& -N& -N& (-N)^2& (-N)^2& (-N)^2\\
1& -N+1& -N& (-N+1)(-N)& (-N+1)^2& (-N)^2\\
&\vdots&\vdots&\vdots&\vdots&\vdots&\\
1&0&0&0&0&0\\
&\vdots&\vdots&\vdots&\vdots&\vdots&\\
1& N-1& N& (N-1)N& (N-1)^2& N^2\\
1& N& N& N^2& N^2& N^2\\
\end{array}
\right)
\left(
\begin{array}\\
a_0\\
a_1\\
a_2\\
a_3\\
a_4\\
a_5\\
\end{array}
\right)
=
\left(
\begin{array}\\
f_{-N, -N}\\
f_{-N+1, -N}\\
\vdots\\
f_{0, 0}\\
\vdots\\
f_{N-1, N}\\
f_{N, N}\\
\end{array}
\right)
\end{align} \tag{1}
それぞれの行列，ベクトルを $\boldsymbol{X}, \boldsymbol{A}, \boldsymbol{F}$ とおいて，
\begin{align}
\boldsymbol{X} \boldsymbol{A} &= \boldsymbol{F}\\
\boldsymbol{X^T}\boldsymbol{X} \boldsymbol{A} &= \boldsymbol{X^T}\boldsymbol{F}\\
\boldsymbol{A} &= (\boldsymbol{X^T}\boldsymbol{X})^{-1}\boldsymbol{X^T}\boldsymbol{F}
\end{align}
で係数ベクトル $\boldsymbol{A}$ が求まります．
SG法では，このフィッティングした曲面の中心点 $(0, 0)$ での値が知りたいのでした．さて対応する値はどうなっているでしょうか．
式（1） を見ると．．．
$$a_0 = f_{0,0}$$
あら簡単．
つまり，$(\boldsymbol{X^T}\boldsymbol{X})^{-1}\boldsymbol{X^T}$ の1行目のみ $\boldsymbol{F}$ にかけてあげれば求まるのです．
あとは $(\boldsymbol{X^T}\boldsymbol{X})^{-1}\boldsymbol{X^T}$ の1行目を $(2N+1)\times (2N+1)$ の行列に戻してあげれば，
任意の $(2N+1)\times (2N+1)$ の領域を2次関数でフィッティングして中心点の値を抽出するフィルタの完成です．
$(\boldsymbol{X^T}\boldsymbol{X})^{-1}\boldsymbol{X^T}$ は $f_{x, y}$ に依存しない定数ですから任意の領域に使い回せるんですね．
それを各画素に適用，まさに畳み込みです．
$f(x,y)$ の次数を上げれば，より複雑なフィッティングも可能となります．
例えば以下は $11\times 11$ の領域を3次関数でフィッティングするフィルタです．
SG法の魅力の一つに，微分値も得ることができる点があります．
この特徴も勿論2次元でも活きていて，例えば $x$ 方向の3次フィッティング平滑化微分フィルタは
$$\left.\frac{\partial f(x, y)}{\partial x}\right|_{x=0, y=0} = a_1$$
であるから，$(\boldsymbol{X^T}\boldsymbol{X})^{-1}\boldsymbol{X^T}$ の2行目を $(2N+1)\times (2N+1)$ の行列に戻してあげれば良いことになります．
いざ実装
意気込むほど難しくはありません．
def SGFilter(m, N):
mask = (np.c_[:m+1] + np.r_[:m+1]).flatten() <= m # そのままだと2m次まで入ってしまうのでマスクを作成
X = np.zeros(((2*N+1)**2, mask.sum()))
for i, j in np.indices((2*N+1, 2*N+1)).T.reshape(-1, 2):
X[i*(2*N+1)+j] = ((i-N)**np.c_[:m+1] * (j-N)**np.r_[:m+1]).flatten()[mask]
C = np.linalg.pinv(X).reshape(-1, 2*N+1, 2*N+1) # (X^T X)^-1 X^T
return C
使ってみる
Lennaちゃん
$3$次，$7\times 7$
C = SGFilter(3, 3)
imf1 = cv2.filter2D(im, -1, C[0])
plt.imshow(imf1)
plt.show()
imf2 = cv2.filter2D(im, -1, C[1])
plt.imshow(imf2[..., 0], cmap='gray')
plt.show()
$4$次，$27\times 27$
C = SGFilter(4, 13)
imf1 = cv2.filter2D(im, -1, C[0])
plt.imshow(imf1)
plt.show()
imf2 = cv2.filter2D(im, -1, C[1])
plt.imshow(imf2[..., 0], cmap='gray')
plt.show()
結論
使えなくはないけど，フィルターサイズが大きくなると明るさの境界の影響が気になる仕上がり．
画像以外ならあるいは．
画像系のNNの初期重みとして使ったらどうなるか少し気になる．