Separable Circular Filter を調査してみた

こんにちは、研究開発室の川口です

少し前に SNS で Circular Filter を Separable にできるという投稿があり気になったので検証してみました。
この記事では Separable Circular Filter について GDC 2018年の EA の Garcia 氏による講演：Circular Separable Convolution Depth of Field “Circular Dofの発表（と同氏の SIGGRAPH 2017 の talk）、およびそのベースとなる理論の説明をしているNiemitalo 氏のブログ記事：Circularly symmetric convolution and lens blurを元に、実装検証を行なった内容の紹介を行います。

0. はじめに

画像処理におけるフィルタにはいろいろな種類があります。
例えば画像の各ピクセルについて、その周辺の正方形領域の平均値を求める Box Filter、ガウス分布に基づいた加重平均を行う Gaussian Filter などがあげられます。

この記事で取り上げる Circular Filter とはピクセル周辺の円形領域の平均値をとるフィルタを指し、Disk Filter ともいいます。

画像処理フィルタはレンダリングにおいてポストプロセスなどに利用され、ボケ表現やノイズ除去など様々な場面で活用されています。
しかし画像処理は決して軽い処理ではなく、ゲームをはじめとするリアルタイム処理では広い範囲にフィルタを適用することはパフォーマンス上簡単ではありません。

Separable とは画像処理フィルタにおいて、広い範囲を効率的フィルタリングできる特性を指します。
この記事では Separable とは何なのか、Circular Filter を Separable にするとはどういうことなのか、説明してその結果を確認してみます。

1. Separable Filter

画像処理フィルタにおける Separable とは、直訳では分離可能、つまり 2D のフィルタを 1D のフィルタ2回に分けて処理できることを示します。
つまり、縦横同時にフィルタリングした結果と、横方向にフィルタリングした結果を縦方向にフィルタリングした結果が一致するようなフィルタのことです。（縦横は順不同です）

単純な Box Filter の Compute Shader 実装のようなものを考えると、

float4 sum = 0.0;
for(int y=-Size; y<=Size; y++)
{
    for(int x=-Size; x<=Size; x++)
    {
        int2 p = pixel - int2(x, y);
        sum += image[p];
    }
}
result[pixel] = sum / ((2 * Size + 1) * (2 * Size + 1));

と記述される 2D のフィルタ処理を次のように、

// first pass
float4 sum = 0.0;
for(int x=-Size; x<=Size; x++)
{
    int2 p = pixel - int2(x, 0);
    sum += image[p];
}
result[pixel] = sum / (2 * Size + 1);

// --------------------------------------------------

// second pass
float4 sum = 0.0;
for(int y=-Size; y<=Size; y++)
{
    int2 p = pixel - int2(0, y);
    sum += image[p];
}
result[pixel] = sum / (2 * Size + 1);

2回の処理で記述できることを意味します。
このとき、first pass と second pass は異なるレンダリングパスとして実行しており、first pass の result が second pass の image になっています。

2D の naive 実装ではテクスチャのアクセスが ((2 * Size + 1) * (2 * Size + 1)) 回と２乗のオーダーで増加するので、ちょっとフィルタ範囲を広げるだけで非常に処理が重くなります。
e.g., Size = {1,2,3,4} → テクスチャアクセス回数 = {3×3 =9, 5×5 =25, 7×7 =49, 9×9 =81}

Separable Filter はパスを分ける分のコストはありますが、テクスチャのアクセスは 2 * (2 * Size + 1) で、線形に増えるので広い範囲をフィルタしたい時に非常に効率的になります。
e.g., Size = {1,2,3,4} → テクスチャアクセス回数 = {3+3 =6, 5+5 =10, 7+7 =14, 9+9 =18}

上記の処理をシェーダで実装して実行した結果も見てみます。
実際の実装は比較検証のための処理などを加えており最適化は行っていません。
Separable 処理はフィルタの方向を指定するパラメータを用意してそれを切り替えながら同じシェーダを2回 Dispatch することで実行しています。

RTX 4080 でフル HD の画像に対して実行した際の処理時間です。
表の括弧の中はテクスチャへのアクセス回数を示します

パスを2回に分けることで固定コストが増えるため Size が小さいときは Naive 実装の方が優位ですが、Size=4 からは Separable のほうが圧倒的に高速になっています。
リアルタイムレンダリングにおける大きなぼかしは、ダウンサンプリングなどを使って最適化するので極端に大きなフィルタサイズを直接使う機会は少ないですが効果は絶大です。
結果画像は省略しますが見た目は全く同じになります。

どんなフィルタでも Separable にできるというわけではなく 2D の処理と 1D の処理2回が等価である必要があります。
Box Filter のほかに Gaussian Filter も Separable であり効率的なフィルタリングとしてよく利用されています。

余談

ちなみに私は Separable という用語を Screen-Space Separable Sub-Surface Scattering (SSSSSS) というユニークな名前の手法で覚えました。
（Screen-Space はわざわざ言わなくても良いので遊び心の強い呼び方ですね）
Separable Bilateral Filter も画像処理でよく聞くと思います。
どちらも Gaussian Filter が Separable なことで成り立つフィルタです。

2. Separable の数学的な説明

Separable Filter が数学的にどうなっているかを考えましょう。

ある2次元のフィルタカーネルを $F(u,v)$ と表します。
カーネル中心の座標系 $u,v$ についての重み関数です。

このとき Separable とは、1次元のフィルタ $G(u)$ を使って、$F(u,v)=G(u)G(v)$ と書くことができるフィルタを言います。

画像処理のフィルタリングは、画像とフィルタカーネルの畳み込み演算で記述されます。
入力画像 $x,y$ における画素値を $I(x,y)$ としたとき、
$$
I_{\text{filtered}}(x,y) = \int\int I(x-u, y-v)F(u,v) \mathrm{d}u \mathrm{d}v
$$
となります。

このときフィルタカーネルが Separable であるならば、
$$
I_{\text{filtered}}(x,y) = \int\int I(x-u, y-v)G(u)G(v) \mathrm{d}u \mathrm{d}v
$$
であり、
$$
I_{\text{filtered}}(x,y) = \int\left[\int I(x-u, y-v)G(u) \mathrm{d}u \right] G(v) \mathrm{d}v
$$
と二重積分を $u$ と $v$ で2回の積分に分けて記述できます。
ここでは内側の積分が横方向だけのフィルタ、外側の積分が縦方向のフィルタを表します。
画像は離散信号なので積分を総和に置き換えて書くと

$$
\begin{align}
I_{\text{filtered}}[x,y]&=\sum\sum I[x-u, y-v]G[u,v]\\
&=\sum\sum I[x-u,y-v] G[u]G[v]\\
&=\sum\left[\sum I[x-u,y-v] G[u]\right]G[v]
\end{align}
$$
です。

次に具体的な例として Gaussian Filter を見てみましょう。
2D の Gaussian Filter は、
$$
F(u,v)=e^{-a(u^2+v^2)}
$$
です。
指数法則から
$$
F(u,v)=e^{-au^2}e^{-av^2}
$$
となり、1D の Gaussian Filter $G(u)=e^{-au^2}$ を使うと、
$$
F(u,v)=G(u)G(v)
$$ と記述できるため Separable であることが分かります。

またこのとき
$$
F(u,v)=F(r), r=\sqrt{u^2+v^2}
$$
と記述できるフィルタを円対称フィルタといます。
円対称な2Dフィルタカーネルはカーネル中心を通るどの 1D 断面でも同じ形状を取ります（中心からの距離だけに依存する関数）。
2次元ですが半径方向に見た1次元の特性だけを考えればよいので扱いやすいです。

3. Separable Circular Filter

Circular Filter は Disk Filterともいい指定された半径内で一様に円対称なフィルタを指します。
言葉で見ると Gaussian Filter に似ていますが、Gaussian はフィルタ中心に重みが円形の山のようになっているのに対して、Circular はフィルタ中心にある半径の中で一様に同じ重みを持つものを言います。
また Gaussian はフィルタ内のどこでも重みがありますが、Circular は円の外は一切重みがありません。
1D のグラフで重みを見ると以下のようになっています。

3.1. 導出

今欲しい Circular Filter は以下のフィルタカーネルとして定義されます。

$$
C(u,v)=
\begin{cases}
1&(u^2+v^2\lt R^2)\\
0&\text{othrewise}
\end{cases}
$$
このカーネルは Separable ではありません。
今回の手法では複素ガウシアンを使うことでこのカーネルを近似します。

まず複素数の 1D Gaussian Filter を考えます。通常のガウシアンの係数を複素数にしたものです。
$$
G(u)=e^{(-a+ib)u^2}
$$

オイラーの公式より
$$
G(u)=e^{(-a+ib)u^2}
=e^{-au^2}\left[\cos(bu^2)+i\sin(bu^2)\right]
$$
です。
2方向でフィルタをかけてみます。
$$
G(u)G(v)=e^{(-a+ib)u^2}e^{(-a+ib)v^2}=e^{(-a+ib)(u^2+v^2)}
$$
これは $r=\sqrt{u^2+v^2}$ を使って書くことができ
$$
G(r)=e^{(-a+ib)r^2}=e^{-ar^2}\left[\cos(br^2)+i\sin(br^2)\right]
$$
となります。
縦横にフィルタをかけた結果が、$r$ についての 1D 複素ガウシアンの式になっており、Separable かつ円対称であることが確認できました。

式の実部と虚部をカーネル形状を中心を通る断面で考えると Gaussian envelop（中心から外側に向かって減衰する形）を持ち振動しています。
フーリエ級数展開が波を足し合わせて任意の関数を表現できるように、この複素ガウシアンに適当な係数をかけて足し合わせると任意の形状を近似できそうに見えます。
この実部と虚部を係数を付けて足し合わせます。
$$
F(r) = A \ \Re(G(r)) + B \ \Im(G(r))
=e^{-ar^2}\left[A \ \cos(br^2) + B \ \sin(br^2) \right]
$$
これを一つのコンポーネントとして複数足し合わせます。

$$
\begin{align}
C(r) &\sim \sum_j F_j(r)\\
&=\sum_j\left[A_j \ \Re(G(r)) + B_j \ \Im(G(r))\right]\\
&=\sum_j e^{-a_jr^2}\left[A_j \ \cos(b_jr^2) + B_j \ \sin(b_jr^2) \right]
\end{align}
$$
1コンポーネントあたり4つの係数（$a_j,b_j,A_j,B_j$）を調整しながら複数を足し合わせていくことで狙いの形状を作ることができます。
これが Separable Circular Filter の定式化です。

係数をどう決めるのかは目的の形状に対する最適化問題になります。詳しくは元ブログ記事を参照ください。
実践的にはブログやGDC資料に記載されている Separable Circular Filter の為に作られた係数を使うことになるので、ここでは最適化問題にまでは踏み込みません。

[補足1] 正確に言うとこれは Circular Filter の定式化ではなく複素ガウシアンを足し合わせてカーネル形状を作れるということの説明です。
複数の複素ガウシアンの実部と虚部の足し合わせが任意の関数を本当に近似できるのかはもう少し議論が必要だと思います。
元のブログの中では
> Arbitrary circularly symmetric shapes can be constructed by a weighted sum of the real parts and imaginary parts of complex phasors of different frequencies.
とあり、円対称な任意の形状を構築できると言っていますが具体的な証明は書かれていません（式から明らかなのかもしれませんが）。
確かに Gaussian envelop を除くとほぼそのままフーリエ級数になりますし、envelop が円対称性のある形状の近似を阻害することはなさそうなので正しそうではあります。
私にはさらに踏み込んだ話は重かったので、ぜひ詳しい方に数学的な議論を行ってほしいです。

[補足2] 近似式 C(r) は Separable ではありません。
複素ガウシアンの実部と虚部を複数コンポーネントにわたって足し合わせる計算は Separable な性質を消しているからです。
（C(r) = Σj Fj(r) から C(r)=g(u)g(v) という分解はできない。）
なので Separable Circular Filter という用語は Separable な Circular Filter ではなく、Separable なフィルタの組み合わせで Circular Filter を近似しているという意味です。
結局 Separable な処理で Circular Filter を行うという目的は達成しており、言葉の使い方が気になる人向けの注意なのであまり気にしなくて大丈夫です。

改めて Separable なフィルタを記述してみます。
1次元の複素ガウシアンフィルタは係数も併せて書くと、
$$
\begin{align}
G(u;a,b)&=G_R(u;a,b)+iG_I(u;a,b),\\
G_R(u;a,b)&=e^{-au^2}\cos(bu^2),\\
G_I(u;a,b)&=e^{-au^2}\sin(bu^2)
\end{align}
$$
です。

Separable に縦横 $u,v$ についてかけると、

$$
\begin{align}
G(u;a,b)G(v;a,b)&=(G_R(u;a,b)+iG_I(u;a,b))(G_R(v;a,b)+iG_I(v;a,b))\\
&=\Re(G(u)G(v))+i\Im(G(u)G(v))\\
\Re(G(u)G(v))&=G_R(u)G_R(v)-G_I(u)G_I(v),\\
\Im(G(u)G(v))&=G_R(u)G_I(v)+G_I(u)G_R(v)
\end{align}
$$
となります。
最後のフィルタの出力では実部と虚部を足し合わます。

$$
F(u,v;a,b,A,B)=A\ \Re(G(u;a,b)G(v;a,b))+B\ \Im(G(u;a,b)G(v;a,b))
$$

文字だらけで非常に読みにくいですが、横方向だけのフィルタ $G(u)$ を画像にかけて、その実部 $G_R(u)$ と 虚部 $G_I(u)$ を持っておき、
それをさらに縦方向のフィルタにかけて一つの出力を作るという流れを数式で書いているだけです。

加えて、必要な近似の精度によってコンポーネントを複数足し合わせます。
$$
\begin{align}
C(u,v)&\sim\sum_j\left(F(u,v;a_j,b_j,A_j,B_j)\right)\\
&=\sum_j\left[A_j\ \Re(G(u;a_j,b_j)G(v;a_j,b_j))+B_j\ \Im(G(u;a_j,b_j)G(v;a_j,b_j)) \right]
\end{align}
$$

3.2. 画像形式

画像としての定式化も書いておきます。
実装上は画像との離散畳み込みになるので、まず実部虚部について横方向に画像をフィルタリングした結果 $H[x,y]=H_R[x,y]+iH_I[x,y]$ を作ります。
$$
\begin{align}
H[x,y] &= \sum_uI[x-u,y]G[u]\\
H_R[x,y] &= \sum_u I[x-u,y]G_R[u]\\
H_I[x,y] &= \sum_u I[x-u,y]G_I[u]
\end{align}
$$
次に縦方向のフィルタをかけます。横方向のフィルタ結果は実部虚部が分かれた画像として入力されます。
$$
\begin{align}
V[x,y]&=\sum_vH[x,y-v]G[v]\\
&=\sum_v\left(H_R[x,y-v]+iH_I[x,y-v])(G_R[v]+iG_I[v])\right)\\
&=\sum_v\left(H_R[x,y-v]G_R[v]-H_I[x,y-u]G_I[v]\right)\\
&\quad+i\left(H_R[x,y-v]G_I[v]+H_I[x,y-v]G_R[v]))\right)
\end{align}
$$
この結果 $V$ は実部虚部を分けて画像として出力する必要は無くフィルタ結果を作るための中間値です。
最終的に出力する画像は係数を付けた足し合わせになります。
$$
\begin{align}
I_{\text{filtered}}[x,y]&=A\ \Re(V[x,y])+B\ \Im(V[x,y])\\
&=\sum_v\left(H_R[x,y-v] (AG_R[v]+BG_I[v]\right)\\
&\quad+H_I[x,y-v]\left(BG_R[v]-AG_I[v])\right)
\end{align}
$$
実装上はカーネル部分を計算した結果が重み配列になっているので $H$ で括った式で書いています。

複数コンポーネントのときは画像の足し合わせになります。
$$
\text{result_image}[x,y]=\sum_j I_{\text{filtered}}^j[x,y]=\sum_jA_j\ \Re(V[x,y])+B_j\ \Im(V[x,y])
$$
2パスの Separable 処理をコンポーネントの数分だけ実行して、その結果の画像を足し合わせていきます。
実装上は2パス目の処理の結果画像を共有して足し合わせていくような感じになると思います。

3.3. 別形式

ところで、元の複素ガウシアンには $a,b$ しか係数が存在しないのに、実部と虚部に分けて書いたときに急に $A,B$ がでてきて違和感はないでしょうか。
基底展開と考えると元の式は基底を示しているだけで、展開する係数として $A,B$ を導入することはおかしくありません。
（数式に起こすときに、同じ文字を使うのかとか、$=$で結んで良いのかとか、かなりあやふやに書いてしまいました）

元のブログ記事では、update としてそのあたりを整理した定式化も書かれています。
1Dの複素ガウシアンに係数 $c+id$ を導入します。
$$
\tilde{G}(u)=(c+id)\ e^{(-a+ib)u^2}
$$
こうすることで
$$
\tilde{G}(u)=e^{-au^2}(c\ \cos(bu^2)-d\ \sin(bu^2))+i\left[ e^{-au^2}(d\ \cos(bu^2)-d\ \sin(bu^2)) \right]
$$
となり、
2D では
$$
\tilde{G}(u)\tilde{G}(v)=(c+id)^2e^{(-a+ib)(u^2+v^2)}
$$
です。
ここで $(c+id)^2=A-iB$ と前の係数で書き換えてみると、
$$
\begin{align}
\tilde{G}(r)&=(A-iB)e^{(-a+ib)r^2}=Ae^{(-a+ib)r^2}-iBe^{(-a+ib)r^2}\\
&=Ae^{-ar^2}\left[ \cos(br^2)+i\sin(br^2) \right]-iBe^{-ar^2}\left[ \cos(br^2)+i\sin(br^2) \right]\\
&=\left[Ae^{-ar^2}\cos(br^2)+Be^{-ar^2}\sin(br^2)\right]\\
&\quad+i\left[ Ae^{-ar^2}\sin(br^2) – Be^{-ar^2}\cos(br^2) \right]
\end{align}
$$
となります。
この実部は先ほど Circular Filter の近似として示した1コンポーネントの式に一致しています。
$$
F(r) = A \ \Re(G(r)) + B \ \Im(G(r))
=e^{-ar^2}\left[A \ \cos(br^2) + B \ \sin(br^2) \right]
$$

つまりこの係数付きの 1D 複素ガウシアンフィルタを Separable に適用して、2pass目は実部だけ出力すればよいということになります。
係数付きのフィルタカーネル $\tilde{G}(u)$ でも Circular Filter の近似を立式できました。

虚部を捨てても良いのかと不思議ですが、実部に重み $A$ を虚部に重み $B$ をそれぞれ与えて足す操作は $A-iB$ をかけてから虚部を捨てるのと同じと見なせます。
そのため複素係数を 1D のカーネル側に入れても同じ結果が得られます。

使う係数が変わるので文献に乗っている値の取扱いに注意は必要ですが、こちらの方が自然にフィルタを記述できていると思います。
元の $(a,b,A,B)$ 係数から $(a,b,c,d)$ 係数へは単純な計算で変換出来るので、元係数の値さえあれば専用の最適化を用意する必要もありません。
元ブログ内に 6-components までの係数も載っています。

GDC 2018 の講演時点ではこの update は入っていなかったようなので、今回はこちらの方法で実装してみました。

なお元のブログ記事には $(a,b,c,d)$ を使ってカーネルに負数が含まれないように最適化した例がのっています。
そうすると品質が上がるようですが今回は扱っていません。
興味があれば試してみても良いかもしれません。

4. 実装

実装は非常に単純で $\tilde{G}(u)=(c+id)\ e^{(-a+ib)u^2}$ を処理する1Dのフィルタ処理を用意するだけです。
2パス目ではコンポーネントの足し合わせが必要です。
係数はコンポーネント数によって変わるので StructuredBuffer で渡すようにしました。

シェーダ実装：

// first pass
void horizontalPass(int2 pixel)
{
    float4 realSum = 0.0;
    float4 imaginarySum = 0.0;
    for (int offset = -R; offset <= R; ++offset)
    {
        int2 p = pixel - int2(offset, 0);
        float3 value = inputImage[p];
        float2 weight = complexWeights[offset + R];
        realSum += value * weight.x;
        imaginarySum += value * weight.y;
    }
    outputReal[pixel] = realSum;
    outputImaginary[pixel] = imaginarySum;
}

// ---------------------------------------------------
// second pass
void verticalPass(int2 pixel)
{
    float4 realSum = 0.0;
    for (int offset = -R; offset <= R; ++offset)
    {
        int2 p= pixel - int2(0, offset);
        float3 horizontalReal = inputReal[p];
        float3 horizontalImaginary = inputImaginary[p];
        float2 weight = complexWeights[offset + R];
        realSum += horizontalReal * weight.x - horizontalImaginary * weight.y;
    }
    // 複数コンポーネントの場合 outputReal には前のコンポーネントの結果を入れておき、
    // 今回の結果と足し合わせる（正規化項もかける）
    float4 accumulated = outputReal[pixel];
    accumulated += realSum / kernel.normalization;
    outputReal[pixel] = accumulated;
}

＊インデックスのクランプなど一部簡略化して載せています。

• inputImage : 入力画像
• outputReal = inputReal: first pass の実部画像
• outputImaginary = inputImaginary: first pass の虚部画像
• complexWeights: フィルタ重みの StructuredBuffer（C++実装参照）
• kernel.normalization: コンポーネント合成のための正規化係数（C++実装参照）

フィルタカーネル本体はアプリケーション側（C++）で作って complexWeights（float2 の StructuredBuffer）に格納します。
係数は元のブログ記事を参照して $(a,b,A,B)$ で持ち、カーネルの重みを作る時に $(a,b,c,d)$ に変換しています。

struct Coefficient
{
    float envelope = 0.0f; // exp((-a+ib)u^2) の -a
    float phase = 0.0f; // cos,sin (bu^2) の b
    float realWeight = 0.0f; // A Re[G(r)] + B Im[G(r)] の A
    float imaginaryWeight = 0.0f; // A Re[G(r)] + B Im[G(r)] の B
};
struct Component
{
    std::vector<std::complex<float>> weights;
};
struct Kernel
{
    int radius;
    float normalization = 1.0f;
    std::vector<Component> components;
};

// 1 component の時の係数
constexpr Coefficient kOneComponent[] = {
    {-0.862325f, 1.624835f, 0.767583f, 1.862321f},
};
// 2 component の時の係数
constexpr Coefficient kTwoComponents[] = {
    {-0.886528f, 5.268909f, 0.411259f, -0.548794f},
    {-1.960518f, 1.558213f, 0.513282f, 4.561110f},
};
// 3 component の時の係数
constexpr Coefficient kThreeComponents[] = {
    {-2.176490f, 5.043495f, 1.621035f, -2.105439f},
    {-1.019306f, 9.027613f, -0.280860f, -0.162882f},
    {-2.815110f, 1.597273f, -0.366471f, 10.300301f},
};
std::vector<Coefficient> coefficientsFor(int componentCount)
{
    if (componentCount == 1)
        return {std::begin(kOneComponent), std::end(kOneComponent)};
    if (componentCount == 2)
        return {std::begin(kTwoComponents), std::end(kTwoComponents)};
    if (componentCount == 3)
        return {std::begin(kThreeComponents), std::end(kThreeComponents)};
}

Kernel makeCircularKernel(
    int radius, // フィルタ半径
    int componentCount, // コンポーネントの数 (1,2,3)
    float supportScale // transition bandwidth
){
    Kernel kernel;
    kernel.radius = static_cast<int>(std::ceil(radius * supportScale)); // 裾を含めるため半径を広げる

    for (const Coefficient coefficient : coefficientsFor(componentCount))
    {
        Component component;
        component.weights.reserve(kernel.radius * 2 + 1);

        // multiplier = c+id = sqrt(A-iB)
        // A,B 係数のときは multiplier は不要
        const std::complex<float> multiplier = std::sqrt(std::complex<float>(coefficient.realWeight, -coefficient.imaginaryWeight));

        for (int offset = -kernel.radius; offset <= kernel.radius; ++offset)
        {
            const float u = static_cast<float>(offset) / static_cast<float>(radius);
            const float uSquared = u * u;
            // exp(-a u^2)
            const float magnitude = std::exp(coefficient.envelope * uSquared);
            // b u^2
            const float angle = coefficient.phase * uSquared;
            // (c+id) e^{(-a+ib) u^2} = (c+id) * e^{-au^2} [cos(bu^2) + isin(bu^2)]
            component.weights.push_back(multiplier * magnitude * std::complex<float>(std::cos(angle), std::sin(angle)));
        }
        kernel.components.push_back(std::move(component));
    }

    // コンポーネント合成のための正規化係数（2Dフィルタ×複数コンポーネントとしてみたときの重みの総和）
    kernel.normalization = 0.0;
    for (const Component& component : kernel.components)
    {
        const std::complex<float> sum = std::accumulate(
            component.weights.begin(), component.weights.end(), std::complex<float>{});
        // sum2D=(Σ1D係数)^2
        const std::complex<float> sum2D = sum * sum;
        // A,B 係数のときは A sum2D.real() + B sum2D.imag() の総和
        kernel.normalization += sum2D.real();
    }
    return kernel;
}

＊非対応のコンポーネント数の扱いなど一部簡略化して載せています。

ここでは C++ サイドで定義した Kernel のセットアップだけ示しています。
半径とコンポーネント数を指定してカーネルを作ります。
これを元に push_constant と StructuredBuffer を作ってシェーダに渡しています。

5. 結果

実行パフォーマンスを見てみましょう。
RTX 4080 でフル HD の画像に対して実行した際の処理時間です。

Naive 実装と比較すると Size=4~8 あたりに分岐点があります。
また Size=8 まではフィルタリング以外のオーバーヘッドが支配的でコストの増え方は緩やかです。
component の数が増えると処理の回数が単純に倍加していくのが処理時間からも分かります。
ちなみに係数 $(a,b,A,B)$ と $(a,b,c,d)$ の両方ケースで試してみましたが、どちらも処理時間はほぼ同じでした。

1～3コンポーネントでのフィルタカーネルの断面形状も載せておきます。

1コンポーネントでは波が大きいですが、2,3コンポーネントではある程度平らな形が得られています。
しかし指定した半径境界で Gaussian envelop による滑らかな形状は残っています。

実際にフィルタとして使うときは、この重みを利用半径内の整数座標でサンプルした値を持つことになりますが、半径の扱いは注意が必要です。
半径の位置でトリミングして不連続な形にするか、裾を含めて滑らかな形にするかで、フィルタリングの品質に影響が出ます。
元ブログ記事ではどれだけ裾を含めるかを transition bandwidth で制御するとしており、リンギングの要因になるとしています。

元のブログ記事は transition bandwidth は0.2としており、実装ではそれに倣って指定半径を 1.2 倍して使いました。
T=ceil(1.2*R) として少しテクスチャアクセス数が増えています。テクスチャアクセス数はそのままでカーネル形状を内側に縮小するような方法もあります。

画像に適用した結果の絵も見てみましょう。

フィルタ形状を見るため 257×257 の中心に1ドットだけ白色を置いた物と、Kodak Lossless True Color Image Suite より kodim 23 (768×512) に適用した結果を乗せます

フィルタ種類	1ドット画像	kodim 23
Input
Box 2D
Gaussian 2D
Circular Naive
Separable 1 component R=16 (T=20)
Separable 2 components R=16 (T=20)
Separable 3 components R=16 (T=20)

1ドットのフィルタ結果ではカーネルのアーティファクトがよく分かります。
一方画像に適用すると目視では判別が付かないくらいになっています（差分画像も用意すればよかったですね）。
このような自然画像では十分な品質は得られていると言えそうです。
なお構造がはっきりした非自然画像ではアーティファクトが目立ちやすかったです。

ランダムな輝点にフィルタを適用した DoF 風の結果も載せます。512×512 の画像を HDR 値（0.0~100.0）で作成して出力時に適当にトーンマップを行っています。

フィルタ種類	ランダム輝点
Input
Box 2D
Gaussian
Circular Naive
Separable 1 component R=64 (T=77)
Separable 2 component R=64 (T=77)
Separable 3 component R=64 (T=77)

Separable Circular Filter の結果は Naive と比べると画質が悪く見え、アーティファクトも気になります。
高輝度値を扱うとカーネルの負数の部分に起因する悪い結果が目立ちやすくなるようです。
一方で Circular Filter は Box や Gaussian とも違うボケの形状が得られるのでその利点は確認できます。
係数の設定や bandwidth の扱いを変えたり、縮小バッファと組み合わせたりすることで、品質の向上やアーティファクトを目立たせないようにできないか検証が必要です。（実装にどこかミスがある可能性もあります）

GDC の講演の中では Low quality として 1 component、High quality として 2 components をゲーム内に採用したと言及されています。（1/4解像度で R=8）

6. まとめ

結果のカーネルそのものを見ると結構模様が出ていて、これを純粋に Circular Filter として利用していいかどうかは議論の余地があります。
一方でパフォーマンスの改善は明確で結果の品質もちゃんと理解していれば十分採用できるとも思います。

実際に GDC の講演では DoF（被写界深度） にこの Separable Circular Filter を採用していて、ゲームのスクリーンショット付きで紹介しています。
ガウシアンフィルタと比較して Circular Filter が有効なのは均一の重みと有限の領域が必要なケースで、まさに DoF は適した使用先と言えるでしょう。

しかし Separable Filter はマスクを含む（例えば、DoF の手前と焦点と奥でフィルタの有無やサイズが不連続に変化する部分）を扱う時に注意が必要です。
例えばフィルタ領域の右下 1/4 が遮蔽されている（フィルタ結果に寄与しない）とします。
2D フィルタではその領域をサンプルしないだけで良いですが、Separable Filter では寄与しない領域が1方向のフィルタで評価されるので、2方向目のフィルタでその寄与しない領域が広がってしまいます。
GDCの講演の中でもこの遮蔽によるアーティファクトについて言及されています。

DoF のために大きなフィルタを高速に処理できるのはメリットですが、カーネル形状の近似によるアーティファクトや不連続なカーネルサイズの変化に対するアーティファクトを許容できるかどうか、ユースケースによってしっかり議論したいところです。

数学的な議論も含めてこの手法はかなり興味深いので、性能と品質のトレードオフをよく考えて活用するチャンスがあればぜひ使ってみたいところです。
特に、Circular ではないような非均一な円対称のプロファイルにフィッティングして、特殊なフィルタリングとして利用するような応用には可能性を感じます。

この記事が読者の皆さんの何かの参考になれば幸いです。