PIVをやってみる．part 3 FFT相互相関法で作成

はじめに

part 1で無料のOpenPIVを試してみたが内部の処理がわからずパラメーターの最適化ができなかった．part 2では直接相互相関法で作成したPIVのプログラムを試したが，上手く渦が計測できた．しかし，処理速度が遅いため，次はFFT相互相関法を試してみる．

FFT相互相関法

計測時間が前の画像を第1画像，後の画像を第2画像とする．画像を一定のサイズに区切って第1画像と第2画像に同サイズの探査領域を設定する．両方を二次元のフーリエ変換を行い，第2画像の方は複素共役をとる．第1画像をフーリエ変換したものと第2画像をフーリエ変換して複素共役をとったものの積を逆フーリエ変換すると相互相関係数分布となる．FFT(Fast Fourier Transform)を用いることができるのでプログラムは簡単にでき，GPUとの相性がよく高速化に適していると言われている．

プログラムの作成，結果

OpenPIVと比較するために前回と同様に画像はPIVChallengeのから1st PIV Challenge (Sept.14-15, 2001, Göttingen, Germany) のCase 3: medium particle density and small particle を使用した．Numpyの関数を用いて二次元のFFTをを行った．img_a_trimとimg_b_trimは同じサイズ，同じ位置で設定して，同じ粒子が探査領域内で計測できるように粒子の移動量を考えて設定する．前の時間の方はnp.conjで複素共役をとっている．最後に逆フーリエ変換した後に複素成分が残ってしまうため．np.realで実数部を使っている．

import numpy as np

img_a_trim_fft = np.fft.fft2(img_a_trim)
img_a_trim_fft = np.conj(img_a_trim_fft)
img_b_trim_fft = np.fft.fft2(img_b_trim)
img_ab_trim_fft = img_a_trim_fft * img_b_trim_fft
img_ab_trim_ifft = np.fft.ifft2(img_ab_trim_fft)
img_ab_trim_ifft = np.real(np.fft.fftshift(img_ab_trim_ifft))

直接相互相関法を用いたpart 2と同様に自作のプログラムでは渦が計測できており，うまくいっている気がする．探査領域を半分にしており4倍の計算負荷がかかるはずだが，一対の画像の計算に3秒程度で早い．2倍の高速化なので今後は並列化などでさらに高速化を目指したい．

終わりに

OpenPIVが上手くいかなかったため，PIVの更なる理解のために自分で直接相互相関法のプログラムを作成して，次にFFT相互相関法を試した．Numpyを用いて相互相関係数を求めてPIVをしたが，渦を計測できており，高速化もできた．次回はこれらのプログラムを用いてBOS法(Background Oriented Schlieren)と呼ばれるシュリーレン法を試してみたい．

参考

PIVハンドブック，可視化情報学会
numpy.fft.fft2 — NumPy v1.24 Manual