Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@matsxxx

動的モード分解を使った背景・動体分離

Last updated at Posted at 2020-10-15

はじめに

　動的モード分解（Dynamic Mode Decomposition: DMD)を使い、動画の背景を分離して、動体の動的モードのみの動画を作成します。
　動的モード分解は、主成分分析とフーリエ変換を合わせたような方法です。動的モード分解の詳細解説は、以前書いた記事にリンクがあるので、リンク先を参照ください。
-https://qiita.com/matsxxx/items/5e4b272de821fb1c11e0

環境

  • Windows 10 home (CPU:Core i7 Mem:16GB)
  • Python(3.7.6)
  • Opencv(4.1.1)
  • Numpy(1.18.5)
  • Scipy(1.5.2)

動的モード分解とは

　簡単に動的モード分解の手順を紹介をします。動的モード分解は、時系列データから遷移行列の固有値、固有ベクトルを求めて、動的モードを作り出すものです。次元と時間方向の特徴に分解できます。次元方向の特徴は、固有ベクトルに、時間方向の特徴は、固有値の複素数に現れます。

　動的モード分解では、遷移行列の固有値、固有ベクトルが、現実的な計算量で求められるよう、実装されています。本記事では、Exact DMDで実装しています。
dmd_describe.png

import scipy.linalg as la
#DMD定義
def dmd(X, Y, truncate=None):#X=X_{1:n-1} Y=X_{2:n}
    u2,sig2,vh2 = la.svd(X, False) 
    r = len(sig2) if truncate is None else truncate
    u = u2[:,:r]
    sig = np.diag(sig2)[:r,:r]
    v = vh2.conj().T[:,:r]
    Atil = np.dot(np.dot(np.dot(u.conj().T, Y), v), la.inv(sig))
    mu,w = la.eig(Atil)
    phi = np.dot(np.dot(np.dot(Y, v), la.inv(sig)), w)#DMD mode
    return mu, phi

動画

　ここの動画のAtriumを利用しました。歩いている人がいる動画です。120frameから170frameを利用しています。歩いている人と背景を分離します。
original.png

import cv2
#画像のパス
video_path = "./atrium_video.avi"
#画像読み込み
cap = cv2.VideoCapture(video_path)

#画像の解像度、フレームレート、フレーム数取得
wid = cap.get(cv2.CAP_PROP_FRAME_WIDTH)#横
hei = cap.get(cv2.CAP_PROP_FRAME_HEIGHT)#縦
fps = cap.get(cv2.CAP_PROP_FPS)#フレームレート
count = cap.get(cv2.CAP_PROP_FRAME_COUNT)#フレーム数
dt = 1/fps#フレーム間秒数
print(f"wid:{wid}", f" hei:{hei}", f" fps:{fps}", f" count:{count}", f"dt:{dt}")

#利用するフレーム
start_frame =120
end_frame = 170

#画像解像度を1/4にする（計算量削減のため）
r = 4

#フレーム抜き出し
cat_frames = []
cap.set(cv2.CAP_PROP_POS_FRAMES,start_frame)
for i in range(end_frame - start_frame):
    ret, img = cap.read()
    if not ret:
        print("no image")
        break
    buf = cv2.cvtColor(cv2.resize(img,(int(wid/r), int(hei/r))), cv2.COLOR_BGR2GRAY).flatten()#
    cat_frames.append(buf)
cat_frames = np.array(cat_frames).T#DMDに利用する画像
cap.release()

動画の動的モード分解と背景・動体分離

　動画を動的モード分解します。背景は、振動・振幅モードが振幅０、周波数０であるものとしています。動体は、背景以外のものとしています。

#DMDの計算
X1 = cat_frames[:,:-1]
X2 = cat_frames[:,1:]
mu, phi = dmd(X1,X2)
omega = np.log(mu)/dt#振動・振幅モード

#背景と動体の判定
bg = np.abs(omega) < 1e-2 #背景の抽出
fg = ~bg #動体

phi_bg = phi[:,bg]#背景の動的モード
phi_fg = phi[:,fg]#動体の動的モード
omega_bg = omega[bg]#背景の振動・振幅モード
omega_fg = omega[fg]#動体の振動・振幅モード

背景成分以外の動的モードの再構成

　以下の式で、動体モードの動画を再構成します。

$$
X_{dmd}^{fg} = \sum_{k=1}^{n}b_k^{fg}\phi_{k}^{fg}\exp(\omega_k^{fg}t)=\Phi^{fg} diag(\exp(\omega^{fg} t))\mathbf{b}^{fg}
$$
$\Phi^{fg}$ は動体の動的モード行列、$\omega^{fg}$は動体の振動・振幅モード、 動体の動的モード行列の疑似逆行列を、右から$\mathbf{b}^{fg}$の動画初期値に行列積をとった行列です。

#再構築
phi_fg_pinv = la.pinv(phi_fg)#疑似逆行列。かなり時間がかかる。メモリが足りない場合は、rを大きくしてください。

X_fg = np.zeros((len(omega_fg), end_frame - start_frame), dtype='complex')
b = phi_fg_pinv @ cat_frames[:,0]#初期値

for tt in range(end_frame - start_frame):
    X_fg[:,tt] = b * np.exp(omega_fg * dt * tt)
X_fg = phi_fg @ X_fg

#輝度調整用
lum_max = np.abs(X_fg.real.max())
lum_min = np.abs(X_fg.real.min())
lum_diff = lum_max + lum_min

#動画作成
fourcc = cv2.VideoWriter_fourcc(*"mp4v")
writer = cv2.VideoWriter("out_dmd_fg.mp4", fourcc, fps, (int(wid/r), int(hei/r)))

for tt in range(end_frame - start_frame):
    a = X_fg[:,tt].real.reshape(int(hei/r), -1)
    a = (a + lum_min)/lum_diff * 255
    a = a.astype("uint8")
    out_img = np.tile(a, (3,1,1)).transpose((1,2,0))#グレースケールを[wid, hei, 3]の行列に変換する
    writer.write(out_img)
writer.release()

結果

　背景がなくなり、人だけの動きとなっていることが確認できます。ただ、人の動きに残像が見えます。利用した動画の画像が、50枚しかないせいもあると思いますが、直線的な動きの分解は苦手なようです。
original　DMD

参考文献

  • Jodoin, J.-P., Bilodeau, G.-A., Saunier, N., Urban Tracker: Multiple Object Tracking in Urban Mixed Traffic, Accepted for IEEE Winter conference on Applications of Computer Vision (WACV14), Steamboat Springs, Colorado, USA, March 24-26, 2014.
  • JN Kutz, SL Brunton, BW Brunton, JL Proctor, "Dynamic mode decomposition: data-driven modeling of complex systems", 2016.
  • Tu, Rowley, Luchtenburg, Brunton, and Kutz, "On Dynamic Mode Decomposition: Theory and Applications", *American Institute of Mathematical Sciences, 2014.
  • J. Grosek and J.N. Kutz. "Dynamic mode decomposition for real-time background/foreground separation in video.", 2014.

全ソースコード

import numpy as np
import scipy.linalg as la
import cv2

#画像のパス
video_path = "./atrium_video.avi"

#DMD定義
def dmd(X, Y, truncate=None):
    u2,sig2,vh2 = la.svd(X, False) 
    r = len(sig2) if truncate is None else truncate
    u = u2[:,:r]
    sig = np.diag(sig2)[:r,:r]
    v = vh2.conj().T[:,:r]
    Atil = np.dot(np.dot(np.dot(u.conj().T, Y), v), la.inv(sig))
    mu,w = la.eig(Atil)
    phi = np.dot(np.dot(np.dot(Y, v), la.inv(sig)), w)#DMD mode
    return mu, phi

#画像読み込み
cap = cv2.VideoCapture(video_path)

#画像の解像度、フレームレート、フレーム数取得
wid = cap.get(cv2.CAP_PROP_FRAME_WIDTH)#横
hei = cap.get(cv2.CAP_PROP_FRAME_HEIGHT)#縦
fps = cap.get(cv2.CAP_PROP_FPS)#フレームレート
count = cap.get(cv2.CAP_PROP_FRAME_COUNT)#フレーム数
dt = 1/fps#フレーム間秒数
print(f"wid:{wid}", f" hei:{hei}", f" fps:{fps}", f" count:{count}", f"dt:{dt}")

#利用するフレーム
start_frame =120
end_frame = 170

#画像解像度を1/4にする（計算量削減のため）
r = 4

#フレーム抜き出し
cat_frames = []
cap.set(cv2.CAP_PROP_POS_FRAMES,start_frame)
for i in range(end_frame - start_frame):
    ret, img = cap.read()
    if not ret:
        print("no image")
        break
    buf = cv2.cvtColor(cv2.resize(img,(int(wid/r), int(hei/r))), cv2.COLOR_BGR2GRAY).flatten()#
    cat_frames.append(buf)
cat_frames = np.array(cat_frames).T
cap.release()

#DMDの計算
X1 = cat_frames[:,:-1]
X2 = cat_frames[:,1:]
mu, phi = dmd(X1,X2)
omega = np.log(mu)/dt

#背景と動体以外判定
bg = np.abs(omega) < 1e-2 #背景の抽出
fg = ~bg #動体の抽出

phi_bg = phi[:,bg]#背景の動的モード
phi_fg = phi[:,fg]#動体の動的モード
omega_bg = omega[bg]#背景の振動・振幅モード
omega_fg = omega[fg]#動体の振動・振幅モード

#動的モードの再構成
phi_fg_pinv = la.pinv(phi_fg)#疑似逆行列。メモリエラーが出る場合はrを大きくする

X_fg = np.zeros((len(omega_fg), end_frame - start_frame), dtype='complex')
b = phi_fg_pinv @ cat_frames[:,0]#初期値

for tt in range(end_frame - start_frame):
    X_fg[:,tt] = b * np.exp(omega_fg * dt * tt)
X_fg = phi_fg @ X_fg

#輝度調整用
lum_max = np.abs(X_fg.real.max())
lum_min = np.abs(X_fg.real.min())
lum_diff = lum_max + lum_min

#動画作成
fourcc = cv2.VideoWriter_fourcc(*"mp4v")
writer = cv2.VideoWriter("out_dmd_fg.mp4", fourcc, fps, (int(wid/r), int(hei/r)))

for tt in range(end_frame - start_frame):
    a = X_fg[:,tt].real.reshape(int(hei/r), -1)
    a = (a + lum_min)/lum_diff * 255
    a = a.astype("uint8")
    out_img = np.tile(a, (3,1,1)).transpose((1,2,0))#グレースケールを[wid, hei, 3]の行列に変換する
    writer.write(out_img)
writer.release()
3
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?