Séparation d'objets en mouvement avec Robust PCA
introduction
Dans l'article précédent (https://qiita.com/matsxxx/items/652e58f77558faecfd23), j'ai montré un exemple de séparation d'objets en mouvement dans une vidéo avec la décomposition en mode dynamique (DMD). Dans cet article, nous allons essayer de déplacer la séparation des objets avec l'analyse en composantes principales robuste (PCA robuste) de l'algorithme PCP (Principal Component Pursuit) [1].
environnement
- Windows 10 home edition
- Python (3.7.6)
- Numpy(1.18.5)
- Opencv(4.1.1)
Robust PCA via PCP
L'ACP robuste est une méthode de décomposition en une structure de bas rang et une structure clairsemée. Dans le cas d'un corps mobile, le fond a une structure de bas rang et le corps mobile a une structure clairsemée. Cette fois, j'essaierai Robust PCA, qui est un algorithme PCP. Décompose la matrice de données originale $ M $ en une matrice de rang bas $ L $ et une matrice creuse $ S $.
M = L + S
PCP se décompose en $ L $ et $ S $ en minimisant le lagrangien étendu suivant $ l $.
l(L, S, Y) = \| L \|_* + \lambda\|S\|_1 +<Y, M - L - S> + \frac{\mu}{2}\|M-L-s\|^2_F
$ L $ est une matrice de bas rang, $ S $ est une matrice creuse, $ Y $ est une matrice multiplicatrice de Lagrange, et $ \ lambda $ et $ \ mu $ sont des paramètres PCP. $ \ <, * > $ Est le produit interne. Les dimensions de chaque variable sont $ M, L, S, Y \ in \ mathbb {R} ^ {n_1 \ times n_2}
\begin{align}
&||A||_* = \mathrm{tr}(\sqrt{A^\mathrm{T}A})= \sum^{\mathrm{min}(n_1,n_2)}_{i=1} \sigma_i\\
&||A||_1 = \underset{1\leq j \leq n_2}{\mathrm{max}}\sum^{n_1}_{i=1}|a_{ij}|\\
&||A||_F = \sqrt{\sum^{n_1}_{i=1}\sum^{n_2}_{j=1}|a_{ij}|^2} = \sqrt{\mathrm{tr}(A^\mathrm{T}A)} = \sqrt{\sum^{\mathrm{min}(n_1,n_2)}_{i=1}\sigma^2_i}
\end{align}
Ce sera. L'indice $ \ mathrm {T} $ représente une matrice transposée. La formule qui décrit la norme est [Description Wikipedia](https://ja.wikipedia.org/wiki/%E8%A1%8C%E5%88%97%E3%83%8E%E3%83%AB% J'ai fait référence à E3% 83% A0). La matrice creuse $ S $ est mise à jour avec l'opérateur de retrait, et la matrice de rang bas $ L $ est mise à jour avec l'opérateur de seuillage de valeur singulière. L'opérateur de retrait $ S_ \ tau $ du paramètre $ \ tau $ est
\begin{align}
&S_\tau(x) = \mathrm{sgn}(x) \mathrm{max}(|x|-\tau, 0)
\end{align}
Il est exprimé comme. L'opérateur de seuillage de valeur singulière $ D_ \ tau $ utilise l'oprateur de retrait $ S_ \ tau $ et la décomposition de valeur singulière.
D_{\tau}(X) = US_\tau(\Sigma)V^\mathrm{T}, \
\mathrm{where} \ X = U\Sigma V^\mathrm{T}
Ce sera. $ U, \ Sigma et V $ sont respectivement le vecteur singulier gauche, la valeur singulière et le vecteur singulier droit de $ X $. Utilisez l'opérateur de retrait $ S_ \ tau $ et l'opérateur de seuillage de valeur singulière $ D_ \ tau $ pour mettre à jour la matrice de rang bas $ L $ et la matrice creuse $ S $. La formule de mise à jour est écrite comme $ X_k $ avant la mise à jour de la variable $ X $ et $ X_ {k + 1} $ après la mise à jour, et les paramètres PCP $ \ tau, \ mu $ sont saisis et écrits.
\begin{align}
&L_{k+1} = D_{1/\mu}(M - S_k - \mu^{-1}Y_k) \\
&S_{k+1} = S_{\lambda/\mu}(M - L_{k+1} + \mu^{-1}Y_k)
\end{align}
Ce sera. Le symbole $ S $ est déroutant, mais il suit la notation de Ref. [1], sauf que la matrice contingente est paraphrasée comme une matrice transposée. De plus, la méthode de définition des paramètres $ \ mu et \ lambda $ de l'opérateur de retrait $ S_ \ tau $ et de l'opérateur de seuillage de valeur singulière $ D_ \ tau $ est erronée dans l'article PCP 2009 placé dans arXiv. , Attention (la version 2009 de arXiv est $ D_ \ mu, S_ {\ lambda \ mu} $, mais elle devrait être $ D_ {1 / \ mu}, S_ {\ lambda / \ mu} Cela ressemble à $). La matrice multiplicatrice de Lagrange $ Y $ est calculée par la méthode étendue de Lagrange.
Y_{k+1} = Y_k + \mu(M - L_{k+1} - S_{k+1})
Et mettre à jour. Les paramètres PCP $ \ lambda $ et $ \ mu $ sont, dans la référence [1], comme la matrice de données $ M \ in \ mathbb {R} ^ {n_1 \ times n_2} $, respectivement.
\begin{align}
&\lambda = \frac{1}{\sqrt{\mathrm{max}(n_1,n_2)}} \\
&\mu = \frac{n_1 n_2}{4\|M\|_1}
\end{align}
Est réglé. Ce qui suit est un résumé de la procédure de calcul du PCP.
\begin{align}
&\mathrm{initialize:} S_0 = 0, Y_0 =0 \\
& \mathrm{while}\ \mathrm{not}\ \mathrm{converged}:\\
&\qquad L_{k+1} = D_{1/\mu}(M - S_k + \mu^{-1}Y_k) \\
&\qquad S_{k+1} = S_{\lambda/\mu}(M - L_{k+1} + \mu^{-1}Y_k) \\
&\qquad Y_{k+1} = Y_k + \mu(M - L_{k+1} - S_{k+1}) \\
&\mathrm{end} \ \mathrm{while} \\
&\mathrm{output:}L,S.
\end{align}
La convergence est déterminée par la norme de Frobenius suivante.
\|M -L - S\|_F \leq \delta\|M\|_F \quad \mathrm{with} \ \delta=10^{-7}
Code PCA robuste de l'algorithme PCP ci-dessus en Python.
import numpy as np
import numpy.linalg as LA
def rpca(M, max_iter=800,p_interval=50):
def shrinkage_operator(x, tau):
return np.sign(x) * np.maximum((np.abs(x) - tau), np.zeros_like(x))
def svd_thresholding_operator(X, tau):
U, S, Vh = LA.svd(X, full_matrices=False)
return U @ np.diag(shrinkage_operator(S, tau)) @ Vh
i = 0
S = np.zeros_like(M)
Y = np.zeros_like(M)
error = np.Inf
tol = 1e-7 * LA.norm(M, ord="fro")
mu = M.shape[0] * M.shape[1]/(4 * LA.norm(M, ord=1))
mu_inv = 1/mu
lam = 1/np.sqrt(np.max(M.shape))
while i < max_iter:
L = svd_thresholding_operator(M - S + mu_inv * Y, mu_inv)
S = shrinkage_operator(M - L + mu_inv * Y, lam * mu_inv)
Y = Y + mu * (M - L - S)
error = LA.norm(M - L - S, ord='fro')
if i % p_interval == 0:
print("step:{} error:{}".format(i, error))
if error <= tol:
print("converted! error:{}".format(error))
break
i+=1
else:
print("Not converged")
return L, S
Vidéo
Article précédent De même, en utilisant l'Atrium de l'ensemble de données dans ici Je vais. J'utilise 120frame à 170frame. Utilisez Robust PCA pour séparer la personne qui marche de l'arrière-plan.
import cv2
def prepare_frames_to_use(video_path="./atrium_video.avi",#Chemin de l'image
start_frame=120,#Première image
end_frame=170,#Dernière image
r=4,#Paramètres de réduction d'image
):
#Chargement d'image
cap = cv2.VideoCapture(video_path)
#Obtenez la résolution d'image, la fréquence d'images, le nombre d'images
wid = cap.get(cv2.CAP_PROP_FRAME_WIDTH)#côté
wid_r = int(wid/r)
hei = cap.get(cv2.CAP_PROP_FRAME_HEIGHT)#Verticale
hei_r = int(hei/r)
fps = cap.get(cv2.CAP_PROP_FPS)#fréquence d'images
count = cap.get(cv2.CAP_PROP_FRAME_COUNT)#Nombre de cadres
dt = 1/fps#Nombre de secondes entre les images
print(f"wid:{wid}", f" hei:{hei}", f" fps:{fps}", f" count:{count}", f"dt:{dt}")
#Extraction de l'image cible
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,(wid_r, hei_r)),
cv2.COLOR_BGR2GRAY).flatten()#
cat_frames.append(buf)
cap.release()
cat_frames = np.array(cat_frames).T
return cat_frames, wid_r, hei_r, fps
Séparation des objets en mouvement par PCA robuste
Utilisez Robust PCA pour trouver des matrices de bas rang et de faible densité. Entrez simplement la matrice de données vidéo dans la fonction rpca.
def detect_moving_object():
#Acquisition des cadres utilisés pour le calcul
frames, wid, hei, fps = prepare_frames_to_use()
#Exécutez un PCA robuste
L, S = rpca(frames)
#Sortie vidéo de structure clairsemée
fourcc = cv2.VideoWriter_fourcc(*"mp4v")
writer = cv2.VideoWriter("sparse_video.mp4",
fourcc, fps, (int(wid), int(hei)))
for i in range(S.shape[1]):
s_buf = S[:,i]
#Réglage de la luminosité
lum_min = np.abs(s_buf.min())
lum_max = np.abs(s_buf.max())
lum_w = lum_max + lum_min
s_buf = (s_buf + lum_min)/lum_w * 255
s_buf = cv2.cvtColor(s_buf.reshape(int(hei), -1).astype('uint8'),
cv2.COLOR_GRAY2BGR)
writer.write(s_buf)
writer.release()
résultat
Si vous regardez la vidéo de la procession clairsemée dans la figure de droite, vous pouvez voir que les gens qui marchent sont clairement séparés.
Références
- [1]: Candes, E. J., Li, X., Ma, Y., and Wright, J. 2011. Robust principal component analysis? J. ACM 58, 3, ` Article 11 (May 2011), 37 pages.
- [2]: 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.
Tout le code source
import numpy as np
import numpy.linalg as LA
def rpca(M, max_iter=800,p_interval=50):
def shrinkage_operator(x, tau):
return np.sign(x) * np.maximum((np.abs(x) - tau), np.zeros_like(x))
def svd_thresholding_operator(X, tau):
U, S, Vh = LA.svd(X, full_matrices=False)
return U @ np.diag(shrinkage_operator(S, tau)) @ Vh
i = 0
S = np.zeros_like(M)
Y = np.zeros_like(M)
error = np.Inf
tol = 1e-7 * LA.norm(M, ord="fro")
mu = M.shape[0] * M.shape[1]/(4 * LA.norm(M, ord=1))
mu_inv = 1/mu
lam = 1/np.sqrt(np.max(M.shape))
while i < max_iter:
L = svd_thresholding_operator(M - S + mu_inv * Y, mu_inv)
S = shrinkage_operator(M - L + mu_inv * Y, lam * mu_inv)
Y = Y + mu * (M - L - S)
error = LA.norm(M - L - S, ord='fro')
if i % p_interval == 0:
print("step:{} error:{}".format(i, error))
if error <= tol:
print("converted! error:{}".format(error))
break
i+=1
else:
print("Not converged")
return L, S
import cv2
def prepare_frames_to_use(video_path="./atrium_video.avi",#Chemin de l'image
start_frame=120,#Première image
end_frame=170,#Dernière image
r=4,#Paramètres de réduction d'image
):
#Chargement d'image
cap = cv2.VideoCapture(video_path)
#Obtenez la résolution d'image, la fréquence d'images, le nombre d'images
wid = cap.get(cv2.CAP_PROP_FRAME_WIDTH)#côté
wid_r = int(wid/r)
hei = cap.get(cv2.CAP_PROP_FRAME_HEIGHT)#Verticale
hei_r = int(hei/r)
fps = cap.get(cv2.CAP_PROP_FPS)#fréquence d'images
count = cap.get(cv2.CAP_PROP_FRAME_COUNT)#Nombre de cadres
dt = 1/fps#Nombre de secondes entre les images
print(f"wid:{wid}", f" hei:{hei}", f" fps:{fps}", f" count:{count}", f"dt:{dt}")
#Extraction de l'image cible
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,(wid_r, hei_r)),
cv2.COLOR_BGR2GRAY).flatten()#
cat_frames.append(buf)
cap.release()
cat_frames = np.array(cat_frames).T
return cat_frames, wid_r, hei_r, fps
def detect_moving_object():
#Acquisition des cadres utilisés pour le calcul
frames, wid, hei, fps = prepare_frames_to_use()
#Exécutez un PCA robuste
L, S = rpca(frames)
#Sortie vidéo de structure clairsemée
fourcc = cv2.VideoWriter_fourcc(*"mp4v")
writer = cv2.VideoWriter("sparse_video.mp4",
fourcc, fps, (int(wid), int(hei)))
for i in range(S.shape[1]):
s_buf = S[:,i]
#Réglage de la luminosité
lum_min = np.abs(s_buf.min())
lum_max = np.abs(s_buf.max())
lum_w = lum_max + lum_min
s_buf = (s_buf + lum_min)/lum_w * 255
s_buf = cv2.cvtColor(s_buf.reshape(int(hei), -1).astype('uint8'),
cv2.COLOR_GRAY2BGR)
writer.write(s_buf)
writer.release()
if __name__ == "__main__":
detect_moving_object()