ラプラシアン・フィルタをVerilog HDLで実装する

FPGA(Nexys4 DDR)を使ったリアルタイムエッジ検出 - メモ置き場
メディアン・フィルタに続きラプラシアンフィルタの解説．
この記事で作成したプロジェクトではラプラシアン・フィルタを使ってエッジ検出をしている*1．ラプラシアン・フィルタの原理とVerilog HDLでの実装について説明する．

ラプラシアン・フィルタの原理

ラプラシアン・フィルタは，注目画素とその周囲の画素の値を使って輪郭(エッジ)の検出を行うフィルタである．ラプラシアンという微分演算子の名前が示すように，画素の微分(2次微分)を計算する．画素の値が大きく変わっている(=何らかのエッジがある)部分では微分の値が大きくなる．よって微分の値と事前に設定した閾値を比較することでエッジの有無を判定できる*2．

OpenCVを使ってラプラシアン・フィルタを施して見た．
f:id:okchan08:20190119135816j:plain:w200

画像では空間方向の情報は画素として離散化されているため，微分は差分として表す．位置 $(x,y)$ の画素の値を $I(x,y)$ と書く． $x$ 方向と $y$ 方向への微分は
\begin{align}
I_{x}(x,y) := \frac{\partial I(x, y)}{\partial x} &= \frac{I(x+\Delta x, y) - I(x,y)}{\Delta x} = I(x+1, y) - I(x,y) \\
I_{y}(x,y) := \frac{\partial I(x, y)}{\partial y} &= I(x,y+1) - I(x,y)
\end{align}
となる．ただし，空間方向のステップ $\Delta x, \Delta y$ を1としている．同様に2次微分を中心差分で計算すると
\begin{align}
I_{xx}(x,y) := \frac{\partial^2 I(x, y)}{\partial x^2} &= I_{x}(x, y) - I_{x}(x-1,y) = \{I(x+1,y) - I(x,y)\} - \{I(x,y) - I(x-1,y) \} \\
&= I(x+1,y) - 2I(x,y) + (x-1,y) \\
I_{xx}(x,y) := \frac{\partial^2 I(x, y)}{\partial x^2} &= I(x,y+1) - 2(x,y) + I(x,y-1) \\
\nabla I(x,y) = I_{xx}(x,y) + I_{yy}(x,y) &= I(x+1,y) + (x-1,y) + I(x,y+1) + I(x,y-1) - 4I(x,y)
\end{align}
のようになる．したがってカーネル $K$ を
\begin{align}
K =
\left(
\begin{array}{ccc}
0 & 1 & 0 \\
1 & -4 & 1 \\
0 & 1 & 0
\end{array}
\right)
\end{align}
と定義して画像 $I(x,y)$ と畳み込みを計算すればラプラシアンによる2次微分が計算できる．これは中心画素を $(x,y)$ とその周囲8個の画素を使った3×3のフィルタ処理になる．

初めは3×3のラプラシアン・フィルタでエッジ検出を行っていたが，ノイズがひどいのかうまく動作しなかった．そこでフィルタサイズを一回り大きくし5×5のラプラシアンフィルタで実装することにした．5×5 ラプラシアン・フィルタのカーネルは差分の取り方によっていくつか種類があるようだが，今回は
\begin{align}
K =
\left(
\begin{array}{ccccc}
-1 & -3 & -4 & -3 & -1 \\
-3 & 0 & 6 & 0 & -3 \\
-4 & 6 & 20 & 6 & -4 \\
-3 & 0 & 6 & 0 & -3 \\
-1 & -3 & -4 & -3 & -1 \\
\end{array}
\right)
\end{align}
を使った*3．

Verilog HDLでの実装

実際のVerilog HDLファイルはこちら．
上記で説明したように，ある注目画素でラプラシアン・フィルタを計算するには25個の画素が必要となる．各画素は，BRAMで実装されたバッファに保持されている．バッファからは1クロックごとに1個の画素しか読み出せないため，画素25個を読み出すには25クロックが必要となる．画素25個を1個ずつ読み出し，読み出しが終わったら畳み込みの計算を行う，といったステートマシンを作成した．例えば12個目の画素を取り出すときは

always@(posedge CLK) begin
case(state)
    ...
    DATA12 : begin
        rd_addr <= addr_out;
        tmp_data12 <= data_in;
        state <= DATA13;
    end
    DATA13 : begin
    ...
endcase
end

のように12個目の画素が保存されたいるBRAMのアドレスをrd_addrに書き込み，データをtmp_data12に保存する．次に13個目の画素を読みに行くために，stateへは次のステートを表すDATA13を代入する．画素の番号は，25個の画素を5×5に並べた時に一番左上を0個目，その右となりを1個目，…と数えている．
画素の読み込みが終わったら畳み込み計算を行う．畳み込み計算では整数の掛け算と加減しかないので単純に

calc_data <= (tmp_data7 + tmp_data11 + tmp_data13 + tmp_data17)*6 + 20*tmp_data14 - (tmp_data2 + tmp_data10 + tmp_data14 + tmp_data22)*4
                       -(tmp_data1 + tmp_data3 + tmp_data5 + tmp_data9 + tmp_data15 + tmp_data19 + tmp_data21 + tmp_data23)*3
                       - tmp_data0 - tmp_data4 - tmp_data20 - tmp_data24;

とする．ただしcalc_dataのビット幅はオーバーフローしないように注意する必要がある．プロジェクトでは画素データは12ビットとしているので，畳み込み計算後の値は16ビット幅で表現できる範囲におさまる．したがってcalc_dataは16ビット幅のレジスタとしている．

最終的な注目画素の値は計算したcalc_dataを閾値と比較し決定する．ここでは $16384 = 2^{14}$ を閾値とした．

assign calc_edge = (|calc_data[15:13]) ? 12'hfff : 12'h0;

calc_dataの14ビット目より上のビットが立っていたら $2^{14}$ より大きいので，0xFFFを出力する．

動作確認

実装したラプラシアン・フィルタをFPGA上で動作させた．
わりと綺麗にエッジが見えている．

www.youtube.com

*1:はじめはソーベル・フィルタを使う予定だったのでプロジェクト名はsobelとなっている

*2:【画像処理】ラプラシアンフィルタの原理・特徴・計算式 | アルゴリズム雑記

*3:これはどこかの文献から引用した値なのだが，引用元をメモするのを忘れてしまった