高速度ビデオでの運動解析プログラムの開発- OpenMP 等を用いた相関計算並列化によるトレース高速化 横山直樹（東海大学総合科学技術研究所）

Development of versatile motion analysis program based on normalized correlation for
high-speed videography - Parallel processing with OpenMP and MPI Naoki YOKOYAMA (Institute of Research and Development, Tokai University)
キーワード：高速度ビデオシステム、運動解析、画像相関
Keywords: High-speed video system, Motion analysis, Image correlation

Motion analysis program for images acquired by high-speed videography was developed andtested. This program was originally developed for Shimadzu’ ultra high-speed video recording system. Thenit was extended to handle general movie AVI file. Fully automatic analysis method of position detection
based on cross-correlation was proposed. Experiments were performed for rotation of cooling fan and butanecombustion. Laser interferometer was utilized to visualize refractive index change caused by butane-aircombustion. Temporal resolutions were 1,000 frames/s to 100,000 frames/s. Cross-correlation calculations
can be very heavy load to computer systems, so any kind of parallel processing should be used if available.
CPU based parallel processing method OpenMP and network based method MPI were discussed this time.
1. はじめに高速度ビデオシステムによって記録された動画を対象とし、物体の運動を自動解析するプログラムを開発し、その有効性を実験的に検証した。特に時間分解能を高めた高速度ビデオシステムにおいては、空間分解能に制限があり、各フレームにおける対象物体の位置の計測精度が低くなりがちであるが、画像相関値を評価関数とし、その極大値が得られる位置を物体位置として把握するアプ
ローチでは、位置情報をサブピクセルの単位で求めることが可能である。
2. VMAP(Versatile Motion Analysis Program for high-speed videography)
このプログラムの主な特徴は以下のとおりである、
 島津製作所の 100 万駒/秒超高速度ビデオカメラで記録されたファイルに関しては、その内部構造までを把握しており、ヘッダ部分から撮影速度等の記録時の主要なパラメータを読み取る。
 上記システム以外で得られた動画であっても、Windows の標準的な AVI (Audio VideoInterface) ファイルであれば解析可能である。この場合は撮影速度等を手動で正しくセットする必要がある。
 物体の位置把握は、相互相関計算により行われる。位置検出は、サブピクセルの精度で実行される。
 変位後の物体位置をさがすパターンは、螺旋である。これをスパイラルサーチと呼ぶ。図示すると Fig.1 のようになる。時間分解能が充分な場合は、効率的なサーチが可能である。

12



テンプレートを回転させながら相関計算を行うことで、物体の動きに回転成分が含まれている場合にも、その角度を検出することができる。
静止画を扱うことができる。時間的に連続する静止画は、それらを一旦 AVI ファイルに連結してから取り扱うことができる。
Fig. 1 Spiral search
運動解析結果の軌跡を、元の動画に重ねた形で動画として出力できる。
マニュアル解析が可能である。
各フレームの画像を、ヒストグラムの等値化アルゴリズムにより強調できる。
ROI の数が多いか、またはその領域が大きいときは、相関計算にかなり時間がかかることが見込まれるが、前者の場合はネットワーク上で複数の計算機を用意することで、MPI (Message Passing Interface) を用いた分散処理が可能である。
相関計算を工夫することで、上記のネットワーク分散処理に加えて単一ノードでの処理の高速化も図ることができる。この場合は CPU ないし GPU により並列処理を行う、これら二つの併用も可能である。

Fig.2 にこのプログラムの操作パネルを示す。

Fig. 2 Operational Panel of VMAP
物体の着目領域である ROI(Region of Interest)は、フレーム毎に相互相関値の極大点を探索することで追跡できる。この場合の正規化相互相関(ZNCC:Zero-mean Normalized Cross-Correlation)は、次の式で与えられる。

13 ( fi  f )  ( g j  g ) ( fi  f )   ( g j  g )22

ここで fi と gj はそれぞれ ROI と探索対象画像上にとったその対応部分上の画素を示す。 f と gは対象領域にわたって計算された画素の輝度の平均値である。Fig.2 の左端中央に示されている
Graph ボタンをクリックするだけで、自動解析された結果は、外部のプログラムを使用してグラフ表示することができる。例を Fig.3 に示す。

260240220Position

200180160
PositionX
PositionY

1401201000.00.5

1.01.5Time

2.02.53.0x10-3

Fig. 3 Typical result by VMAP motion analysis (plot by Igor)
3．実験および議論
3.1 ブタンガスの燃焼
レーザ干渉計を用い、ブタンガスの燃焼による光路長変化を可視化し、十分な時間分解能で記録した。撮影速度は、100,000 fs/s までを用いた。動画から抜き出した静止画の例を Fig.4 に示す。

Fig. 4 Still Images from butane combustion by laser interferometer 5,000 fs/s

Fig.4 に見られるような干渉縞の動きを自動解析するためには、上述の正規化相関はそのままで14は適用できない。サーチ方向を限定することで、部分的に自動解析が可能となる。今回はマニュア
ルで縞の拡がりを解析した。Fig.5 に時間に対する最外郭の縞の空間的拡がりの変化を示す。縦軸は画素単位である。水平方向と垂直方向の差は、電極の存在によるものと考えられる。

Expansion of outmost fringes

3530
Horizontal
Vertical25
20
151050.0

0.5

1.0

1.5

2.0
Time [s]2.5

3.0

3.5x10-3

Fig. 5 Experimental result of Butane-air combustion
3.2 冷却ファンの運動解析
つぎに、回転運動も伴う場合の解析が可能かどうかを検証するために、PC 用の高速冷却ファンを対象として実験を行った。ファンの設定回転速度は最高毎秒 3,800 r.p.m.であった。結果の例はFig.3 に示したものである。単なる並進運動でなく、対象が回転する場合は、相関計算のためのテンプレートの回転も考慮に入れる必要がある。そうすることにより、角度が大きくなった時点でも、適切な位置の検出ができる。Fig.6 にその例を示す。左側(a)が回転補正あり、右側(b)が回転補正なし（右端付近でマッチ失敗）である。

Fig. 6 Effects of rotational compensation (a)without compensation(b)

回転補正はテンプレートを微少角回転させたものをアフィン変換によって求め、本来のテンプレ

15

ートでなくこの変換後のものと対象画像上の領域との相関値を計算することで行う。この性質を利用することで、ROI の位置だけでなく回転角度の計測も可能である。Fig. 7 には角度も含めた解析結果の例を示す。時間-角度のプロットに直線を当てはめ、その勾配を計算することでファンの回転速度が得られる。結果は 3,830 r.p.m.であった。こうした解析は、データをエクスポートしたプロットプログラム側で、そのプログラムのマクロ等を用いて実行することが可能である。

250

Position

200
PositionXPositionY
Angle

1501005000.00.5

1.0

1.5

2.0

2.5

3.0x10-3

Time
Fig. 7 Analysis result for rotation of cooling fan @3,800 r.p.m.
3.3 複数の気泡の同時トラッキング
正しく撮影速度を設定することにより、この運動解析プログラムを、一般の AVI ファイルにも適用可能である。Fig.8 に示すのは、島津製作所の超高速度ビデオシステムより高い空間分解能を持った、Photron 社製 FASTCAM-Ultima(512x512)で記録した動画ファイルを解析した結果である。
Fig.8 の左の 8 個のターゲットの気泡の位置をフレーム毎に示したのが右の結果である。前述のように解析結果を動画として書き出すことも可能である。

Fig. 8 Tracking of multiple objects (bubbles in water)

16

Fig.8 のような多数の気泡を対象とする場合、相関計算の負荷が大きくなることが予想されるので、数百個の ROI が存在する場合も考慮し、なんらかの負荷分散が必要となる。
3.4 負荷分散のストラテジ
対象部位が多数存在するような場合、負荷分散の方法は原則的に下記のものが考えられる。

１） フレームによる分割（時間軸での分割）
これは実際には、前のフレームにおける物体位置が後のフレームのスタート位置として影響するので、順次位置を求めることが必要なので、かなり困難である。
２） 位置による分割（空間軸での分割）
ROI によって分割する場合は実現容易であるが、ROI の個数が少なく、さらにそのサイズが大きい場合に負荷分散の効果があまり大きくならないことが予想される。つまり ROI の数以上には分割できないということである。その場合でも、相関計算ルーチン自体を並列化することで、マルチコアを利用した単一ノードでの高速化は可能である。
空間軸での分割を用いて、ネットワーク分散処理により負荷軽減ができることは後で MPIに関連して述べるが、ここではまず単一ノードでの相関計算の高速化をテンプレートの形状の工夫と並列計算により目指す。
3.5

ROI のサイズおよび形状による相関計算の高速化
Pixel count v.s. Proc. time
3.5

3Proc. time in sec

2.5
Athlon 1.2G
Athlon 1.2G 384MB
X6800@2.93GHZ
Athlon 5200+
Pentium3 3.2G
Corei7 940@3666
Corei7 965@3342
Athlon x2 6000+

2

1.5

10.5

0050

100

150

200

250300350

AREA size in pixel

Fig.9 Relation between pixel count and processing time
Fig.9 に各 ROI の画素単位での大きさと処理時間の関係を、最近のものを含めた種々のマシンで計測した結果を示す。対象画像と 8 個の ROI は Fig.8 と同様である。画素の数と処理時間はほぼリニアな関係で、CPU の速度やメモリの読み出し速度によりその傾斜が異なるグラフとなっている。直17線の Y 切片がゼロになっていないのは、ディスク I/O 等の時間のためと思われるが、それもまたマシンによって異なる結果となっている。
次に ROI（テンプレート）と対象部位との過不足性とトレース結果の関係について述べる。Fig.10
はテンプレートが余剰な場合、ちょうど外接する場合、内接する場合について、計算時間とトレース結果を示したものである。計算時間は左から 3.35 秒、1.56 秒、0.49 秒で、トレース結果は左二つの場合はほぼ同一で正しい結果が得られている。

3.35

1.560.49350

300

250200150100500380390400410420430440

Fig.10 ROI size’s effects to traceability

結局矩形のテンプレートを用いる限り、対象部位にちょうど外接するサイズが最適なテンプレートということがわかった。しかし内部の均一な部分の計算は省ける可能性がある。そこで主に対象部位の輪郭領域に物体情報があると考え、輪郭部のみを相関計算に使うアプローチを試してみた。結果の一部を Fig.11 に示す。本来の 195 画素のうちエッジ情報を持つ 39 画素に対してのみ計算することで、計算時間が 0.36 秒から 0.18 秒と半減し、かつ相関値そのものも大18きく、より良好なトレース結果が得られた。

Fig. 11 Traceability results by weighted templates
3.6 OpenMP による並列処理での相関計算高速化
スレッド数と処理時間
7060

処理時間

50
core2 duo T7800 2.6GHz
Core 2 Extreme X6800 2.93GHz
P4 3.2GHz
Core i7 965
Linux icc
Linux gcc
Core 2 duo 2.2GHz
P4 3.2GHz

40302010002

468101214

スレッド数

Fig.12 Processing time v.s. numbers of threads (OpenMP)

19

Fig.12 には 32 個の ROI を対象とし、複数の相関計算をスレッドに分配することで高速化を図った
場合での、処理時間とスレッド数との関係をプロットしたものである。もちろん CPU やコンパイラの種類にも依存するが、その CPU のネイティブなコア数やスレッド数の場合に処理時間が最小化されることがわかった。たとえば 2009 年時点で最新の Intel の Core i7 の場合は、ハイパースレッドと呼ばれる機構により 1 コアあたり 2 個のスレッドが使えるので、スレッド数８で処理時間が最小となっている。
3.7

MPI による相関計算高速化

次に複数のホストが高速ネットワークでつながれている環境で、上述の OpenMP と同様な並列化の戦略で相関計算を高速化することを目指した。Fig.13 は横軸にランク数（ノード数にコア数をかけた同時に実行される処理の数で、この場合は１から６４にわたる）を、縦軸に秒単位の処理時間を取ったもので、筑波大学、東京大学、京都大学の三大学が共同で開発した、PC クラスター型の分散並列型スーパーコンピュータである T2K システム上で計測した。ゲストアカウントで処理を試行した関係で最大のランク数が 64 であった。本来の同時利用可能な最大単位は 512 ノード(8,192コア)である。

MPI effects on T2K80
node 170

node 2
node 360
Processing time [sec]node 450403020100010203040506070

Total numbers of Rank

Fig. 13 Total number of Rank v.s. Processing time at T2K guest account 256 ROIs

処理時間はランク数に逆比例して減少しており、並列処理のためのオーバーヘッドがほとんどない20理想的な傾向をしている。その理由は負荷(ROI 256 個)が十分大きいためであると考えられる。
5. 結論
テンプレートのエッジ情報に着目し、エッジ成分が高い画素のみを相関計算に使用することで、計算の高速化だけでなくより良い追跡が可能であることがわかった。また OpenMP を用いることにより、複数の相関計算を並行して走らせることができ、それによってコア数程度の高速化が可能となることがわかった。ただしエッジ情報抽出の際のしきい値などを自動的に決定するまでにはいたっていない。従って、多数の気泡ないし粒子を対象とし、重み付きテンプレートと相関計算を併用して運動解析を行う場合には、自動的にテンプレートを生成することが必要であり、今後の課題である。

参考文献[1] N. Yokoyama: Position detection by cross-correlation with linear scan camera(one-dimensional high-speed video camera), Int. J. Japan Soc. Prec. Eng., Vol.
33, No.1 (Mar. 1999) 57-58[2] 横山直樹,高速度撮影とフォトニクスに関する総合シンポジウム 2007 講演論文集,2007, pp.37-40 (2007)[3] MPI http://www.mpi-forum.org/[4] OpenMP http://openmp.org/

21