OpenPIV を用いて流体の運動を計測してみた．
前置き
OpenPIV とは
OpenPIV とは，オープンソースの PIV ライブラリで，Matlab, Python, C++ にて利用できる．今回は Python 版を使う．少なくとも Python のライブラリはまだ β 版である．
PIV とは
PIV とは， Particle Image Velocimetry の頭文字をとったもので，日本語では粒子画像流速測定法という．画像処理を用いて流体中に混入した粒子の移動をとらえることで，流体の運動を計測する手法である．
PIV については以下のサイトにわかりやすくまとめられているので，詳しくはそちらを参照してください．
PIV 事情
企業や研究機関で PIV を行う際は，専用の商用ソフトウェアを用いるのが一般的だと思われる．しかし，高額な PIV システムと一緒に購入する必要があったり，処理内容がブラックボックスだったりといったことが問題となる場合がある（ソフトウェアが PIV システムのコントローラを兼ねていたり，PIV 自体それほど複雑な処理を必要とはしないので，多くの場合問題にはならない）．
例えば，気軽に PIV を行いたいという場合や，解析ソフトウェアではかゆいところに手が届かない場合，処理内容を正確に理解する必要がある場合はソフトウェアを自作することになる．とは言っても 0 から作るのは大変なので，ソースコードが確認可能なオープンソースのプログラム（OpenPIV）を利用する．
OpenPIV の setup 手順
setup 手順はこちらを参考にしてください．私が試した当時より，
$ pip install cython numpy
$ pip install openpiv --pre
私が試した時点では**サンプルコードに誤り**があった．インストールを行った上でサンプルコードを実行したときにエラーになる場合はサンプルコードの誤りの可能性がある．
参考ページ：OpenPIVをMacにインストールする
マニュアル
開発中のため完全ではないが，Open PIV の使用方法はこちらの文書で確認できる．
既存ソフトウェアとの性能比較
以降に私が研究に使用している解析ソフトウェア（PIV システムの制御，データベース管理を兼ねたもの）で解析したものと OpenPIV で解析したものを比較した結果を示す．
比較方法
既存のソフトウェアは，画像の撮影から解析までの一連の流れを管理できるもの．既存のソフトウェアのデータベースから画像を出力することで解析結果の比較を行った．
Interrogation Area のサイズやオーバーラップ幅などの条件は揃えるようにした．
解析対象
下図のような流れ場を対象に 2 次元 PIV 計測を行った．その他条件や，流れ場が発達しているかどうかなどについてはここでは議論しない．
計測領域：
撮影された画像：
実際には 1 pixel あたりの輝度階調が $2^{12}$ 程度のデータを粒子の移動前・移動後の 2 時刻分使用して以下の解析を行った．
比較① 画像相関
流体中の粒子が撮影領域に流入・流出するため，撮影領域の端の部分は相関が取られにくく，正しく計測されにくい．この時点では OpenPIV を利用した結果の方が誤ベクトル（誤った相関が取られたと考えられる部分）が多い．
比較② 誤ベクトル除去
画像相関のみ行った時点での誤ベクトルをほとんど除去することができた．特に空間平均値は既存の商用ソフトウェアと同等の解析結果を得ることができた．
パラメータを自分で入力することになるし，頑張ればソースコードを覗くこともできるので，画像相関や誤ベクトル除去の処理について理解を深めるきっかけになるのではないかと思う．
使用したソースコード
下記では解析対象の画像データが bmp ファイルになっているが，実際には 1 pixel あたりの輝度階調が $2^{12}$ 程度のデータを使用した．frame_a, frame_b がそれぞれ numpy ndarray 形式の 2 次元配列になっていれば処理を実行可能．
開発中のため完全ではないが，Open PIV の使用方法はこちらで確認できる．
import matplotlib.pyplot as plt
import numpy as np
import openpiv.filters
import openpiv.process
import openpiv.scaling
import openpiv.tools
import openpiv.validation
# read image data
frame_a = openpiv.tools.imread( 'pic01.bmp' ) # before
frame_b = openpiv.tools.imread( 'pic02.bmp' ) # after
# correlation
u, v, sig2noise = openpiv.process.extended_search_area_piv(
frame_a.astype(np.int32),
frame_b.astype(np.int32),
window_size=32, #Interrogation Area サイズの 1 辺の長さ [pix]
overlap=16, #Interrogation Area のオーバーラップ幅 [pix]
dt=0.04, #2 フレーム間の時間 [s]
search_area_size=128, #相関をとる範囲のサイズの 1 辺の長さ [pix]
sig2noise_method='peak2peak' #後に誤ベクトル除去に使う S/N 比の取得方法
)
# invalid vector validation
u, v, mask = openpiv.validation.global_val( #最大値と最小値を使った誤ベクトル除去
u,
v,
(-200, 500), #u の最大値，最小値
(-100, 100) #v の最大値，最小値
)
u, v, mask = openpiv.validation.sig2noise_val( #相関の強度を使用した誤ベクトル除去
u,
v,
sig2noise, #S/N 比
threshold = 1.3 #しきい値
)
u, v = openpiv.filters.replace_outliers( #除去されたベクトルの補間
u,
v,
method='localmean', #周囲のベクトルの平均値を使用して保管する
max_iter=5, #反復回数
kernel_size=2 #カーネルサイズ
)
# set [pixcel -> meter] factor
x, y = openpiv.process.get_coordinates(
image_size=frame_a.shape,
window_size=32,
overlap=16
)
x, y, u, v = openpiv.scaling.uniform(
x,
y,
u,
v,
scaling_factor = 22400 # 1 m あたり何 pixel か
)
# plot data
plt.imshow(u,cmap="jet")
plt.title(r'$U$ (OpenPIV)')
plt.xlabel(r'$x$ [grid]')
plt.ylabel(r'$y$ [grid]')
plt.clim(0,0.01)
plt.colorbar()
fig.tight_layout()
plt.show()
# save
openpiv.tools.save(x,y,u,v, 'expt_001.txt')