Python/OpenCVでPIV計測!粒子移動をベクトル化する

  • このエントリーをはてなブックマークに追加

撮影した映像から粒子の流れや物体の移動量を可視化するPIVに興味を持ちました。市販のPIV計測ソフトは高価なので、ここではPython/OpenCVを使って我流のPIV計測コードを書いてみた方法を紹介します。

こんにちは。wat(@watlablog)です。ここではPython/OpenCVを使って粒子や物体の移動をベクトル化するPIVコードを書いてみます

本記事の目的

無料でそれっぽい可視化をする

PIV(Particle Image Velocimetry)とは、粒子画像流速測定法を意味する計測技術の事です。おそらくメーカ勤務の実験屋さんとかはよく耳にする技術だと思います。

PIV計測は、一般に高速度カメラを使ってある程度速い現象を可視化する時に使われます。そのため通常はPIV計測システムとして高速度カメラ、レンズ、煙発生装置、トレーサ粒子、レーザーシート発生装置、暗幕、PC…と一緒におまけのようにPIV解析ソフトを購入します。

するとあら不思議、1千万円を超える金額が必要になってしまいます。

ソフトはシステムに最適化されており、様々な機能がついていると思います。もちろん品質の高い本格的な計測のためには必要な内容なのですが、ライトな用途で使いたい時にはソフトの習得にかかる工数が腰を重くしてしまいます。

ちょっと手持ちの動画から動きを可視化したり変位量を出したいだけだから高いソフトは必要無いかな。

本記事ではそんなライトな用途、何ならスマホやWebカメラで撮影した動画から無料でPIV解析を行う事を目的とします。

目標

Advertisements

オープンソースで完結させる

無料でPIV計測コードを書くため、当然無料で利用可能なプログラミングのみで実現させます。今回はPython3.7とOpenCV(opencv-python 4.2.0.34)を使っています。

Python

この記事をご覧の方はPythonについては十分な知識を持っていると想定します。ただ、もしこれからPythonを始める、まだインストールしていないという方は以下の記事を参考にして下さい。

僕は宗教上の理由(?)からAnacondaを使っていない事に注意が必要ですが、Windowsであれば「Pythonインストール方法とAnacondaを使わない3つの理由」、Macであれば「macOSにPython3をインストールする方法をまとめてみた」にインストール方法をまとめています。

Pythonは他にももっと便利な環境構築方法がありますが、僕はこちらが一番初心者向けと思いました(人によります)。

OpenCV

OpenCVは画像処理をしている人なら誰もが知っていると思われる有名なライブラリです。C++ベースの高速な処理が可能で、Pythonバインディングで使うOpenCVはNumpy形式のデータがそのまま使えます。

この辺の事は「Pythonで画像処理!OpenCVのメリットとインストール法」に記載しましたが、当ブログの画像処理カテゴリはほぼOpenCVで構成されています。興味のある方は是非ご覧になって下さい。

参考)OpenPIVというライブラリもある

どうやら世の中にはOpenPIVなるライブラリもあるようです。ただ2021年1月現在はベータ版のみしか無いみたいです。

ちょっと使ってみましたが、情報が少なく、パラメータがブラックボックスなので、僕にとっては高価なPIVソフトを使うのと似たような感覚でした。今後開発が成熟してユーザが増え、情報が沢山出てきたら使ってみようと思います。

今後に期待しているため、以下に公式ページをメモしています。

公式ページ:OpenPIV

まずは理想状態の粒子流れで検証できるように作る

コードは検証しながら作ります。しかしながら当ブログは趣味で書いているため、サンプルの動画がありませんでした(自宅にPIV計測システムがあれば良かったのですが…無念)。

そこで事前に「Python/PIVの検証用に管内の粒子流れ動画を作ってみた」という記事でPIVで計測するようなサンプル動画をPythonで作る所から始めました。

以下の動画がその例です。この動画はあえて速度分布のイメージを\(x=0\)にセットしていますが、上記記事のコードを使えばこの速度分布の形や最大値位置、流れの速さといったパラメータを変更する事が可能です。

速度分布の異なる粒子流れ動画

こんな感じ。

最大流速位置を変更した粒子流れ動画

まずは粒子画像を使ってコードが適切に機能するかを確認する目標を立てました。

粒子動画以外で機能するか確認する

今回のPIV解析コードは画像相関法(相関係数を使ったテンプレートマッチング手法)を使って構築します。

そのため粒子画像で無くてもテンプレートマッチングが使えそうな画像ならば機能するはずです。

本記事では基本的なコード検証を粒子画像を使って行い、応用編として最後に他の画像への適用を確認する事も目標とします。

それでは早速以下からコードの説明と結果を紹介します。

Python/OpenCVでPIV解析をするコード

事前準備

ファイル構成と用意する画像

本記事のコードはメイン.pyファイル実行フォルダ下に「input」フォルダを作成し、その中にPIV解析用の動画を画像形式で、さらに連番のファイル名を付けて置いておく事が必要です。

ファイル構成

Python/PIVの検証用に管内の粒子流れ動画を作ってみた」で作成したサンプルファイルならば既にこの形式でデータが出力されるので、そのまま使用する事が可能です。

その他の動画を使う場合で、動画を画像ファイル化する方法は「Python/OpenCVで動画から静止画を間引いて抽出する方法」に記載していますので是非参考にしてみて下さい。

全コード

色々説明する前に、全コードと結果を先に紹介します。

以下のコード内piv関数が自作のPIV解析コードとなります。そのままでも動くように各変数にはデフォルト値を設定しています。

このコードを実行するとpiv_outというフォルダが自動生成されて、PIV解析後の画像が蓄積されます。inputフォルダに2枚の画像しか入れていなければ単純に1回のPIV解析として使えますが、複数枚の画像があれば後で動画にする事も可能です。

threshold=30という所は誤ベクトル処理の閾値です。後ほど説明しますが、流速の速い映像にはこの値を大きくする必要があります。

分布流れの可視化結果

中心流速が最も速く、端の方が遅いという代表的な管内流れを模した映像で試した結果を以下に示します。

cmap='jet'で色付けしているため、赤いほど速い流れ、青いほど遅い流れを示しており、狙い通り中心の方が速いという結果が得られました。

ハーゲンポアズイユ流れのようなもののPIV可視化結果

各高さ方向の平均流速をサンプル動画を作成した時の理論流速分布に重ねてみた図を以下に示します。

このように今回のデフォルト設定(overlap=0)では分布の位置がずれている事がわかりますが、overlap=12とする事で計測点が増え、さらに山の位置が改善される結果を得ました。

これはパターンマッチング時のウィンドウ走査の分解能を変更している事になりますが、overlapを増やすと計算数もどんどん増えていくので、処理が重くなってしまう事に注意が必要です。

デフォルト流れの精度検証結果

オフセット流れの可視化結果

今度は上の方の流速が速い流れの可視化をしてみた結果を以下に示します。
こちらも狙い通り。

オフセット流れの可視化結果
オフセット流れの精度検証結果

一様流れの可視化結果

管内に流速分布が無く、一様に流速が発生している場合の可視化結果を以下に示します。
画像中全域で同じ色のベクトルを得ました。ただ、どうやら流入付近は計算が苦手のようです。

一様流れの可視化結果

一様流れの場合の速度分布比較は以下。このような流れの場合は特にoverlap設定をしていてもしていなくても、精度的にはそんなに変わらないようですね(定量的で無くて申し訳ないですが…)。

PIV解析後の画像をGIF動画化するコード

先ほどの全コードではPIV解析部分のみにフォーカスしたかったため、GIF動画作成部分は別.pyにしていました。

以下のコードでpiv_outフォルダにある複数連番画像をpiv-movie.gifと動画にする事が可能です。参考までに。

コード解説

ストライド

画像内を走査していく処理で、ウィンドウをずらしていく幅をストライドと呼ぶそうです(CNNを学習した時に初めて知りました…スライドじゃないのか…)。

コード内では以下の文でストライド数を計算しており、オーバーラップに対応させています。

図解すると以下です。今回ウィンドウ幅\(w_{w}, h_{w}\)は同じですが、まずは幅方向にoverlapを考慮しながら走査をしていきます。

幅方向のストライド

幅方向終端まで到達した後は縦方向にウィンドウをずらしますが、その際にも同一のoverlapを考慮してストライドさせます。

高さ方向のストライド

この1画像に対するストライド走査がforループのjとkの部分に相当します。

画像走査処理をアニメーションで示したのが以下。オーバーラップ0だと隙間なくタイルが敷き詰められていくようなものです。

オーバーラップ0の場合の画像走査動画

横方向のみですが、オーバーラップを考慮すると以下のようになります。ウィンドウが重なっている様子が確認できると思います。

オーバーラップを考慮した画像走査

テンプレートマッチング

今回のPIVコードのメインはテンプレートマッチングで、以下の部分になります。

手法自体は相関係数を用いるNCC(Normalized Cross Correlation)を使っており、式は以下です。

\[ NCC(d_{x},d_{y})=\frac{\sum_{x=0}^{W-1}\sum_{y=0}^{H-1}\left [ I(x+d_{x},y+d_{y})T(x,y) \right ]}{\sqrt{\sum_{x=0}^{W-1}\sum_{y=0}^{H-1}\left [ I(x+d_{x},y+d_{y})\right ]^{2}}\sqrt{\sum_{x=0}^{W-1}\sum_{y=0}^{H-1}\left [ T(x,y)\right ]^{2}}} \]

しかし、この部分はOpenCVに任せています。詳しくは「OpenCVのテンプレートマッチングで変形量を算出する方法」の記事をご覧下さい。

やっている事は流れの前後で、ウィンドウで抽出した部分がどこからどこに行っているかを相関係数を使って計算しているのみです。

下図は流れの前後でウィンドウがしっかりと追従している(流れを追えている)例を示した図です。これを画像走査で画像全域で行います。

テンプレートマッチングの例

誤ベクトル除去

テンプレートマッチングに失敗すると、算出したベクトルが物理的にあり得ない長さになる場合があります。

通常は相関係数がいくつ以上のデータしか使わない…といった処理をするようですが、本コードではthreshold以上のベクトルを0とする事で対応しました(相関係数的には高くてもテンプレートマッチングが失敗しているものもあったので)。

誤ベクトル除去の有無を比較した図が以下です。今回はベクトルの色をベクトルの振幅最大値で自動的に付けているので、除去しないと色の分布も偏ってしまいます。

ベクトルプロット

最終的なベクトルのプロットは以下のコードのようにmatplotlibのarrowで行っています。vsfはベクトルに対するスケールファクターです。このarrowの中の変数をいじると、矢印の見え方が変わってくるため、画像サイズ等に合わせて使い分けてみて下さい。

ベクトル長さが0の時はプロットしないので、上で書いた誤ベクトル除去されたデータは反映されません。

主要なコードはこんな感じです。全体的に全コード内にコメントを残しましたので、詳細はそちらを参照して下さい。

応用:粒子動画以外への適用例

以下のGIF動画は「Python/OpenCVで動画速度(FPS)を変えて再保存!」とかで使っていた動画に対して、本記事で作成した自作PIVコードを適用してみた例です。

粒子動画以外でも、テンプレートマッチングさえうまくいけば動画内の物体運動量を測定する用途にも使えそうです。

粒子動画以外へのPIV適用例

実際にPIVコードを使う場合は1ピクセル当たりの長さ[mm]を使ってベクトルの校正をする必要があると思います。

今回はその部分まで作り込みませんでしたが、校正をする際は画像サイズの変化、カメラの角度、計算に使用したスケールファクターといった所に注意が必要です。

まとめ

本記事ではOpenCVを使ってPIVコードを自作してみました。

もちろん本格的なPIV計測で使用する市販のソフトと比べると、プログラムの効率も悪いと思うので計算速度も遅いと思います。

しかしこんなに簡単なコードでも誤ベクトル除去やベクトルプロットができ、それなりに流れが可視化できる、さらに粒子動画以外でも適用が簡単にできるというのは驚きました。

ライトな用途であればPythonで十分かなと思います。

今後は本記事のコードを応用させて物体運動の計測に特化したものも作ってみようと思います。

簡易PIVコードがかけました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメント

  1. 萬矢敦啓 より:

    突然のコントで失礼します.御社のサイトを頻繁に閲覧させて頂き,勉強させて頂いているものです.
    特にこのページのPIVが大変参考になりまして,私自信が開発しているプログラム(以後プログラムA)でも一部使用させて頂いております.その開発プログラムを何らかの形で第三者に配布をしたいと考えているのですが,著作権に関して確認したことがあり,本コメントを書いております.

    上記のような条件において,御社に対して何かロイヤリティが発生しますでしょうか?またプログラムAの関連するdefで御社のサイトに記載されているcopyrightを記述することで,著作権に対して,その権利を尊重していると考えることができますでしょうか?ちなみに私自信は無償で配布する予定です.

    著作権に関して勉強中で,非常に稚拙な質問で恐縮です.ご教授頂ければ幸甚です.よろしくお願いします.

    1. wat より:

      いつもご訪問ありがとうございます。
      当ブログのコードは全て非営利目的、かつ趣味で書いています。
      そのため、ご自由にコピーしてご活用や配布をして頂いて構いません。
      特にコードにURLやcopyrightを記載する必要もありません。

      但し、当ブログのコードはバグの保証や商用利用の際の不都合当は配慮していません(今回は無償とのことですが念のため)。
      使用している外部ライブラリによっては配布に関する条項がある可能性があるため、OpenCV等の公式ページ、OSS(Open Source Software)の考え方をご確認されることをおすすめ致します。

      以上、今後とも記事をご活用頂けると幸いです。

      1. 萬矢敦啓 より:

        ご丁寧な回答ありがとうございました.よく理解することができました.

  2. RM より:

    突然のコメント失礼いたします.私は現在理系大学で研究をしている学生です.いつも貴方のサイトを拝見し,様々なことを勉強させていただいています.
    現在,PIVの計測をパイソンを用いて行おうと考えているのですが,ベクトルマップの数値データをCSVファイルに書き出す手順で苦戦しております.このページを参考にベクトルマップを描写することはできましたがどうしても数値データの抽出の方法がわかりません.もしよろしければ,ベクトルマップを数値データとして抽出する方法をご教授願えないでしょうか.
    プログラミングについてまだまだ初心者であるため,初歩的な質問となってしまうことをお詫び申し上げます.何卒よろしくお願いいたします.

    1. wat より:

      ご訪問ありがとうございます。
      また、いつも参考にして頂き恐縮です。

      ベクトルマップを数値データとしてcsv化したいとのことですが、誤ベクトル除去も含めた最終的なベクトル成分毎の数値自体は
      「# 長さ0以外の場合に図にベクトル(dx, dyにそれぞれスケールを乗じた後のベクトル)を描画」部分に記載のようにcoordinatesに集約されています。

      例えばax1.arrow()内のx, y, dx, dyはそれぞれベクトルの始点と終点を意味していますので、vsfといった余計なスケールファクターをかけない状態の生データを記録していけば目的を達成できそうです。

      csvへの保存ですが、Pandasの.to_csv()を使うのが個人的には簡単だと思います。
      ↓参考
      https://watlab-blog.com/2021/04/17/csv-dft/

      ただ、保存フォーマットはcoodinatesをそのまま文字列や数値として記録するか、それとも成分毎にファイルを分けるのか…といった欲しい形に整形する必要があると思います。

      1. RM より:

        ご丁寧で分かりやすい回答,ありがとうございます.よく理解することができました.参考にさせていただきます.

コメントを残す

*