Без кейворда

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

本記事の目的

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

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

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

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

目標

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

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

OpenCV

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

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

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

まずは理想状態の粒子流れで検証できるように作る粒子動画以外で機能するか確認する

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

事前準備ファイル構成と用意する画像

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

全コード

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

import cv2 import glob import numpy as np from matplotlib import pyplot as plt # PIV解析をする関数 def piv ( dir = 'input' , out_dir = 'piv_out' , wsize = 32 , overlap = 0 , threshold = 30 ) : path_list = sorted ( glob . glob ( os.path . join ( * [ dir , '*' ] ) ) ) # ファイルパスをソートしてリストする # 全画像ファイルをリストしてPIV計算を実施 for i in range ( len ( path_list ) - 1 ) : vectors_amps = [ ] # ベクトルの大きさ情報を格納する空配列 coordinates = [ ] # 座標情報を格納する空配列 correlation_coef = [ ] # 相関係数情報を格納する空配列 # グレースケール画像で2枚(i, i+1)の画像を読み込み img1 = cv2 . imread ( path_list [ i ] , 0 ) img2 = cv2 . imread ( path_list [ i + 1 ] , 0 ) h , w = img1 . shape h2 , w2 = img2 . shape # ファイルサイズが揃っていない場合はエラー if h != h2 or w != w2 : print ( 'Error:img(i).shape != img(i+1).shape.' ) print ( 'Align image size.' ) w_st = int ( w / ( wsize - overlap ) ) # 幅のストライドを計算 h_st = int ( h / ( wsize - overlap ) ) # 高さのストライドを計算 # 画像を走査しながら処理をする（h_st, w_stのサイズでループ） for k in range ( h_st - 1 ) : for j in range ( w_st - 1 ) : # テンプレートマッチング部分------------------------------------------------------------- # 抽出範囲の点（左上、右上）を計算 p_h1 = k * ( wsize - overlap ) p_h2 = p_h1 + wsize p_w1 = j * ( wsize - overlap ) p_w2 = p_w1 + wsize # パターンマッチングから移動先座標を計算 template = img1 [ p_h1 : p_h2 , p_w1 : p_w2 ] # img[i]からテンプレート画像を抽出 method = cv2 . TM_CCOEFF _ NORMED # NCC(正規化相互相関係数)を選択 res = cv2 . matchTemplate ( img2 , template , method ) # img[i+1]に対するテンプレートマッチングの結果 min_val , max_val , min_loc , max_loc = cv2 . minMaxLoc ( res ) # 最小値と最大値、その位置を取得 # テンプレート画像の中心座標を計算 before_w = int ( p_w1 + ( p_w2 - p_w1 ) / 2 ) # 探査前のx座標 before_h = int ( p_h1 + ( p_h2 - p_h1 ) / 2 ) # 探査前のy座標 after_w = int ( max_loc [ 0 ] + wsize / 2 ) # 探査後のx座標 after_h = int ( max_loc [ 1 ] + wsize / 2 ) # 探査後のy座標 dx = after_w - before _ w # x増分値 dy = after_h - before _ h # y増分値 vector_amp = np . sqrt ( dx * * 2 + dy * * 2 ) # ベクトルの大きさ coordinate = [ before_w , before_h , dx , dy ] # 座標と増分値セット vectors_amps . append ( vector_amp ) coordinates . append ( coordinate ) correlation_coef . append ( max_val ) # フォントの種類とサイズを設定する。 plt . rcParams [ 'font.size' ] = 14 plt . rcParams [ 'font.family' ] = 'Times New Roman' plt . rcParams [ 'xtick.direction' ] = 'in' plt . rcParams [ 'ytick.direction' ] = 'in' # グラフの上下左右に目盛線を付ける。 fig = plt . figure ( ) ax1 = fig . add_subplot ( 111 ) ax1 . yaxis . set_ticks_position ( 'both' ) ax1 . xaxis . set_ticks_position ( 'both' ) # 背景画像の用意と表示設定 ax1 . imshow ( img2 , cmap = 'gray' ) # 背景画像 cm = plt . get_cmap ( 'jet' ) # カラーマップ vsf = 2 # ベクトル表示のスケールファクタ # 誤ベクトル処理（一度に一定ピクセル（threshold）以上のベクトルは長さを0にする） for m in range ( len ( vectors_amps ) ) : print ( m , '/' , len ( vectors_amps ) - 1 ) if vectors_amps [ m ] >= threshold : coordinates [ m ] [ 2 ] = 0 coordinates [ m ] [ 3 ] = 0 vectors_amps [ m ] = 0 # 誤ベクトル処理後のベクトルをMin-Max正規化(cmapでベクトルに色を付けるためだけの変数)

norm_vectors = ( vectors_amps - np . min ( vectors_amps ) ) / ( np . max ( vectors_amps ) - np . min ( vectors_amps ) )

# 画像にベクトルをプロットする for n in range ( len ( vectors_amps ) ) : print ( n , '/' , len ( vectors_amps ) - 1 ) # 進捗表示 # 長さ0以外の場合に図にベクトル(dx, dyにそれぞれスケールを乗じた後のベクトル)を描画 if vectors_amps [ n ] != 0 : ax1 . arrow ( x = coordinates [ n ] [ 0 ] , y = coordinates [ n ] [ 1 ] , dx = vsf * coordinates [ n ] [ 2 ] , dy = vsf * coordinates [ n ] [ 3 ] , width = 0.01 , head_width = 10 , color = cm ( norm_vectors [ n ] ) ) # レイアウトタイト設定 fig . tight_layout ( ) # out_dirで指定したフォルダが無い時に新規作成 if os.path . exists ( out_dir ) : os . mkdir ( out_dir ) # 画像保存パスを連番で準備 path = os.path . join ( * [ out_dir , str ( "" . format ( count ) ) + '.png' ] ) plt . savefig ( path ) # 画像作成の進捗を表示 print ( count , '/' , len ( path_list ) - 1 ) # PIV解析の関数を実行 piv ( dir = 'input' , out_dir = 'piv_out' ) 分布流れの可視化結果

各高さ方向の平均流速をサンプル動画を作成した時の理論流速分布に重ねてみた図を以下に示します。このように今回のデフォルト設定(overlap=0)では分布の位置がずれている事がわかりますが、overlap=12とする事で計測点が増え、さらに山の位置が改善される結果を得ました。これはパターンマッチング時のウィンドウ走査の分解能を変更している事になりますが、overlapを増やすと計算数もどんどん増えていくので、処理が重くなってしまう事に注意が必要です。

オフセット流れの可視化結果今度は上の方の流速が速い流れの可視化をしてみた結果を以下に示します。こちらも狙い通り。一様流れの可視化結果管内に流速分布が無く、一様に流速が発生している場合の可視化結果を以下に示します。画像中全域で同じ色のベクトルを得ました。ただ、どうやら流入付近は計算が苦手のようです。一様流れの場合の速度分布比較は以下。このような流れの場合は特にoverlap設定をしていてもしていなくても、精度的にはそんなに変わらないようですね（定量的で無くて申し訳ないですが…）。 PIV解析後の画像をGIF動画化するコード create gif from PIL import Image import glob # GIFアニメーションを作成 def create_gif ( in_dir , out_filename ) : path_list = sorted ( glob . glob ( os.path . join ( * [ in_dir , '*' ] ) ) ) # ファイルパスをソートしてリストする imgs = [ ] # 画像をappendするための空配列を定義 # ファイルのフルパスからファイル名と拡張子を抽出 for i in range ( len ( path_list ) ) : img = Image . open ( path_list [ i ] ) # 画像ファイルを1つずつ開く imgs . append ( img ) # 画像をappendで配列に格納していく # appendした画像配列をGIFにする imgs [ 0 ] . save ( out_filename , save_all = True , append_images = imgs [ 1 : ] , optimize = False , duration = 100 , loop = 0 ) # GIFアニメーションを作成する関数を実行する create_gif ( in_dir = 'piv_out' , out_filename = 'piv-movie.gif' ) コード解説ストライド w_st = int ( w / ( wsize - overlap ) ) # 幅のストライドを計算 h_st = int ( h / ( wsize - overlap ) ) # 高さのストライドを計算図解すると以下です。今回ウィンドウ幅\(w_, h_\)は同じですが、まずは幅方向にoverlapを考慮しながら走査をしていきます。幅方向終端まで到達した後は縦方向にウィンドウをずらしますが、その際にも同一のoverlapを考慮してストライドさせます。この1画像に対するストライド走査がforループのjとkの部分に相当します。画像走査処理をアニメーションで示したのが以下。オーバーラップ0だと隙間なくタイルが敷き詰められていくようなものです。横方向のみですが、オーバーラップを考慮すると以下のようになります。ウィンドウが重なっている様子が確認できると思います。テンプレートマッチング template matching # パターンマッチングから移動先座標を計算 template = img1 [ p_h1 : p_h2 , p_w1 : p_w2 ] # img[i]からテンプレート画像を抽出 method = cv2 . TM_CCOEFF _ NORMED # NCC(正規化相互相関係数)を選択 res = cv2 . matchTemplate ( img2 , template , method ) # img[i+1]に対するテンプレートマッチングの結果 min_val , max_val , min_loc , max_loc = cv2 . minMaxLoc ( res ) # 最小値と最大値、その位置を取得手法自体は相関係数を用いる NCC ( N ormalized C ross C orrelation)を使っており、式は以下です。

しかし、この部分はOpenCVに任せています。詳しくは「OpenCVのテンプレートマッチングで変形量を算出する方法」の記事をご覧下さい。やっている事は流れの前後で、ウィンドウで抽出した部分がどこからどこに行っているかを相関係数を使って計算しているのみです。下図は流れの前後でウィンドウがしっかりと追従している（流れを追えている）例を示した図です。これを画像走査で画像全域で行います。

誤ベクトル除去 error vector elimination # 誤ベクトル処理（一度に一定ピクセル（threshold）以上のベクトルは長さを0にする） for m in range ( len ( vectors_amps ) ) : print ( m , '/' , len ( vectors_amps ) - 1 ) if vectors_amps [ m ] >= threshold : coordinates [ m ] [ 2 ] = 0 coordinates [ m ] [ 3 ] = 0 vectors_amps [ m ] = 0 誤ベクトル除去の有無を比較した図が以下です。今回はベクトルの色をベクトルの振幅最大値で自動的に付けているので、除去しないと色の分布も偏ってしまいます。ベクトルプロット vector plot # 画像にベクトルをプロットする for n in range ( len ( vectors_amps ) ) : print ( n , '/' , len ( vectors_amps ) - 1 ) # 進捗表示 # 長さ0以外の場合に図にベクトル(dx, dyにそれぞれスケールを乗じた後のベクトル)を描画 if vectors_amps [ n ] != 0 : ax1 . arrow ( x = coordinates [ n ] [ 0 ] , y = coordinates [ n ] [ 1 ] , dx = vsf * coordinates [ n ] [ 2 ] , dy = vsf * coordinates [ n ] [ 3 ] , width = 0.01 , head_width = 10 , color = cm ( norm_vectors [ n ] ) ) 主要なコードはこんな感じです。全体的に全コード内にコメントを残しましたので、詳細はそちらを参照して下さい。

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

まとめ

本記事ではOpenCVを使ってPIVコードを自作してみました。もちろん本格的なPIV計測で使用する市販のソフトと比べると、プログラムの効率も悪いと思うので計算速度も遅いと思います。しかしこんなに簡単なコードでも誤ベクトル除去やベクトルプロットができ、それなりに流れが可視化できる、さらに粒子動画以外でも適用が簡単にできるというのは驚きました。ライトな用途であればPythonで十分かなと思います。今後は本記事のコードを応用させて物体運動の計測に特化したものも作ってみようと思います。

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

SNSでもご購読できます。

wat

機械工学を専攻し大学院を修了後、技術系の職に就き日々実験やシミュレーションを使う仕事をしています。このブログでは初心者が科学技術プログラムを作れるようになることを目標に、学習結果を記録していきます。著書「いきなりプログラミングPython（翔泳社）」 ※資格・計算力学技術者 1級 2級 (振動) ・計算力学技術者 2級 (熱流体) ・基本情報技術者・G検定・AI実装検定 S級 A級 ※Information ・Amazonのアソシエイトとして、当メディアは適格販売により収入を得ています。・当メディアはGoogleによるユーザーに合わせた自動広告が表示されます。

突然のコントで失礼します．御社のサイトを頻繁に閲覧させて頂き，勉強させて頂いているものです．特にこのページのPIVが大変参考になりまして，私自信が開発しているプログラム（以後プログラムA）でも一部使用させて頂いております．その開発プログラムを何らかの形で第三者に配布をしたいと考えているのですが，著作権に関して確認したことがあり，本コメントを書いております．上記のような条件において，御社に対して何かロイヤリティが発生しますでしょうか？またプログラムAの関連するdefで御社のサイトに記載されているcopyrightを記述することで，著作権に対して，その権利を尊重していると考えることができますでしょうか？ちなみに私自信は無償で配布する予定です．著作権に関して勉強中で，非常に稚拙な質問で恐縮です．ご教授頂ければ幸甚です．よろしくお願いします．

いつもご訪問ありがとうございます。当ブログのコードは全て非営利目的、かつ趣味で書いています。そのため、ご自由にコピーしてご活用や配布をして頂いて構いません。特にコードにURLやcopyrightを記載する必要もありません。但し、当ブログのコードはバグの保証や商用利用の際の不都合当は配慮していません（今回は無償とのことですが念のため）。使用している外部ライブラリによっては配布に関する条項がある可能性があるため、OpenCV等の公式ページ、OSS(Open Source Software)の考え方をご確認されることをおすすめ致します。以上、今後とも記事をご活用頂けると幸いです。

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

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

ご訪問ありがとうございます。また、いつも参考にして頂き恐縮です。ベクトルマップを数値データとして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をそのまま文字列や数値として記録するか、それとも成分毎にファイルを分けるのか…といった欲しい形に整形する必要があると思います。

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

コメントを残すコメントをキャンセル

カテゴリー

📎📎📎📎📎📎📎📎📎📎

Без кейворда

本記事の目的

目標

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

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

まとめ

SNSでもご購読できます。

wat

コメント

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

Без кейворда

本記事の目的

目標

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

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

まとめ

SNSでもご購読できます。

wat

コメント

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル