Visual Basic でOpenCV(28) - フーリエ変換

Visual Basic とOpenCVを用い、画像に対する離散フーリエ変換を行う例を紹介します。

離散フーリエ変換
CCvFunc.vb
実行

離散フーリエ変換

画像のフーリエ変換は、画像をsin、cos の要素に分解する処理です。空間ドメインから周波数ドメインへ変換を行います。画像圧縮などでは、人間の目が高周波に鈍感なことを利用し、高周波成分を削除することによってデータ量を低減します。
ここで扱うフーリエ変換およびその逆変換は、正確には離散フーリエ変換（DFT）および逆離散フーリエ変換（IDFT）です。画像に対しDFT を実施し、結果を可視化してパワースペクトルを画像として表示します。
フォームに画像ファイルをドロップすると、入力画像とDFT のパワースペクトルを表示します。

CCvFunc.vb

以降に、実際にDFT を行う、CCv の派生クラスCCvFunc を示します。

Private Function mat2Dft(src As Mat) As Mat
    Using gray = New Mat()
        Cv2.CvtColor(src, gray, ColorConversionCodes.BGR2GRAY)
        Dim srcReal As New Mat()

        Dim dftRows As Integer = Cv2.GetOptimalDFTSize(gray.Rows)
        Dim dftCols As Integer = Cv2.GetOptimalDFTSize(gray.Cols)

        ' srcRealの左隅に入力の値が格納されている、サイズを
        ' DFTへ最適し拡張された部分は0で埋める。
        Cv2.CopyMakeBorder(gray, srcReal, 0,
                           dftRows - gray.Rows, 0, dftCols - gray.Cols,
                                        BorderTypes.Constant, Scalar.All(0))

        ' real と image(0)とマージして複素画像(行列)へ
        Dim f32 As New Mat()
        srcReal.ConvertTo(f32, MatType.CV_32F)
        Dim planes As Mat() = {f32, Mat.Zeros(srcReal.Size(), MatType.CV_32F)}
        Dim complex As New Mat()
        Cv2.Merge(planes, complex)  ' complexはrealとimaginaryを持ったMat

        ' DFTの実行、結果は複素数、log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
        Dim dft As New Mat()
        Cv2.Dft(complex, dft)          ' dftはDFTの結果
        Return dft
    End Using
End Function

mat2Dft メソッドは入力画像src をDFT し、結果を返します。返すMatオブジェクトは、複素数（2 チャンネル）で型はfloat（CV_32F）です。サイズはDFT に適したサイズです。Cv2.CvtColor メソッドで入力画像src をグレイスケールに変換し、Mat オブジェクトgray へ格納します。次に、Cv2.GetOptimalDFTSize メソッドでDFT を高速に処理できるサイズをdftRows とdftCols へ求めます。求めたサイズでMatオブジェクトsrcReal を生成し、左隅に入力画像をコピーし、余白のボーダーに0 を設定します。
DFT 処理を行うためにはMat オブジェクトの型は浮動小数点、そしてチャンネル数が2 でなければなりません。そこで、srcReal オブジェクトのConvertTo メソッドで、srcRealオブジェクトの型を浮動小数点へ変換しf32 へ格納します。次に、この変換したMat オブジェクトと値が0 のMat オブジェクトで、Mat オブジェクト配列のplanes を作成します。これをCv2.Merge メソッドで2 チャンネルのMat オブジェクトcomplex へ変換します。このcomplex をCv2.Dft メソッドに与えDFT 変換を行います。
最後に、DFT 変換したMat オブジェクトdft を呼び出し元へ返します。
dft2dispMat メソッドは、DFT 処理されたMat オブジェクトを表示できるように変換します。

Private Function dft2dispMat(complex As Mat) As Mat
    Dim planes = New Mat(1) {}

    Cv2.Split(complex, planes)     ' planes(0) = real, (1) = imaginary

    ' dst(x, y) = sqrt(pow(src1(x, y), 2) + pow(src2(x, y), 2))
    Dim magnitude As New Mat()
    Cv2.Magnitude(planes(0), planes(1), magnitude) ' planes(0) = DFT magnitude

    magnitude += Scalar.All(1)         ' 表示用に各ピクセル値に1.0を加算
    Cv2.Log(magnitude, magnitude)      ' 対数へ

    Dim srcDFT As New Mat()
    Cv2.Normalize(magnitude, magnitude, 0, 1, NormTypes.MinMax)
    magnitude.ConvertTo(srcDFT, MatType.CV_8U, 255, 0) ' 表示用DFT

    Return srcDFT
End Function

本メソッドが返すMat オブジェクトは1 チャンネルで、型はCV_8Uです。引数で受け取ったMat オブジェクトcomplex にはDFT を行った結果が格納されています。complex は、実数部と虚数部の2 チャンネルを持った複素数です。まず、Cv2.Split メソッドで実数部と虚数部をplanes 配列に分離します。次にCv2.Magnitude メソッドでx およびy 配列の対応する要素から形成される2D ベクトルの大きさを計算します。x はplanes(0)、y はplanes(1) に格納されます。Cv2.Magnitude メソッドの処理を以降に式で示します。

この結果をそのまま正規化して線形スケール輝度で表示すると、わずか256 諧調であることと低周波成分と高周波成分に大きな差がある場合、高周波成分が黒で潰れてしまいます。一般的に画像は低周波の部分に集中しますので、高周波部分の観察が困難になります。このため、見やすくなるように対数スケールを採用します。なお、Cv2.Log メソッドへ与える値が1.0 以下ではマイナスの値になるため、以降の正規化で面倒が起きないように、本来の値に1.0 を加算してlog（自然対数）します。その後、Cv2.Normalize メソッドで0.0 ～ 1.0へ正規化します。最小値と最大値が0.0 ～ 1.0 とは限りませんので入力Mat オブジェクトの最小値と最大値を探したのちに正規化した方が見やすくなることが考えられます。しかし、このような方法を採用すると得られる画像の輝度が原画像に依存します。ここでは異なる画像を相対的に観察することも考慮し、0.0 ～ 1.0 の固定値を採用します。単純にひとつの画像の周波数成分を観察したいときは、先に説明した方法を採用するのも良いでしょう。最後に、この正規化したMat オブジェクトの型を、Mat.ConvertTo メソッドでCV_8U へ変更します。その際に0.0 ～ 1.0 を0 ～ 255 へスケール変換します。この結果をフォームに表示します。これらの正規化やスケールはDFT 的な意味はなく、単にパワースペクトルの視認性がよくなるように加えた処理です。

swapDft メソッドは、Mat オブジェクトの象限入れ替えを行います。

Private Function swapDft(src As Mat) As Mat
    Dim swap As Mat = src.Clone()

    Dim cx As Integer = swap.Cols \ 2     ' center
    Dim cy As Integer = swap.Rows \ 2     ' center

    Dim j1 As New Mat(swap, New Rect(0, 0, cx, cy))       ' Top-Left
    Dim j2 As New Mat(swap, New Rect(cx, 0, cx, cy))      ' Top-Right
    Dim j3 As New Mat(swap, New Rect(0, cy, cx, cy))      ' Bottom-Left
    Dim j4 As New Mat(swap, New Rect(cx, cy, cx, cy))     ' Bottom-Right

    Dim jPart As New Mat()
    j1.CopyTo(jPart)       ' swap 1 <-> 4
    j4.CopyTo(j1)
    jPart.CopyTo(j4)
    j2.CopyTo(jPart)       ' swap 2 <-> 3
    j3.CopyTo(j2)
    jPart.CopyTo(j3)

    Return swap
End Function

以降に入れ替えの様子を示します。DFT を行った場合、結果の中心部が高周波で、四隅が低周波です。一般的にパワースペクトルは中心を低周波として表現しますので、象限を入れ替えます。以降に象限を入れ替える様子を図で示します。

これまで紹介したメソッドを呼び出すのがDoCvFunction メソッドです。このメソッドはForm1 から呼び出されます。以降に、ソースリストを示します。

Public Function DoCvFunction() As Bitmap
    Dim mDFT As Mat = mat2Dft(mSrc)         ' image to DFT
    Dim dft8u As Mat = dft2dispMat(mDFT)    ' DFT to display image
    mDst = swapDft(dft8u)                   ' swap: 1 <-> 4, 2 <-> 3

    Return OpenCvSharp.Extensions.BitmapConverter.ToBitmap(mDst)
End Function

入力画像が格納されているmSrc を引数に、mat2Dft メソッドを呼び出します。これによって、入力画像をDFT した結果がmDFT に格納されます。次に、このmDFT を引数に、dft2dispMat メソッドを呼び出します。これによって、dft8u に表示形式に変換した結果が格納されます。このまま表示すると象限が入れ替わっています。そこでswapDft を呼び出し、mDst に象限を入れ替えた結果を得ます。最後に、このmDst をBitmap オブジェクトへ変換して、呼び出し元へ戻します。