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

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} . Chaque norme est||\cdot||_Est la norme de trace,||\cdot|| _1$Est la norme L1,||\cdot||^2_FEst la norme de FrobeniusA\in\mathbb{R}^{n_1 \times n_2}Composante de singularité de\sigmaensuite

\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. original.png

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

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

Recommended Posts

Séparation d'objets en mouvement avec Robust PCA
PCA avec Scikit-learn
Séparation d'arrière-plan / objet en mouvement à l'aide de la décomposition en mode dynamique
Moyenne mobile avec numpy
Régression linéaire robuste avec scikit-learn