2021.01.31

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

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

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

本記事の目的

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

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

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

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

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

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

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

目標

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

無料で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という所は誤ベクトル処理の閾値です。後ほど説明しますが、流速の速い映像にはこの値を大きくする必要があります。

import cv2
import os
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):
    count = 0
    path_list = sorted(glob.glob(os.path.join(*[dir, '*'])))  # ファイルパスをソートしてリストする

    # 全画像ファイルをリストしてPIV計算を実施
    for i in range(len(path_list)-1):
        count += 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.')
            break

        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):
            pass
        else:
            os.mkdir(out_dir)

        # 画像保存パスを連番で準備
        path = os.path.join(*[out_dir, str("{:05}".format(count)) + '.png'])

        # 画像を保存する
        plt.savefig(path)

        # 画像作成の進捗を表示
        print(count, '/', len(path_list) - 1)

    return

# PIV解析の関数を実行
piv(dir='input', out_dir='piv_out')

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

import cv2

import os

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):

count = 0

path_list = sorted(glob.glob(os.path.join(*[dir, '*']))) # ファイルパスをソートしてリストする

# 全画像ファイルをリストしてPIV計算を実施

for i in range(len(path_list)-1):

count += 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.')

break

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):

pass

else:

os.mkdir(out_dir)

# 画像保存パスを連番で準備

path = os.path.join(*[out_dir, str("{:05}".format(count)) + '.png'])

# 画像を保存する

plt.savefig(path)

# 画像作成の進捗を表示

print(count, '/', len(path_list) - 1)

return

# PIV解析の関数を実行

piv(dir='input', out_dir='piv_out')

分布流れの可視化結果

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

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

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

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

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

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

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

一様流れの可視化結果

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

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

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

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

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

from PIL import Image
import os
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')

from PIL import Image

import os

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')

コード解説

ストライド

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

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

w_st = int(w / (wsize - overlap))                    # 幅のストライドを計算
h_st = int(h / (wsize - overlap))                     # 高さのストライドを計算

1 2	w_st = int(w / (wsize - overlap)) # 幅のストライドを計算 h_st = int(h / (wsize - overlap)) # 高さのストライドを計算

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

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

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

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

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

テンプレートマッチング

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

# パターンマッチングから移動先座標を計算
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)  # 最小値と最大値、その位置を取得

# パターンマッチングから移動先座標を計算

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(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とする事で対応しました（相関係数的には高くてもテンプレートマッチングが失敗しているものもあったので）。

# 誤ベクトル処理（一度に一定ピクセル（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

# 誤ベクトル処理（一度に一定ピクセル（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

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

ベクトルプロット

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

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

# 画像にベクトルをプロットする
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]))

# 画像にベクトルをプロットする

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]))

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

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

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

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

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

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

まとめ

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

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

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

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

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

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

萬矢敦啓より:

2022年4月25日 6:19 PM

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

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

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

返信
1. wat より:
  
  2022年4月26日 7:16 PM
  
  いつもご訪問ありがとうございます。
  当ブログのコードは全て非営利目的、かつ趣味で書いています。
  そのため、ご自由にコピーしてご活用や配布をして頂いて構いません。
  特にコードにURLやcopyrightを記載する必要もありません。
  
  但し、当ブログのコードはバグの保証や商用利用の際の不都合当は配慮していません（今回は無償とのことですが念のため）。
  使用している外部ライブラリによっては配布に関する条項がある可能性があるため、OpenCV等の公式ページ、OSS(Open Source Software)の考え方をご確認されることをおすすめ致します。
  
  以上、今後とも記事をご活用頂けると幸いです。
  
  返信
  1. 萬矢敦啓より:
    
    2022年5月24日 7:18 PM
    
    ご丁寧な回答ありがとうございました．よく理解することができました．
    
    返信
RM より:

2022年12月2日 3:55 PM

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

返信
1. wat より:
  
  2022年12月3日 10:05 PM
  
  ご訪問ありがとうございます。
  また、いつも参考にして頂き恐縮です。
  
  ベクトルマップを数値データとして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 より:
    
    2022年12月5日 8:59 AM
    
    ご丁寧で分かりやすい回答，ありがとうございます．よく理解することができました．参考にさせていただきます．
    
    返信

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

本記事の目的

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

目標

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

Python

OpenCV

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

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

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

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

事前準備

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

全コード

分布流れの可視化結果

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

一様流れの可視化結果

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

コード解説

ストライド

テンプレートマッチング

誤ベクトル除去

ベクトルプロット

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

まとめ

コメント

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

最近の投稿

アーカイブ

カテゴリー

本記事の目的

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

目標

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

Python

OpenCV

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

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

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

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

事前準備

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

全コード

分布流れの可視化結果

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

一様流れの可視化結果

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

コード解説

ストライド

テンプレートマッチング

誤ベクトル除去

ベクトルプロット

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

まとめ

SNSでもご購読できます。

コメント

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

最近の投稿

アーカイブ

カテゴリー

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