Patent Publication Number: US-9892520-B2

Title: Method and system for segmenting pedestrian flows in videos

Description:
FIELD OF THE INVENTION 
     This invention relates generally to computer vision, an more particularly segmenting pedestrians in flows as depicted in videos. 
     BACKGROUND OF THE INVENTION 
     The segmentation of motion flows within dense crowds of pedestrians in videos is an essential tool in crowd safety and crowd control applications. Videos of crowded scenes can exhibit complex crowd behaviors even under normal situations. For example, crowd flows in large congested areas, such as train stations, can initially appear chaotic. However, is often the case that low dimensional dynamical structures exist in the flow that is desirable to identify and segment from unstructured flows. Moreover, the automatic segmentation of independent crowd flows aids in the monitoring and prediction of hazardous situations in crowded environments. 
     Of particular interest is the detection and estimation of crowd flows using motion information extracted from the videos. Motion vectors can be determined using optical flow estimation applied to texture, that is pixel intensities, in videos, or the motion vectors can be extracted directly from a bitstream. The bitstream can encoded using any of the well known coding standards, e.g., MPEG, H.264, HEVC, etc. 
     Considering pedestrians in a crowded scene as particles in a flow, the motion vectors of a video frame correspond to observations of the velocities of particles at a time instance in the flow. Processing motion vectors, instead of the video texture, protects the privacy of individuals observed in surveillance videos. 
     U.S. Pat. No. 8,773,536 discloses a method for detecting independent motion in a surveillance video. The method constructs a linear system from macroblocks of the video by comparing texture gradients with motion vector. Independent flows are detected when the motion is flagged as a statistical outlier relative to the linear system. 
     U.S. Pat. No. 8,358,806 discloses a method for crowd segmentation in a video using shape indexing. Background differencing is performed on the video to identify a foreground silhouette shape. Approximate numbers and positions of people are determined by matching the foreground silhouette shape against a set of predetermined foreground silhouette shapes. 
     Dynamical System Modeling 
     When the density of a crowd of pedestrians is high, the motion of individuals in the crowd can be modeled as a fluid flow. One such model that is commonly used for crowd analysis is the Hughes model, see: Hughes, “A continuum theory for the flow of pedestrians,” Transportation Research Part B: Methodological, Volume 36, Issue 6, Pages 507-535, July 2002. 
     Hughes models a crowd flow as a function of density ρ(x, y, t), and velocities (u(x, y, t) and v(x, y, t)) as 
                           ∂     ∂   t       ⁢     ρ   ⁡     (     x   ,   y   ,   t     )         +       ∂     ∂   x       ⁢     (       ρ   ⁡     (     x   ,   y   ,   t     )       ⁢     u   ⁡     (     x   ,   y   ,   t     )         )       +       ∂     ∂   y       ⁢     (       ρ   ⁡     (     x   ,   y   ,   t     )       ⁢     v   ⁡     (     x   ,   y   ,   t     )         )         =   0     ,           (   1   )               
where u(x, y, t) and v(x, y, t) are the respective velocities in the horizontal and vertical directions of every spatial point (x, y) and time t. A Greenshields model can also be used to relate the density and velocity fields in crowd modeling. The Greenshields model is
 
                       u   ⁡     (     x   ,   y   ,   t     )       =       u   _     ⁡     (     1   -       ρ   ⁡     (     x   ,   y   ,   t     )         ρ   _         )         ,     
     ⁢       v   ⁡     (     x   ,   y   ,   t     )       =       v   _     ⁡     (     1   -       ρ   ⁡     (     x   ,   y   ,   t     )         ρ   _         )         ,           (   2   )               
where ū and  v  are system parameters determining maximal velocities in the horizontal and vertical directions, and  ρ   is a maximal density in the scene.
 
     The solution to (1) results in a crowd density map ρ, and the velocity fields (u, v) for all (x, y, t) that satisfy given initial and boundary conditions. Although the dimensionality of the differential equations governing the evolution of ρ and (u, v) can be infinite, it is often the case that the flows exhibit low dimensional behavior. 
     A low dimensional state variable at time t is x(t), for which an observable vector y(t)=G(x(t)) corresponds to stacking of the density and velocity fields for all positions x and y at time t. The function G is a mapping from the low dimensional manifold on which x evolves to a space of observables. Then, the solution to (1) determines the transient response and stability of the corresponding dynamical system generally characterized by
 
{dot over ( x )}( t )= F ( x ( t )),  (3)
 
where F(•) is some mapping in the low dimensional manifold on which the dynamical system evolves. For discrete time systems, the dynamical system evolution is characterized by
 
 x   k+1   =F ( x   k ),  (4)
 
where k is a time index.
 
     The Koopman Operator and Dynamic Mode Decomposition 
     The Koopman operator is a linear operator K that satisfies
 
 G ( F ( x   k ))= KG ( x   k )
 
   y   k+1   =Ky   k .  (5)
 
     Although this dynamical system is nonlinear and evolves a finite dimensional manifold, the Koopman operator is linear and infinitly dimensional. Spectral analysis of the Koopman operator can be used to decompose the flow in terms of Koopman modes and associated Koopman eigenvalues that determine the temporal behavior of the corresponding Koopman mode. 
     A dynamic mode decomposition (DMD) can be used to estimate the Koopman modes. The DMD has been used in fluid dynamics as a data-driven and equation-free method for identifying system dynamics. Consider data matrices 
     
       
         
           
             
               
                 
                   
                     
                       Y 
                       1 
                     
                     = 
                     
                       [ 
                       
                         
                           
                             ❘ 
                           
                           
                             ❘ 
                           
                           
                             
                                 
                             
                           
                           
                             ❘ 
                           
                         
                         
                           
                             
                               y 
                               0 
                             
                           
                           
                             
                               y 
                               1 
                             
                           
                           
                             … 
                           
                           
                             
                               y 
                               
                                 m 
                                 - 
                                 1 
                               
                             
                           
                         
                         
                           
                             ❘ 
                           
                           
                             ❘ 
                           
                           
                             
                                 
                             
                           
                           
                             ❘ 
                           
                         
                       
                       ] 
                     
                   
                   ; 
                   
                     
                       Y 
                       2 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             
                               ❘ 
                             
                             
                               ❘ 
                             
                             
                               
                                   
                               
                             
                             
                               ❘ 
                             
                           
                           
                             
                               
                                 y 
                                 1 
                               
                             
                             
                               
                                 y 
                                 2 
                               
                             
                             
                               … 
                             
                             
                               
                                 y 
                                 m 
                               
                             
                           
                           
                             
                               ❘ 
                             
                             
                               ❘ 
                             
                             
                               
                                   
                               
                             
                             
                               ❘ 
                             
                           
                         
                         ] 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     The DMD determine the best fit matrix K that satisfies the relation
 
 Y   2   ≈KY   1 .  (7)
 
     The eigenvectors and eigenvalues of K approximate the Koopman modes and Koopman eigenvalues. Herein, the terms Koopman modes and DMD modes are used interchangeably. 
     SUMMARY OF THE INVENTION 
     The embodiments of the invention provide a method and system for segmenting pedestrian flows in a scene based on motion vectors in an video acquired of a scene by a video camera. The motion vectors are treated as observables of an underlying dynamical system that defines the flows in the scene as organized in an input data matrix. The method determines a low rank Koopman operator that characterizes the dynamical system using the motion vectors. Then, the spectrum of the Koopman operator is analyzed to segment independent flows in the scene. 
     In one embodiment of the invention, a low rank Koopman operator is determined by solving a nuclear norm minimization problem constrained by a one step prediction equation of a dynamic mode decomposition (DMD). The effect of noise can be minimized by selecting an appropriate regularizer for the distribution of the noise. For example, a Frobenius norm or Euclidean norm regularizer can be used for Gaussian distributed noise. A sum of absolute values or one-norm regularizer can be used for sparse impulsive noise. In the case both types of noise exist, two regularizers can be used during the optimization. 
     One embodiment of the invention performs a singular value decomposition (SVD) of the Koopman operator. Row space (left) singular vectors corresponding to the nonzero singular values are then used to cluster blocks in the video as segments corresponding to independent flows. In another embodiment, the eigenvectors of the Koopman operator are determined, and the input data matrix is projected onto the eigenvectors. The eigenvectors that have magnitude projection coefficients greater than a threshold are then used to cluster blocks and segment the independent flows. 
     Another embodiment of the invention characterizes the segmented pedestrian flows by identifying a dynamical regime according to eigenvalues of the eigenvectors that describe the dynamics for each flow. The segmented flow data are projected onto the eigenvectors of the Koopman operator and the eigenvalues that are associated with the eigenvectors that have magnitude projection coefficients greater than a threshold are t used to define the dynamical regime. The dynamical regime characterizes the type of the flow as, for example, stable or unstable. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a schematic of a system and method for segmenting pedestrians flows in a video acquired of a scene according to the embodiments of the invention; 
         FIG. 2  is a flow diagram of a method for determining Koopman modes according to embodiments of the invention; 
         FIG. 3  is a flow diagram of a method for segmenting independent flows using the Koopman modes according to embodiments of the invention; and 
         FIG. 4  is a flow diagram of a method for characterizing independent segmented pedestrian flows. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     As shown in  FIG. 1 , the embodiments of our invention provide a system and method for segmenting pedestrians flows in a video acquired of a scene. The embodiments can use motion vectors extracted from an encoded bitstream, although the motion vectors can be obtained by other means. The motion vectors are treated as observables of an underlying dynamical system composed of several independent dynamics that characterize the flows in the scene. The method determines a low rank Koopman operator that characterizes the dynamical system using the motion vectors. The spectrum of the Koopman operator is then analyzed using singular value decomposition (SVD) or eigenvalue decomposition. The independent crowd flows are determined by clustering the dominant singular vectors or the dominant eigenvectors, depending on the decomposition that is used. 
     System and Method Overview 
     As shown in  FIG. 1 , a video  101  is acquired of pedestrians in a scene  102  with a video camera  103 . The video includes a set of frames. Each frame includes a set of blocks that typically encode motion vectors. Each frame k  104  is processed by the method operating in a processor  100  connected to memory and the camera by input/output buses as known in the art. The video is in a form of an encoded bitstream. 
     As an advantage, an encoded bitstream is usually compressed by at least two orders of magnitude or more. Hence, the bitstream takes much less storage and time to process. In addition, the method can operate at the block level so that the bitstream does not need to be decode to obtain pixels. However, the method can operate in a camera and obtain the motion vectors by other means. 
     The method extracts  110  motion vectors  105  from the video, e.g., blocks in the bitstream, and constructs a data matrix  111  by stacking the motion vectors as described below in greater detail. 
     A low-rank Koopman operator  121  is determined  120  from the data matrix  111 . The spectrum of a low-rank Koopman operator is analyzed to determine dominant Koopman modes  131 . The Koopman modes can than be used to segment  300  blocks in frames of the video to obtain independent flows  109 . 
       FIG. 2  shows the steps  200  used to analyze the spectrum and determine the Koopman modes  131 . The steps include the SVD  210  to obtain singular vectors  211  or an eigenvalue decomposition to determine eigenvectors, or both. Step  240  determines projection coefficients, which are then use to select  250  eigenvectors. The singular vectors or selected eigenvectors, or both, can then be used to determine the Koopman modes  131 . 
       FIG. 3  shows the steps  200  of segmenting the pixels, which takes as input the Koopman modes  131 . Motion vectors are clustered using k-means clustering  310  or graph spectral clustering  320  to obtain the independent flows  109 . 
       FIG. 4  shows the steps  400  of characterizing each segmented flow by determining a dynamical regime of the flow. The dynamical regimes essentially characterizes the flows. The steps take as input the Koopman eigenvectors  221 , the data matrix  111  and the segmented flows  109 . A subset of rows  401  of the data matrix  111  are extracted using the segmented flows  109 . The subset of rows  401  are projected  410  onto the Koopman eigenvectors  221  to determine projection coefficients  411 . An eigenvalue selection step  420  chooses the eigenvalues that correspond to the eigenvectors that have projection coefficients with magnitudes larger than a threshold. The dynamical regime of the flow or flow characteristics  421  is identified by the selected eigenvalues. 
     We now describe the details of the above steps. 
     Motion Vectors and Observable Data Matrix 
     The horizontal and vertical motion vectors of all spatial blocks in a video frame k  104  are u k  and v k , respectively. We construct the data matrix z k =[u k , v k ]  111  by stacking the horizontal and vertical motion vectors in the data matrix  111 . Then, the temporal evolution of the motion vectors can be modeled, using the Koopman operator K, by fitting the data to a one step prediction model
 
 z   k   =Kz   k−1 .  (8)
 
     Notice that (8) ensures that the same operator K models the temporal evolution of the motions vectors u k  and v k . 
     If we have m+1 observations of motion vectors from m+1 video frames, then we can determine the operator K using a dynamic mode decomposition (DMD) by forming the data matrices Z 1  and Z 2  similar to (6) and determining the operator K ls  that achieves the least squares fit, i.e., 
                             K   ls     =       ⁢     arg   ⁢           ⁢       min   K     ⁢       1   2     ⁢              Z   2     -     KZ   1            F   2             ,                 =       ⁢       Z   2     ⁢     Z   1   †         ,                 (   9   )               
where the superscript † indicates the Moore-Penrose pseudo-inverse of a matrix.
 
     If the motion vectors correspond to the true motion of pedestrians in the video, then the operator K ls  characterizes the full dynamics of the flow in the system. However, compressed domain motion vectors are generally very noisy and often contain motion vectors that do not correspond to the true motion. The noise depends on the specific rate-distortion optimization by the video encoder. 
     Low-Rank DMD 
     A noise-free velocity observations y k  measure true motion of pedestrians in the scene. and the motion vectors z k  correspond to noisy observations, such that
 
 z   k   +e   k   =y   k ,  (10)
 
where e k  is additive noise. Let E 1 , E2 and Y 1 , Y 2  be as in (6) for e k  and y k , respectively. Then, the Koopman operator corresponding to the noise-free dynamics is given by a total least squares constraint
 
 Y   2   =KY   1  
 
   Z   2   +E 2 =K ( Z   1   +E   1 ).  (11)
 
     However, the problem (11) is non-convex and ill posed because only Z 1  and Z 2  are observed, and K, E 1  and E 2  are unknowns. To remedy this situation, we invoke the following prior conditions in our problem. Highly complex systems can be accurately modeled by low-order dynamics. This translates to the prior condition that the operator K modeling the evolution of a noise-free system has a low rank. 
     The second prior condition derives from the definition of the matrices E 1  and E 2 . If we denote by I m−1  the identity matrix of size (m−1)×(m−1), and let    0 =[0|I m−1 ] T  and    0 =[I m−1 |0] T  be selection operators that respectively remove the first column and the last column of a matrix with m columns, then E 1  and E 2  satisfy a relation
 
 E   1     0   =E   2     0 .  (12)
 
     We define the low-rank DMD with total least squares constraints (LR-DMDtls) as 
                           min     K   ,     E   1     ,     E   2         ⁢          K        *       +       γ   2     ⁢     (              E   2          F   2     +            e   0          2   2       )         ,     
     ⁢         subject   ⁢           ⁢   to   ⁢           ⁢     Z   2       +     E   ⁢           ⁢   2       =     K   ⁡     (       Z   1     +     E   1       )         ,   and     ⁢     
     ⁢           E   1     ⁢     0       =       E   2     ⁢     0         ,             (   13   )               
where e 0  is the first column in E 1 , γ is a regularization parameter, and ∥K∥ *  is the nuclear norm equal to the sum of the singular values of a matrix K.
 
     If the number of motion blocks in a video frame is N, then the size of the operator K is N×N, which can become difficult to store and determine for high resolution videos. Therefore, we replace K with rank r factors Lε   N×r  and Rε   N×r , such that K=LR T , to reduce the computational complexity. We use a nuclear norm proxy that replaces the nuclear norm of a matrix with an average of the square of the Frobenius norms of its low rank factors 
     
       
         
           
             
               
                 
                   
                     
                        
                       K 
                        
                     
                     * 
                   
                   = 
                   
                     
                       inf 
                       
                         L 
                         , 
                         
                           
                             R 
                             ⁢ 
                             
                               : 
                             
                             ⁢ 
                             K 
                           
                           = 
                           
                             LR 
                             T 
                           
                         
                       
                     
                     ⁢ 
                     
                       1 
                       2 
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             
                                
                               L 
                                
                             
                             F 
                             2 
                           
                           + 
                           
                             
                                
                               R 
                                
                             
                             F 
                             2 
                           
                         
                         ) 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
     Consequently, the low-rank DMD can be formulated as 
                           min     L   ,   R   ,     E   1     ,     E   2         ⁢       1   2     ⁢     (            L        F   2     +          R        F   2       )         +       γ   2     ⁢     (              E   2          F   2     +            e   0          2   2       )         ,     
     ⁢         subject   ⁢           ⁢   to   ⁢           ⁢     Z   2       +     E   2       =       LR   T     ⁡     (       Z   1     +     E   1       )         ,   and     ⁢     
     ⁢           E   1     ⁢     0       =       E   2     ⁢     0         ,             (   15   )               
which we solve using an alternating direction method of multiplier (ADMM).
 
     The problem defined in (15) addresses the case where the noise is Gaussian distributed. If the measurements Z include outlier sparse noise, then we use the following formulation to solve the problem 
                         min     L   ,   R   ,     E   1     ,     E   2     ,     S   1     ,     S   2         ⁢       1   2     ⁢     (            L        F   2     +          R        F   2       )         +       γ   2     ⁢     (              E   2          F   2     +            e   0          2   2       )       +     λ   ⁢            S   2          1         ,     
     ⁢         subject   ⁢           ⁢   to   ⁢           ⁢     Z   2       +     E   2     +     S   2       =       LR   T     ⁡     (       Z   1     +     E   1     +     S   1       )         ,     
     ⁢         E   1     ⁢     0       =       E   2     ⁢     0         ,     
     ⁢         S   1     ⁢     0       =       S   2     ⁢     0         ,           (   16   )               
where the variables S 1  and S 2  represent the outlier noise in the data, and γ and λ are regularization parameters. The solution to (16) can also be obtained using the ADMM procedure.
 
     Flow Segmentation 
     In order to segment  300  different flows having independent dynamics, we evaluate the spectrum of the estimated low-rank Koopman operator K=LR T . Here, we consider different scenarios for performing the segmentation. 
     Segmentation by Singular Value Decomposition 
     In a first approach, we can perform The SVD  210  of the Koopman operator K to extract row space singular vectors U r    211 , such that K=U r Σ r V r   T . The matrix U r  is an N×r matrix, where r is the rank of the operator K. We treat the rows of U r  as feature vectors for every pixel in the video frame to which we apply k-means clustering  310 , as shown in  FIG. 3 , to segment the pixels into a predetermined number of clusters n. Then, the clustered pixels can be used to determine the spatial segmentation of independent flows  109  in the scene  102  as observed in the video  101 . 
     Alternatively, we can apply spectral clustering to the rows of the matrix U r . Here, we normalize the rows of U r  to produce a matrix Ū r  and construct an affinity matrix A=|Ū r Ū r   T | 2 , where the exponent is applied element-wise in this case. We then apply the following steps to determine the principal eigenvectors of the corresponding graph Laplacian.
 
 A=A −diag( A )
 
 D =diag( A· 1)
 
 L   G   =D   −1/2   AD   −1/2 ,  (17)
 
where 1 is the all one vector, the first diag(·) extracts only the diagonal entries of A, and the second diag(A·1) sums the rows of A and places the result on the diagonal of a matrix.
 
     Next, we determine the SVD of L G  to extract the first n row space singular vectors U G . The first singular vector of U G  is all ones so we remove this vector and then normalize the rows of the remaining vectors in U G . 
     Finally, we apply k-means clustering  310  to the rows of U G  to determine the independent flows  109 . While the spectral clustering approach delivers better results than applying k-means clustering directly to the rows of U r , it requires more large to construct a large N×N affinity matrix. 
     Segmentation by Eigenvalue Decomposition 
     In a second approach, we determine an eigenvalue decomposition  220  of K to obtain the eigenvectors U e    221 , and eigenvalues Λ e , such that K=U e Λ e U e   † . Recall, the eigenvectors of K estimate the Koopman modes of the dynamical system. From the data matrix Z 2 , we identify the active modes by projecting Z 2  onto U e  and averaging magnitudes of the projection coefficients in a vector 
             α   =       1   m     ⁢              U   e   †     ⁢     Z   2            ·   1.             
The active modes are the columns of U e  that correspond to the entries in α with a magnitude greater than a threshold.
 
     The matrix containing the active modes is U α . We perform the segmentation by applying k-means clustering to the normalized rows of U α . The pixels identifying with each cluster correspond to an independent flow in the video. 
     Flow Characterization 
     Each segmented flow can be characterized according to a dynamical regime, such that, the type of the flow is labeled as, for example, stable or unstable. The dynamical regime depends on locations of the eigenvalues of K that contribute most to the motion in that flow. To determine the eigenvalues, a subset of rows Z f  of the data matrix associated with a segmented flow f are projected onto the eigenvectors U e  to determine projection coefficients 
               α   f     =       1   m     ⁢              U   e   †     ⁢     Z   f            ·   1.             
The eigenvalues corresponding to the entries in α f , with magnitudes larger than a threshold, are associated with the flow.
 
     The locations of the eigenvalues in the complex plane are used to determine the regime of the flow. For example, stable flows have eigenvalues located on a unit circle, while unstable flows have eigenvalues located outside the unit circle. Moreover, transient flows that gradually decay have eigenvalues located inside the unit circle. Also, the complex phase of the eigenvalues indicate oscillating frequencies of the motion in a flow. 
     Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.