A14 Bionic GPU での配列の合計計算が高速であることを活かして、二値化された画像の重心をリアルタイムに計算する方法を紹介します。

＜完成イメージ＞
※赤丸が重心。二値画像の黒い部分が円（10円玉）の場合、円の中心が重心になっていることがわかります。
　

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

iOS14以降であれば実行できますが(二値化・クロージングの動作確認はできますが)、GPUでの重心計算確認にはiPhone12(Apple A14)以降が必要です。

■紹介する要素技術

CIFilterを使った二値化とクロージング
GPUを使った画像モーメントの計算

１．二値化とクロージング

今回の処理に必要とする画像は「一塊の物体（連結成分）」です。次のようなものです。

物体の重心を知りたいので処理がしやすい白黒画像
重心を計算したいので注目する物体が一塊になっている（穴がない）画像

前者を得る方法がグレースケール化・二値化で、後者を得る方法がクロージングです。

1-1）グレースケール化・二値化

グレースケール化

二値化をする前に画像をグレースケールに変換する必要があります。
CoreImageにはカラー画像をグレースケール化する次のようなCIFilterが用意されています。

① 白黒写真フィルムを模倣するフィルタ[CIPhotoEffectTona、CIPhotoEffectNoir、CIPhotoEffectMono]
② RGBチャネルの一番明るい画素を採用するフィルタ[CIMaximumComponent]
③ RGBチャネルの一番明暗い画素を採用するフィルタ[CIMinimumComponent]
④ RGBの各チャネルに対して重みをつけて色を変換できる（輝度にも変換できる）フィルタ[CIColorMatrix]

ここではアートっぽく加工されている必要性はないので①は使用しません。
明るい画像から暗い物体を把握したいので②を使います（撮影条件や二値化の閾値の決め方次第では③でも良いと思います。④は少しコードが増えるので見送ります）。

CIImageExtension

let grayscaleFilter = CIFilter.maximumComponent()
grayscaleFilter.inputImage = self
guard let grayscaleImage = grayscaleFilter.outputImage else { return nil }

ちなみに、わざわざカラー画像をグレースケールにするくらいなら、カメラの設定で最初からグレースケールを取得できないものか、と考えましたが次のコードで背面カメラがサポートするフォーマットを確認したところ420v、420f、x420のようにYCbCrっぽいものしかありませんでした(420vは動画用に輝度情報が狭いもの、420fはフルレンジの輝度、x420は不明だが'420'とあるのでカラー？）。

let videoDevice = AVCaptureDevice.DiscoverySession(deviceTypes: [.builtInWideAngleCamera],
                                                   mediaType: .video,
                                                   position: .back).devices.first
videoDevice?.formats.forEach { format in
    print(format.formatDescription)
}

二値化

グレースケール画像の明るい部分・暗い部分を白・黒の二値に変換します。
二値化の方法にはいくつかありますが、CIFilterとして二値化用に用意されているものは次の２つです（他に見当たりませんでした）。

① CIColorThreshold（手動で閾値を指定）
② CIColorThresholdOtsu（判別分析法により自動で閾値を決めてくれる）

②はヒストグラムから自動で白・黒の山の分離度が大きくなるように閾値を自動でを決めてくれるものです。今回は、撮影状態により物体だけをうまく抽出したい・微調整したいので①を使っています。
ちなみにCIColorThresholdOtsu は内部でMPSHistogramを利用しているようなのでCIColorThresholdより処理は重いと思われます。

CIImageExtension

let thresholdFilter = CIFilter.colorThreshold()
thresholdFilter.inputImage = grayscaleImage
thresholdFilter.threshold = threshold
guard let thresholdImage = thresholdFilter.outputImage else { return nil }

1-2）クロージング

二値化をしただけだと、対象の物体を一塊として見なすことができない場合があります。
次の二値画像は１０円玉とペンを撮影したもので、黒い領域に白い穴が目立ちます。この穴を塞いで一塊に近づける必要があります。

＜二値画像＞

膨張と収縮（モルフォロジー演算）

本記事では黒を物体として扱う前提としており、物体が黒く１つの塊としてみなせるようにする必要があります。　※従って白が背景・穴という前提です。

そこで、黒い部分を大きくする処理が膨張です。

＜黒い領域 の膨張処理結果＞

１回目	２回目	３回目

膨張は白い画素が黒い画素と接する場合、その画素を黒に書き換える処理です。この収縮処理により物体の穴など白い部分が小さくなったり消滅したりします。ただ、物体が大きくなってしまうことに注意が必要です。
次に、黒い部分（白い画素と接する部分）を白に書きかえる処理を収縮と呼びます。膨張後に収縮を行うと物体が元の大きさと同じ程度に戻ります。

＜膨張処理後の収縮処理結果＞

１回目	２回目	３回目

少しずつ黒い領域が小さくなっています。当然ながら収縮をしても消滅した穴は復元しません。
この、膨張→収縮の流れをクロージングと呼びます。

CIImageExtension

// 膨張・収縮処理（クロージング）
let closingDilateFilter = CIFilter.morphologyMinimum()
closingDilateFilter.inputImage = thresholdImage
closingDilateFilter.radius = Const.morphologyFilterRadius
guard let closingDilateImage = closingDilateFilter.outputImage else { return nil }

let closingErodeFilter = CIFilter.morphologyMaximum()
closingErodeFilter.inputImage = closingDilateImage
closingErodeFilter.radius = Const.morphologyFilterRadius
guard let closingErodeImage = closingErodeFilter.outputImage else { return nil }

黒い部分を膨張（白い部分を収縮）するフィルタがCIFilterのCIFilter.morphologyMinimum()で、黒い部分を縮小するフィルタがCIFilter.morphologyMaximum()です。

なお、プログラムの中で膨張・収縮のフィルタは１回ずつしか適用していませんが、CIFilterの内部(GPU処理)は各々３回のフィルタ（3x3カーネルっぽい）処理を実施しています。

＜GPUキャプチャ結果＞

この3x3フィルタの回数はCIFliterに設定しているradiusプロパティの大きさによって変わります。
（ちなみに、このGPUキャプチャ結果を見るとグレースケール化と二値化は１つの処理として結合されているようです）

補足：オープニング処理

クロージングは穴を塞ぐ(小さな白い領域を消す)処理でしたが、小さな黒い領域を消す方法としてオープニングがあります。処理としては収縮→膨張の流れです。画像上のポツポツとした小領域・ノイズを消すのに良いのですが、今回試した環境（背景色が一様）ではあまり効果がなかったので、実装上はコメントアウトしています。

２．モーメント特徴

2-1）計算方法

前述で求めた黒い画像領域の重心を画像のモーメントを使って求めます。
画像のモーメントMと、x軸上の重心、y軸上の重心は次の計算式で求めます¹。

モーメントM_{pq} = \sum{f(i, j)}\; i^p i^q \\
x軸上の重心 = \frac{M_{10}}{M_{00}} \\
y軸上の重心 = \frac{M_{01}}{M_{00}} \\

ここでf(i, j)は画素値が黒なら1、白なら0として計算します。
つまり、M00は黒画素の数、M10は黒画素のx座標値の合計、M01は黒画素のy座標値の合計 となり、これらを算出すれば重心が計算できます。
注目画素のみの値を使った単純な合計を計算するだけなので、A14 Bionic GPUで追加されたsimd_sum命令で計算しやすいと言えるかもしれません。

2-2）GPU側の処理

モーメントの算出はGPUとCPUで分担します。役割分担は次の通りです。

　①GPUは1,024画素を１つの合計単位として全画素分を計算
　②CPUは上記のGPUの結果を受け取り合計を計算

今回のサンプルは、448x448画素の動画を処理しています。1フレームあたりの画素数は 200,704 です。
①のGPUの計算は1,024個単位で合計するので、196個の合計値が得られます。
②のCPUの計算はその196個の合計を合計、ということになります。

以下は、GPU側のコードです。

shader.metal

kernel void group_max(texture2d<half, access::read> inTexture [[ texture(0) ]],
                      device MomentElement* output_array [[ buffer(1) ]],
                      uint2 position [[thread_position_in_grid]],
                      uint2 group_pos [[threadgroup_position_in_grid]],
                      uint simd_group_index [[simdgroup_index_in_threadgroup]],
                      uint thread_index [[thread_index_in_simdgroup]])
{
    // 1Thread Groupに32のSIMD groupがあるのでその各々の計算結果を格納するメモリを確保
    threadgroup MomentElement simd_sum_array[32];

    // 各スレッドのinput値は`0(黒)`画素を処理対象とする
    bool input = inTexture.read(position).r;
    int binary = !input;
    
    // simd group毎のモーメントを求める
    int i = position.x;
    int j = position.y;
    int m00 = simd_sum(binary);
    int m10 = simd_sum(binary * i);
    int m01 = simd_sum(binary * j);
    // 合計を一時保存
    simd_sum_array[simd_group_index].m00 = m00;
    simd_sum_array[simd_group_index].m10 = m10;
    simd_sum_array[simd_group_index].m01 = m01;

    // Thread Group内のすべてのスレッドの合計の計算＆一時保存を待つ
    threadgroup_barrier(mem_flags::mem_threadgroup);
    // thread group毎のモーメントを求める
    // 1つのSIMD Groupにで32スレッドあるので、１つのSIMD Groupのみで処理。
    if (simd_group_index == 0) {
        output_array[group_pos.y * threadGroupXsize + group_pos.x].m00 = simd_sum(simd_sum_array[thread_index].m00);
        output_array[group_pos.y * threadGroupXsize + group_pos.x].m10 = simd_sum(simd_sum_array[thread_index].m10);
        output_array[group_pos.y * threadGroupXsize + group_pos.x].m01 = simd_sum(simd_sum_array[thread_index].m01);
    }
}

スレッドグループ毎にm00、m10、m01を計算しています。
プログラムの流れは前回の記事『Apple A14 で追加されたGPUの命令で100万個の数値を合計』と同様なので、詳しくはそちらを参照ください。

2-3）CPU側の処理

CPU側はGPU側で途中まで計算された合計値の合計を行い、あとは計算式に従い重心を計算するだけです。

MomentCalculator.swift

var m00 = 0
var m10 = 0
var m01 = 0
for element in self.momentElements {
    m00 += Int(element.m00)
    m10 += Int(element.m10)
    m01 += Int(element.m01)
}
var centerGravityX = 0
var centerGravityY = 0
if m00 != 0 {
    // 重心
    centerGravityX = m10 / m00
    centerGravityY = m01 / m00
    // CIImageはy軸反転しているので元に戻す
    centerGravityY = imageSize - centerGravityY
}

一点、m00がゼロ（つまり画像に黒いところがない）場合だけ注意が必要です。

参考URL

wikipedia 画像のモーメント ↩