【Python/OpenCV】パーティクルフィルタで物体追跡

Python版OpenCVでパーティクルフィルタを実装し、物体追跡する方法をソースコード付きで解説します。

パーティクルフィルタで物体追跡
サンプルコード
サンプルコードの解説
関連ページ

パーティクルフィルタで物体追跡

パーティクルフィルタ(Particle filter)とは、確率分布による時系列データの予測手法です。
粒子フィルターや逐次モンテカルロ法とも呼ばれ、画像処理分野では物体追跡に利用されています。
今回はPython + OpenCV + NumPyでパーティクルフィルタを実装し、赤色の物体を追跡してみました。

【画像処理】パーティクルフィルタによる物体追跡の原理

この記事では、画像処理におけるパーティクルフィルタによる動体追跡の原理について解説します。

実行例

サンプルプログラムの実行結果です。赤色の物体を追跡しています。

サンプルコード

サンプルプログラムのソースコードです。
このプログラムは、動画内の特定の色を追跡し、その結果をリアルタイムで表示するものです。各ステップでパーティクルフィルタの処理を行い、追跡対象の位置を推定します。

サンプルコードの解説

上記コードの各部分の詳細を説明します。

了解しました！では、各部分をさらに詳細に解説します。

1. 必要なライブラリのインポート

import cv2
import numpy as np

cv2: OpenCVライブラリで、画像処理やコンピュータビジョンの機能を提供します。例えば、画像の読み込み、表示、色空間の変換などが可能です。
numpy: 数値計算ライブラリで、配列操作を効率的に行います。例えば、画像データの操作や数値計算に使用します。

2. 色範囲の判定

def is_target(roi):
    return (roi <= 30) | (roi >= 150)

is_target関数は、指定された色範囲（Hue値が30以下または150以上）に該当するかを判定します。これは、追跡対象の色を特定するために使用されます。例えば、赤色の物体を追跡する場合、赤色のHue値がこの範囲に入ることが多いです。

3. 最大ブロブの中心座標を算出

def max_moment_point(mask):
    label = cv2.connectedComponentsWithStats(mask)
    data = np.delete(label[2], 0, 0)
    center = np.delete(label[3], 0, 0)
    moment = data[:, 4]
    max_index = np.argmax(moment)
    return center[max_index]

max_moment_point関数は、マスク画像から最大のブロブ（連結成分）の中心座標を計算します。
- cv2.connectedComponentsWithStats: 画像のラベリングを行い、各ブロブの統計情報を取得します。
- np.delete: 最初のラベル（背景）を削除します。
- np.argmax: 最大の面積を持つブロブのインデックスを取得します。
- center[max_index]: 最大ブロブの中心座標を返します。

4. パーティクルの初期化

def initialize(img, N):
    mask = img.copy()
    mask[is_target(mask) == False] = 0
    x, y = max_moment_point(mask)
    w = calc_likelihood(x, y, img)
    ps = np.ndarray((N, 3), dtype=np.float32)
    ps[:] = [x, y, w]
    return ps

initialize関数は、画像内の追跡対象の色を持つ領域をマスクし、最大ブロブの中心座標を初期パーティクルとして設定します。
- mask[is_target(mask) == False] = 0: 追跡対象外の色を持つピクセルを0に設定します。
- max_moment_point(mask): 最大ブロブの中心座標を取得します。
- calc_likelihood(x, y, img): 尤度を計算します。
- ps = np.ndarray((N, 3), dtype=np.float32): パーティクル格納用の配列を生成します。
- ps[:] = [x, y, w]: パーティクル用配列に中心座標と尤度をセットします。

5. リサンプリング

def resampling(ps):
    ws = ps[:, 2].cumsum()
    last_w = ws[ws.shape[0] - 1]
    new_ps = np.empty(ps.shape)
    for i in range(ps.shape[0]):
        w = np.random.rand() * last_w
        new_ps[i] = ps[(ws > w).argmax()]
        new_ps[i, 2] = 1.0
    return new_ps

resampling関数は、パーティクルの重みに基づいて新しいパーティクルを生成します。
- ps[:, 2].cumsum(): パーティクルの重みの累積和を計算します。
- np.random.rand() * last_w: ランダムな重みを生成します。
- ps[(ws > w).argmax()]: 累積重みがランダムな重みを超える最初のパーティクルを選択します。
- new_ps[i, 2] = 1.0: 新しいパーティクルの重みを1.0に設定します。

6. パーティクルの位置予測

def predict_position(ps, var=13.0):
    ps[:, 0] += np.random.randn((ps.shape[0])) * var
    ps[:, 1] += np.random.randn((ps.shape[0])) * var

predict_position関数は、パーティクルの位置をランダムに少しずらします。
- np.random.randn((ps.shape[0])) * var: 正規分布に従ってランダムな値を生成し、位置をずらします。

7. 尤度の計算

def calc_likelihood(x, y, img, w=30, h=30):
    x1, y1 = max(0, x-w/2), max(0, y-h/2)
    x2, y2 = min(img.shape[1], x+w/2), min(img.shape[0], y+h/2)
    x1, y1, x2, y2 = int(x1), int(y1), int(x2), int(y2)
    roi = img[y1:y2, x1:x2]
    count = roi[is_target(roi)].size
    return (float(count) / img.size) if count > 0 else 0.0001

calc_likelihood関数は、指定された座標を中心とする領域内の追跡対象の色の割合を計算します。
- max(0, x-w/2): 領域の左上のx座標を計算します。
- min(img.shape[1], x+w/2): 領域の右下のx座標を計算します。
- roi[is_target(roi)].size: 領域内の追跡対象の色のピクセル数を計算します。
- (float(count) / img.size): 尤度を計算します。

8. パーティクルの重み付け

def calc_weight(ps, img):
    for i in range(ps.shape[0]):
        ps[i][2] = calc_likelihood(ps[i, 0], ps[i, 1], img)
    ps[:, 2] *= ps.shape[0] / ps[:, 2].sum()

calc_weight関数は、各パーティクルの尤度に基づいて重みを計算し、正規化します。
- calc_likelihood(ps[i, 0], ps[i, 1], img): 各パーティクルの尤度を計算します。
- ps[:, 2] *= ps.shape[0] / ps[:, 2].sum(): 重みを正規化します。

9. 観測

def observer(ps, img):
    calc_weight(ps, img)
    x = (ps[:, 0] * ps[:, 2]).sum()
    y = (ps[:, 1] * ps[:, 2]).sum()
    return (x, y) / ps[:, 2].sum()

observer関数は、重み付き平均を計算してパーティクルの推定位置を返します。
- calc_weight(ps, img): パーティクルの重みを計算します。
- (ps[:, 0] * ps[:, 2]).sum(): 重み付きのx座標の合計を計算します。
- (ps[:, 1] * ps[:, 2]).sum(): 重み付きのy座標の合計を計算します。
- (x, y) / ps[:, 2].sum(): 重み付き平均を計算します。

10. パーティクルフィルタ

def particle_filter(ps, img, N=300):
    if ps is None:
        ps = initialize(img, N)
    ps = resampling(ps)
    predict_position(ps)
    x, y = observer(ps, img)
    return ps, int(x), int(y)

particle_filter関数は、パーティクルフィルタの主要な処理をまとめた関数です。
- initialize(img, N): パーティクルを初期化します。
- resampling(ps): パーティクルをリサンプリングします。
- predict_position(ps): パーティクルの位置を予測します。

11. メイン関数

①パーティクルの初期化

ps = None

ps: パーティクルを格納する変数です。最初は None に設定されています。

②動画ファイルのキャプチャ

cap = cv2.VideoCapture("C:/github/sample/python/opencv/dataset/videos/red_marker.mp4")

cv2.VideoCapture: 指定されたパスの動画ファイルを読み込みます。cap は動画キャプチャオブジェクトです。

③メインループ

while cv2.waitKey(30) < 0:

cv2.waitKey(30): 30ミリ秒待機し、キー入力を待ちます。キーが押されるとループを終了します。

④フレームの読み込み

ret, frame = cap.read()

cap.read(): 動画からフレームを読み込みます。ret はフレームが正常に読み込まれたかどうかを示すブール値で、frame は読み込まれたフレームです。

⑤色空間の変換

hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV_FULL)
h = hsv[:, :, 0]

cv2.cvtColor: フレームをBGR色空間からHSV色空間に変換します。hsv は変換後の画像で、h はそのHueチャンネルです。

⑥SとVチャンネルの二値化

ret, s = cv2.threshold(hsv[:, :, 1], 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
ret, v = cv2.threshold(hsv[:, :, 2], 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
h[(s == 0) | (v == 0)] = 100

cv2.threshold: S（彩度）とV（明度）チャンネルを二値化します。大津の手法を用いて自動的に閾値を決定します。
h[(s == 0) | (v == 0)] = 100: SまたはVが0のピクセルをHue値100に設定します。これにより、追跡対象外の領域を除外します。

⑦パーティクルフィルタの適用

ps, x, y = particle_filter(ps, h, 300)

particle_filter: パーティクルフィルタを適用し、パーティクルの位置を更新します。ps は更新されたパーティクル、x と y は推定された追跡対象の中心座標です。

⑧パーティクルの範囲チェック

if ps is None:
    continue
ps1 = ps[(ps[:, 0] >= 0) & (ps[:, 0] < frame.shape[1]) & (ps[:, 1] >= 0) & (ps[:, 1] < frame.shape[0])]

ps is None: パーティクルが初期化されていない場合は次のループに進みます。
ps1: 画像の範囲内にあるパーティクルのみを抽出します。

⑨パーティクルの描画

for i in range(ps1.shape[0]):
    frame[int(ps1[i, 1]), int(ps1[i, 0])] = [0, 0, 200]

frame[int(ps1[i, 1]), int(ps1[i, 0])] = [0, 0, 200]: パーティクルの位置を赤色で塗りつぶします。

⑩追跡対象の描画

cv2.rectangle(frame, (x-20, y-20), (x+20, y+20), (0, 0, 200), 5)

cv2.rectangle: 推定された追跡対象の中心座標を囲む赤い矩形を描画します。

⑪結果の表示

cv2.imshow('Result', frame)

cv2.imshow: 処理結果を表示します。

⑫リソースの解放

cap.release()
cv2.destroyAllWindows()

cap.release(): 動画キャプチャオブジェクトを解放します。
cv2.destroyAllWindows(): すべてのOpenCVウィンドウを閉じます。