iOS で物体の向きをリアルタイムに計算①【2D FFTでパワースペクトル画像生成】

物体の向きを周波数領域で判定してみます。

＜完成イメージ＞
　

これだけでは画像の中での物体位置を考慮した判定はできませんが、物体が強いエッジをもっている場合であれば判定できそうです。
この記事では上記のようにパワースペクトルを表示するまでを記載します。
向きを求める方法は次の記事に記載予定です。
（周波数領域での周波数成分が物体のエッジと直交ので、角度別扇状の成分の総和の比較で判定できると考えてます）。

サンプルコードはGithubに置いています。

パワースペクトル画像を作成するステップ

1. 画像の周囲を暗くする
2. カラー画像を8bitグレースケールに変換
3. 二次元フーリエ変換
4. 複素数を画像として扱えるように加工
5. （見た目のため）エイリアシング?!を生成
6. 画素値配列をCGImageに変換

注：目的が「向きを把握する」なので(フーリエ逆変換したりしないので)荒っぽい操作が含まれているかもしれません。

1. 画像の周囲を暗くする

画像をそのままFFTにかけるとフラットな画像であっても縦横に周波数成分が現れてしまうので、方向を検出する目的としてはうれしくありません。そこで、入力画像の周囲を暗くすることでこの影響を抑制します。

下図は白い平面を写した時のパワースペクトルです。

画像をそのままFFT	周囲を暗くしてからFFT

縦横線が抑えられていることがわかりました。

入力画像に対して周囲を暗くする方法としてはCIFilterのCIGaussianGradientが使えます。

CIImageExtension.swift

    // 画像の中央だけ円状に残して黒く塗りつぶす
    func gaussianGradient(radius: Float) -> CIImage? {
        // 中央だけ白いフィルターを生成
        let radialMask = CIFilter.gaussianGradient()
        radialMask.center = .init(x: extent.origin.x + (extent.width / 2),
                                  y: extent.origin.y + (extent.height / 2))
        radialMask.radius = radius
        radialMask.color0 = CIColor(cgColor: UIColor.white.cgColor)
        radialMask.color1 = CIColor(cgColor: UIColor.black.cgColor)
        
        // 自身と上記のフィルタ生成画像を乗算
        let multiply = CIFilter.multiplyCompositing()
        multiply.inputImage = self
        multiply.backgroundImage = radialMask.outputImage
        
        return multiply.outputImage
    }

CIFilter.gaussianGradient()でCIGaussianGradientをインスタンス化し、画像中央からradiusの大きさだけマスクを作成します。これを入力画像に乗算(CIMultiplyCompositing)することで、画像の真ん中だけ残る画像が作れます。
ちなみに、CIFilterが用意している同様のフィルタCIRadialGradientでも試したところ、CIGaussianGradientよりも縦横成分が残るような結果になりました。

2.カラー画像を8bitグレースケールに変換

FFTで処理したいのは１チャネルのグレースケール画像です。
ここで、キャプチャ画像を次のようにCIFilterでグレースケール化してそれを単純にCGImageに変換しただけだとRGBのチャネルとアルファチャネルが残ったままとなります。

let grayscaleFilter = CIFilter.maximumComponent()
grayscaleFilter.inputImage = gradientImage
guard let grayscaleImage = grayscaleFilter.outputImage else { return }
        
let baseImage = ciContext.createCGImage(grayscaleImage, from: grayscaleImage.extent)

そこで生成するピクセルフォーマットを次のように指定することでカラー画像から一気に１チャネルの画像が生成できます。

// カラー画像をグレースケールに変換する
guard let baseImage = ciContext.createCGImage(gradientImage,
                                            from: gradientImage.extent,
                                            format: .L8,
                                            colorSpace: CGColorSpaceCreateDeviceGray()) else { return }

ただ、カラー→グレースケールの変換式は不明です。.L8の説明が「An 8-bit-per-pixel, fixed-point pixel format in which the sole component is luminance.」なので、輝度であることは確かだと思いますが、計算式はわかりませんでした。
もう一点、後々、FFTで処理する前にFloatに変換するので、ここで.L8の代わりに.Lfを指定してFloat型の配列にしてしまうこともできます。ただ、ここでFloatにするとFFTの変換結果でなぜかノイズが強く出て見づらいのでUInt8に変換しておきます。

3.二次元フーリエ変換

２次元高速フーリエ変換のコードは、Appleの記事『Perform Fourier Transform on 2D Real Data』と同様です。UInt8配列をFloat配列に変換したり、毎回vDSP_create_fftsetup()を行わないようにしたりする程度の変更をしています。

// UInt8配列をFloat配列に変換
var pixelFloat = [Float](repeating: 0, count: Const.imageDataSize)
vDSP.convertElements(of: pixelUInt8, to: &pixelFloat)

// 略

pixelFloat.withUnsafeBytes { imageDataPointer in
    _ = [Float](unsafeUninitializedCapacity: complexElementCount) { realBuffer, _ in
        _ = [Float](unsafeUninitializedCapacity: complexElementCount) { imagBuffer, _ in
            // 画素値配列を複素数配列に変換
            var splitComplex = DSPSplitComplex(
                realp: realBuffer.baseAddress!,
                imagp: imagBuffer.baseAddress!)
            vDSP_ctoz([DSPComplex](imageDataPointer.bindMemory(to: DSPComplex.self)),
                      2,
                      &splitComplex,
                      1,
                      vDSP_Length(complexValuesWidth * complexValuesHeight))
            // ２次元フーリエ変換
            vDSP_fft2d_zrip(fftSetup,
                            &splitComplex,
                            1,
                            0,
                            countLog2n,
                            countLog2n,
                            FFTDirection(kFFTDirection_Forward))
            // （略）後述
        }
    }
}

4.複素数を画像として扱えるように加工

FFTの出力splitComplexは複素数なので、これを画素値として使えるように加工します。

// 複素数を振幅スペクトル（絶対値をとって振幅の大きさ）に変換
vDSP.absolute(splitComplex, result: &amplitudeSpectrum)
// 振幅スペクトルをパワースペクトル[dB]に変換
powerSpectrum = vDSP.amplitudeToDecibels(amplitudeSpectrum, zeroReference: 1)

まず、absolute()で絶対値をとって振幅の大きさに変換します。この絶対値は大きな値なので、これをそのまま画素値にしてしまうと振幅の振れ幅が大きすぎて画像が砂嵐のようになってしまいます。そこで次にのようにamplitudeToDecibels()でデシベルに変換します¹。

例

var amplitudeSpectrum: [Float] = [10,100,1000]  // 20, 40, 60 に変換される
var powerSpectrum = vDSP.amplitudeToDecibels(amplitudeSpectrum, zeroReference: 1)
// amplitudeToDecibels()はiOS13から利用できるメソッド。

ここでのFFTの出力は右中央に高周波成分、左上・左下に低周波成分が出力されるため、左中央が低周波成分になるように上下を入れ替えます。

powerSpectrum.withUnsafeMutableBufferPointer { pointer in
    let p1 = UnsafeMutablePointer(pointer.baseAddress!)
    let p2 = p1.advanced(by: complexElementCount / 2)
    vDSP_vswap(p1, 1,
               p2, 1,
               vDSP_Length(complexElementCount / 2))
}

vDSP_vswap()を使って配列の要素を入れ替えます。

5.（見た目のため）エイリアシング?!を生成

　本記事は調べながら書いているのですがAppleの記事をベースに実装すると、左右対称の周波数成分の出力（エイリアシング）が得られません。vDSP_fft2d_zrip()は入出力で同じサイズのメモリを使い、出力結果で使えるのは半分（実部と虚部で一組）なので、エイリアシングがあるとするならさらに半分が鏡映っぽく出力されるはずですが、そのような出力は得られません。
　物体の方向を知る上では不要ですが、パワースペクトル画像として見かけるあの見た目が欲しいので邪魔なはずのエイリアシングを作ります。本サンプルではpowerSpectrumの配列をreverse()して作成しています。

6.画素値配列をCGImageに変換

画素値の配列から画像(CGImage)を作ります。輝度情報だけなのでグレースケール画像を作ることになります。作り方は前回の記事『画素値の配列からCGImageを作る』で解説しているので興味のある方は参照ください。
一点、CGImage与える画素値を次のようにFloat型のからUInt8型に変換して得ています。

// Float配列をUInt8配列に変換
vDSP.convertElements(of: values,
                     to: &uIntPixels,
                     rounding: .towardZero)

参考URL

• Apple サンプルコード：Halftone Descreening with 2D Fast Fourier Transform
• Apple ドキュメント：Data Packing for Fourier Transforms
• タコさんブログ

1. この記事の目的は「方向を知る」なので単位をdBにせず「絶対値を500で割って使う」でも良いのかもしれません。 ↩