はじめに
気づけば以下の記事を投稿したのはもう5年も前です。
全くの無名だったKuwaharaフィルターも、今となっては3DモデリングソフトBlenderに搭載されるほど出世しました。
当時大学生だった私も今年で社会人3年目。時の流れは早いものです。
お仕事では光学を扱っていることもあり、最近は波数空間と仲良くなりました。
波数空間の積は実空間の畳み込みに対応するんですよ。面白いですよね、知ってましたか？そんなの知ってる？そっか。。。
そんなこんなでKuwaharaフィルターを波数空間の性質を活かすことで（見かけ上）ループをなくして効率的に計算できないかと考えていたところ、結果的にcv2.filter2Dだけで実装できることが判明しましたのでご紹介したいと思います。
また、この手法ではKuwaharaフィルターをより一般的に拡張したものを計算量の変化なく取り扱うことができます。
Kuwaharaフィルターの定義を振り返る
Kuwaharaフィルターは、画像の平滑化（ノイズ除去）とエッジ保存の両立を目的とした非線形フィルタで、
あるピクセルの左上・右上・左下・右下の4つのエリア（下図abcd）でそれぞれ分散を計算し、分散が最も小さいエリアの平均値をそのピクセルの値とするという演算でした。
数式で表すと、
$$
I'(x_0, y_0)=\mu_{k^*}, \quad \text{where} \quad k^*=\arg\min_k\left(\frac{1}{|R_k|}\sum_{(i, j) \in R_k} \left( I(x_0+i, y_0+j) - \mu_k \right)^2 \right)
$$
となります。ここで $\mu_k$ は領域 $k$ の平均値です。
ノイズをやエッジを含むエリア、特に注目ピクセルと異なる（輝度の差がある）領域を含む場合は分散が大きくなりがちなので、境界を飛び越えづらく、エッジが保存されるわけですね。
Kuwaharaフィルターの問題点（またあるときは利点）
Kuwaharaフィルターを使用するとブロックアーティファクトと呼ばれる、元の画像に存在しない四角っぽい形状が発生することが多々あります。
これは領域が離散的かつ四角形であることに由来する現象です。
これが筆のタッチのように見えて、絵画のような出力を生み出す要因の一つとなっているわけですが、もともとのフィルターの目的である平滑化という観点では新たなノイズを生み出しているわけで好ましくありません。
このアーティファクトは領域の形状をより滑らかに変化させたり、サンプリングする点を4領域から増やしたりすることで解消できます。
如何にして高速化するか
Kuwaharaフィルターの計算が重いのは、全てのピクセルについて周囲のピクセルの分散の計算を行い、それらを比較し採用する領域を決定し、平均を計算する...といった局所的かつ複数ステップの計算があることが原因です。
GPUを使用してこれらを並列に処理してしまえばかなりの高速化が可能ですが、Pythonでも簡単に高速に扱いたい！
これら手順のうち最も効率化が難しいのが、局所的な分散の比較、そしてそのインデックスを使用しての領域選択です。
ここで登場するのが、深層学習ではおなじみのSoftmaxです。
Softmaxは以下の式で表現される計算で、「最大値を1に近く、それ以外を0に近くする」といった出力を得ることができます。やっていることとしては結構力技で、値間の差を指数関数的に拡大して再スケールしているだけなんですが。
$$
\text{softmax}(z_i) = \frac{e^{\beta z_i}}{\sum_{j=1}^{K} e^{\beta z_j}}
$$
例えば [6, 3, 1, 1] を入力すると [0.94, 0.05, 0.01, 0.01] のような出力が得られます。
$\beta\rightarrow\infty$ でone-hotに漸近します。
これを重みとして平均値に掛けてやって足し合わせれば、実質的に分散最大のインデックス部分の平均値のみ採用したことになります。
このままだと最大インデックスになってしまうので、マイナスをかけて代入します。
つまり、
$$
I'(x_0, y_0)=\sum_k\mu_{k}\frac{e^{-\beta V_k}}{\sum_k e^{-\beta V_k}}=\frac{\sum_k\mu_{k}e^{-\beta V_k}}{\sum_k e^{-\beta V_k}}
$$
ここで $V_k$ は領域 $k$ の分散です。
これだけの計算で非局所的だった argmin が、比較計算なしで近似可能となりました。
argmin を伴わない局所的な和は移動平均フィルターで実装可能ですので、以下の4つが得られればKuwaharaフィルターの結果が計算できます。
$$
1.\quad\mu_k = \frac{1}{|R_k|}\sum_{(i, j) \in R_k} I(x_0+i, y_0+j)
$$
$$
2.\quad V_k = \left[\frac{1}{|R_k|}\sum_{(i, j) \in R_k} I^2(x_0+i, y_0+j)\right] - \mu_k^2
$$
$$
3.\quad \sum_k\mu_{k}e^{-\beta V_k}
$$
$$
4.\quad \sum_k e^{-\beta V_k}
$$
2式は $V[X]=E[(X-\mu)^2]=E[X^2]-\mu^2=E[X^2]-E[X]^2$ に由来します。
そして、これらは全て cv2.filter2D のみで計算可能です！
カーネルの任意性
もともとのKuwaharaフィルターと同等の結果を得るためには、奇数サイズの定数カーネルと、同サイズの四隅のみ1のカーネルを使用します。
例えばカーネルの1辺=3の場合、
を使用すれば良いです。
しかし、カーネル$K$を円形にすれば前述のブロックアーティファクトを回避できますし、カーネル$K_2$の値を持つ箇所を増やしてやれば候補エリア数を左上・右上・左下・右下の4エリアに限らず増やすことができます。
例えば $K, K_2$ ともに上のようなカーネルを使用すれば、円形領域内のあらゆる点を中心とした半径3のエリアが候補になります。
実装
import numpy as np
import cv2
def g_kuwahara(im, s=3, beta=30, *, k=None, k2=None):
if k is None:
assert s % 2 == 1, "s must be odd"
s //= 2
r = np.linalg.norm(np.mgrid[-s:s+1, -s:s+1], axis=0)
k = np.clip(1-(r-s), 0, 1)
k /= k.sum()
k2 = k2 if k2 is not None else k
a = cv2.filter2D(im, -1, k)
b = cv2.filter2D(im**2, -1, k)
v = b - a**2
if v.ndim == 3:
v = v.sum(axis=-1, keepdims=True)
ev = np.exp(-beta * abs(v / v.max()))
c = cv2.filter2D(a * ev, -1, k2)
d = cv2.filter2D(ev, -1, k2)
c, d = np.atleast_3d(c, d)
return np.squeeze(c / d)
グレースケール／カラー、カーネル指定ありなし等オプションを用意していたり、オーバーフローを少しでも避けるために各所に正規化処理があったりするので少し余計な行が多いですが、コアは以下のこれだけです。
a = cv2.filter2D(im, -1, k)
b = cv2.filter2D(im**2, -1, k)
v = b - a**2
ev = np.exp(-beta * abs(v / v.max()))
c = cv2.filter2D(a * ev, -1, k2)
d = cv2.filter2D(ev, -1, k2)
filtered_im = c / d
Kuwaharaフィルターを自身で実装したことのある方にしか伝わらないかもしれませんが、驚きの短さです。
ちなみにこの実装では、カーネル指定無しの場合円形カーネルを使用するようにしています。
また、なるべくone-hotになるようにデフォルトでは $\beta$ を大きく設定しているので、小さく指定するか入力画像の値が0~1でないとオーバーフローします。
使ってみる
元画像
Kuwaharaフィルター同様の結果になるカーネル
円形カーネル
絵画っぽさではやっぱりオリジナルKuwaharaフィルターが一番ですね。
円形カーネルは移動平均フィルターの階調の滑らかさ＋メディアンフィルターのエッジ保持といったバランスで、これはこれで使い道がありそうです。
おわりに
Softmaxを導入することにより、高速に実装されたcv2.filter2DでKuwaharaフィルターの計算が可能となりました。
これでも普通の畳み込みフィルターの4倍+αの計算量ですが、数百x数百ピクセル程度の画像であればリアルタイムにフィルター結果を表示するUIをPythonで構築できる程度には高速化できました。
あなたのソフトにも気軽に組み込んでみてください。